常微分方程初值问题的MATLAB解法.doc

用Matlab求常微分方程(ODE)的初值问题(IVP)本节考虑一阶常微分方程的数值求解问题，包括算法公式及编程问题。对一阶常微分方程组的初值问题可以通过引入列向量化成类似的形式其中另外一个高阶常微分方程的初值问题也可以通过变换化成维微分方程组：我们在设计算法时通常先对一维一阶常微分方程进行，然后再将这个算法写成适合求解的向量形式，并以向量形式来进行编程。1、 非刚性系统与刚性系统 当时微分方程组变成如果系数矩阵特征值满足，对应的特征向量为,则具有通解其中为的一个特解。如果中的矩阵满足就称微分方程组是刚性的或坏条件的或病态的。s称为刚性比。对于一般的一阶常微分方程组，如果多元向量值函数对中向量的分量的Jacobi矩阵的特征值对于求解区间上的任何满足式，则称微分方程组为刚性的。2、解法器及调用格式解法器适用类型何时使用使用算法ode45NonstiffMost of the time. This should be the first solver you try.Runge-Kutta (4,5) 格式ode23NonstiffIf using crude error tolerances or solving moderately stiff problems.Runge-Kutta (2,3) 格式ode113NonstiffIf using stringent error tolerances or solving a computationally intensive ODE file.Adams-Bashforth-Moulton PECE 格式ode15sStiffIf ode45 is slow because the problem is stiff.ode23sStiffIf using crude error tolerances to solve stiff systems and the mass matrix is constant.ode23tModerately StiffIf the problem is only moderately stiff and you need a solution without numerical damping.ode23tbStiffIf using crude error tolerances to solve stiff systems.有如下三种调用方法，其中solver是上面三个解法器的名子。(1) t,y = solver(odefun,tspan,y0)；(2) t,y = solver(odefun,tspan,y0,options)；(3) t,y = solver(odefun,tspan,y0,options,p1,p2.)；它们输入参数是：(1) odefun:是一个函数句柄，指向初值问题中的函数，它的基本形式是dydt = f(t,y)，其中dydt,y是列向量，而t是标量。(2) tspan：是一个向量，用于指定求解区间。要求这个向量的分量是有序的，或者单调上升或者单调下降。(3) y0：是初值问题中的初始列向量。(4) options：用于对解法器缺省的求解参数进行设置。如果不想改变缺省值，可以用空矩阵来代替。MATLAB的ODE解法器设计了一个结构变量用于设置解法器的各项公共参数，比如精度，步长等。其中已经定义了一些缺省值，用无参数的odeset命令可以列出所有的公共参数及其缺省值： AbsTol: positive scalar or vector 1e-6 绝对误差 RelTol: positive scalar 1e-3 相对误差 NormControl: on | off 是否用向量范数控制误差，缺省是按每个分量进行误差控制。 OutputFcn: function OutputSel: vector of integers Refine: positive integer Stats: on | off InitialStep: positive scalar 初始步长 MaxStep: positive scalar 最大步长，缺省为求解区间的十分之一。 BDF: on | off 是否利用向后差分来求Jacobi矩阵。 MaxOrder: 1 | 2 | 3 | 4 | 5 Jacobian: matrix | function 对刚性方程组的解法器可以提供相应的Jacobi矩阵信息 JPattern: sparse matrix Vectorized: on | off Mass: matrix | function MStateDependence: none | weak | strong MvPattern: sparse matrix MassSingular: yes | no | maybe InitialSlope: vector Events: function 每个解法器可以设置的参数如下表:Parametersode45ode23ode113ode15sode23sode23tode23tbRelTol, AbsTol, NormControlOutputFcn, OutputSel, Refine, StatsEventsMaxStep, InitialStepJacobian, JPattern, Vectorized-MassMStateDependence MvPatternMassSingular-InitialSlope-MaxOrder, BDF-(5) 输入参数p1,p2. 等等是可选的参数,并且由解法器传递给函数odefun.输出参数的含义是：(1) t：节点列向量。(2) y：在节点处的解的近似值数组，第i行对应于t的第i行那个节点处的近似解。即输出是如下形式：节点ty1(t)y2(t)ym(t)例题1、分别用步长来求初值问题在区间上的近似解，并与准确解进行图形上的比较。解：利用如下代码，注意其中作为输入参数的函数的代码的编写方法，以及ODE解法器ode45的调用方法。function zzzh=0.1;y0=1;tspan=0:h:8;t,y = ode45(f,tspan,y0);t,yhold on;plot(tspan,y,-k,tspan,fu(tspan),-b);title(Solution of IVP in eg1);xlabel(time t);ylabel(solution u);legend(numerical solutons,true solutions);function dudt=f(t,u)dudt=u-2*t./u;function u=fu(t)u=sqrt(1+2*t);由于结果较多我们只画出它们的图形如图1。从图形中可以明显看出，当时间变量较大时比如所求的近似解严重偏离了真解。这说明初值问题在时刻0有较好的稳定性，但是在时刻4以后可能稳定性就不太好了，为些我们利用构造一个初值问题：它与初值问题有相同的解。将上述代码适当改动来求解初值问题，并画出图形如图2，发现还是在初始时刻4附近有较好的稳定性，所以我们在求解时不能希望在很大的区间上求解。 图1 图2 例2、求解二阶Van der Pol 微分方程所构成的初值问题在区间上的近似解。其中参数是一个正数，当较大时这个方程是刚性的。解：这是一个高阶方程，要用变量代换化成向量形式的一阶方程组形式，令则化成方程组再令，则的一阶向量形式为。下列代码中编写了向量值函数，由还有一个可选参数故要注意解法器的调用形式。输出矩阵u的每一行代表解向量在每个节点处的近似函数值。对于本题来说输出矩阵u的第二列相当于方程解函数导数的近似值。可以对进行试验解法器ode45的求解结果及图像。function zzz %ode13.mmu=10;mustr=num2str(mu);h=0.1;u0=2,0;tspan=0:h:20;t,u = ode45(f,tspan,u0,mu);t,ufigure;hold on;plot(t,u(:,1),-k,t,u(:,2),-b);title(strcat(例2中的初值问题 mu=,mustr);xlabel(time t);ylabel(solution u);legend(数值解,导函数数值解);hold off;function dudt=f(t,u,mu)dudt=u(2),mu*(1-u(1)2)*u(2)-u(1);表1：当mu=0.5时的几个解节点tu(t)u(t)02.000000.10001.9905-0.18550.20001.9638-0.34450.30001.9223-0.48140.40001.8681-0.60080.50001.8026-0.70680.60001.7271-0.80340.70001.6422-0.89400.80001.5484-0.98130.90001.4459-1.06771.00001.3349-1.15461.10001.2150-1.24431.20001.0858-1.3386 如果再增大，则ode45解法将无法求出较好的数值解，原因是当较大时方程组的Jacobi矩阵模最大特征值与模最小特征比值比较大，或者说刚性比大。在时刻0处，当时为：，其特征值为：-0.0033，-299.9967，刚性比约为90909，非常大，这种方程不能用显式的四阶RK方法来求解，必须用某种隐式方法才能求解。下面代码是利用ode23s解法器来求解Van der Pol方程当时的解，由于刚性比较大，对应绝对值较小特征值的解分是变化较缓慢，故应当采用很大的求解区间，才能观察到解的变化情况，因此取求解区间为。同学们自己也可以使用较小的区间来试验，看看数值解的模拟图像能否得到真解的性态。由于数据比较多，故只给出数值解模拟图像。 function zzzmu=100;mustr=num2str(mu);h=0.1;u0=2,0;tspan=0:h:5*mu;t,u = ode15s(f,tspan,u0,mu);t,u;figure;hold on;plot(t,u(:,1),-b);title(strcat(例2中的初值问题 mu=,mustr);xlabel(time t);ylabel(solution u);legend(数值解);hold off;function dudt=f(t,u,mu)dudt=u(2),mu*(1-u(1)2)*u(2)-u(1);例3、求解如下非刚性初值问题，它可以描述刚性物体不受外力作用时的运动状态。解：由于是非刚性问题故可以使用ode45来求解，这个例子说明怎样设置解法器的options选项中的参数。具体代码如下：function zzzoptions = odeset(RelTol,1e-4,AbsTol,1e-4 1e-4 1e-5);tspan=0,12;y0=0,1,1;T,Y = ode45(rigid,0 12,0 1 1,options);T,Yplot(T,Y(:,1),-,T,Y(:,2),-.,T,Y(:,3),.)title(例题3中的图);xlabel(时间 t);ylabel(解 y(t);legend(数值解y_1(t),数值解y_2(t),数值解y_3(t);function dy = rigid(t,y)dy = zeros(3,1); % a column vectordy(1) = y(2) * y(3);dy(2) = -y(1) * y(3);dy(3) = -0.51 * y(1) * y(2); 解的数值模拟图像为：练习1、求解 例1中的van der Pol 微分方程，其中，并到互联网上去搜索这个微分方程是可以描述哪些问题？练习2、 用上面的所有方法求解如下初值问题的解（1）（2）（3）（4）练习3、 求解 例1中的van der Pol 微分方程，并到互联网上去搜索这个微分方程是可以描述哪些问题？其中：练习4、 Lorenz吸引子是在浑沌世界中经常用到的常微分方程组，它实际上描述了容器内的液体被加热后在容器内的流动情况，可以用于描述大气的对流情况。其中均为常数。试求解如下吸引子问题，并作出数值模拟图像。（1）取及初值（2）取，自已确定不同的初值进行试验。练习5、 Rssler吸引子问题。1976年瑞士数学家Otto Rssler发现了另一个具有吸引子的常微分方程他采用了参数以及，他想找到一种刻画浑沌世界的最小微分方程，但是他导出的这给微分方程比Lorenz的微分方