[MATLAB]常微分方程数值求解(ode45 ode15s ode23 solver)

2024-07-11 19:16| 来源: 网络整理| 查看: 265

本次讨论取材于中南大学《科学计算与matlab语言》内容为常微分方程数值解法

常微分方程数值求解的一般概念常微分方程数值求解函数刚性问题常微分方程数值求解的一般概念

求解常微分方程初值问题就是寻找函数y(t)使之满足如下方程：在这里插入图片描述所谓其数值解法，就是求y(t)在离散结点t0处的函数近似值y0的方法

单步法：在计算Yn+1时只用到前一步的Yn，因此在有了初值之后就可以逐步往下计算，其代表是龙格–库塔（Runge-Kutta）法.多步法：在计算Yn+1时，除了用到前一步的值Yn之外,还要用到Yn-p(p=1,2,…,k,k>0)的值，即前面的k步，其代表就是亚当斯（Adams）法。常微分方程的求解函数

函数调用格式为：

[t,y]=solver(filename,tspan,y0,option)

其中，t和y分别给出时间向量和相应的数值解；solver为求常微分方程数值解的函数；filename是定义f(t,y)的函数名，该函数必须返回一个列向量；tspan形式为[t0,tf]，表示求解区间；y0是初始状态向量;Option是可选参数,用于设置求解属性.

常微分方程数值求解函数的统一命名格式：

odennxx

其中ode是Ordinary Differential Equation的缩写，是常微分方程的意思；nn是数字，代表所用方法的阶数；xx是字母，用于标注方法的专门特征。在这里插入图片描述

>> f=@(t,y) (y^2-t-2)/4/(t+1); >> [t,y]=ode23(f,[0,10],2); >> y1=sqrt(t+1)+1; >> plot(t,y,'b:',t,y1,'r'); >>

在这里插入图片描述例2 已知一个二阶线性系统的微分方程为：取a=2，绘制系统的时间响应曲线和相平面图。令x2=x,x1=x’，则得到系统的状态方程：

f=@(t,x) [-2,0;0,1]*[x(2);x(1)]; [t,x]=ode45(f,[0,20],[1,0]); subplot(2,2,1);plot(t,x(:,2)); subplot(2,2,2);plot(x(:,2),x(:,1));

在这里插入图片描述

刚性问题

有一类常微分方程，其解的分量有的变化很快，有的变化很慢，且相差悬殊，这就是所谓的刚性问题(Stiff) 例3 假如点燃一个火柴，火焰球迅速增大直至某个临界体积，然后，维持这一体积不变，其原因是火焰球内部燃烧耗费的氧气和从球表面所或氧气达到平衡。在这里插入图片描述主要观看λ值得大小，是否会真正影响问题模型！

>> lamda=0.01; >> f=@(t,y) y^2-y^3; >> tic;[t,y]=ode45(f,[0,2/lamda],lamda);toc Elapsed time is 0.034000 seconds. >> disp(['ode45计算的点数' num2str(length(t))]); ode45计算的点数157 >>

时间已过0.034000秒 ode45计算的点数157

>> lamda=1e-5; >> f=@(t,y) y^2-y^3; >> tic;[t,y]=ode45(f,[0,2/lamda],lamda);toc Elapsed time is 5.690000 seconds. >> disp(['ode45计算的点数' num2str(length(t))]); ode45计算的点数120565 >>

不同的λ值，式子所用的秒数不同，问题瞬间变成刚性，直接换用ode15s，再去计算！直接秒出答案！

>> tic;[t,y]=ode15s(f,[0,2/lamda],lamda);toc Elapsed time is 0.348000 seconds. >> disp(['ode15s计算的点数' num2str(length(t))]); ode15s计算的点数98 >>

在这里插入图片描述

>> f=@(t,x) [-8/3*x(1)+x(2)*x(3);-10*x(2)+10*x(3);-x(2)*x(1)+28*x(2)-x(3)]; >> [t,x]=ode45(f,[0,50],[0,0,10^(-10)]); >> plot(t,x)

【本文地址】

公司简介

联系我们