常微分方程初值问题的matlab解法

《常微分方程初值问题的matlab解法》由会员分享，可在线阅读，更多相关《常微分方程初值问题的matlab解法（10页珍藏版）》请在金锄头文库上搜索。

1、1用 Matlab 求常微分方程(ODE) 的初值问题(IVP)本节考虑一阶常微分方程* MERGEFORMAT (1.1)00(,) ufttT的数值求解问题，包括算法公式及编程问题。对一阶常微分方程组的初值问题* MERGEFORMAT (1.2)0111122 220012 ()(,) (,)()mmmmutuftutTft tLLL可以通过引入列向量 化成类似* MERGEFORMAT (1.1)的形式0,ufr* MERGEFORMAT (1.3)00(,)ufttTr其中* MERGEFORMAT 11012220012()()(,),)()()(,)mmmmmttftuuut f

2、tttftLrrrMM(1.4)另外一个高阶常微分方程的初值问题* MERGEFORMAT (1.5)() (1)0(1)00, ()mmuftutTu L也可以通过变换 化成维微分方程组：1123,mu* MERGEFORMAT (1.6)2312(,)mmuuftuLL我们在设计算法时通常先对一维一阶常微分方程* MERGEFORMAT (1.1)进行，然后再将这个算法写成适合求解* MERGEFORMAT (1.3)的向量形式，并以向量形式来进行编程。1、 非刚性系统与刚性系统 当 时微分方程组* MERGEFORMAT (1.3)变成(,)()ftuAtrr* MERGEFORMAT

3、(1.7)()uAtrr如果系数矩阵 特征值 满足 ，对应的特征向量为 ,则* mR12,mL12mL12,mvrrL2MERGEFORMAT (1.7)具有通解* MERGEFORMAT (1.8)1()()jmtjutcevrr其中 为* MERGEFORMAT (1.7) 的一个特解。如果* MERGEFORMAT (1.7)中的矩阵满足()tr* MERGEFORMAT (1.9)j1()Re0,22 ()e)1 msL就称微分方程组* MERGEFORMAT (1.7) 是刚性的或坏条件的或病态的。s 称为刚性比。对于一般的一阶常微分方程组* MERGEFORMAT (1.3)，如果

4、多元向量值函数 对* (,)fturMERGEFORMAT (1.4)中向量 的分量的 Jacobi 矩阵ur* MERGEFORMAT (1.10)1122112()()mmmutffuuJtffurLMOL的特征值 对于求解区间上的任何 满足* MERGEFORMAT (1.9)式，则称微分方程12(),()mttL0,tT组* MERGEFORMAT (1.3)为刚性的。2、解法器及调用格式解法器 适用类型 何时使用 使用算法ode45 Nonstiff Most of the time. This should be the first solver you try.Runge-Kut

5、ta (4,5) 格式ode23 Nonstiff If using crude error tolerances or solving moderately stiff problems.Runge-Kutta (2,3) 格式ode113 Nonstiff If using stringent error tolerances or solving a computationally intensive ODE file.Adams-Bashforth-Moulton PECE 格式ode15s Stiff If ode45 is slow because the problem is s

6、tiff.ode23s Stiff If using crude error tolerances to solve stiff systems and the mass matrix is constant.ode23t Moderately Stiff If the problem is only moderately stiff and you need a solution without numerical damping.ode23tb Stiff If using crude error tolerances to solve stiff 3systems.有如下三种调用方法，其

7、中 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:是一个函数句柄，指向初值问题* MERGEFORMAT (1.3)中的函数 ，它的基本形式是 dydt (,)ftur= f(t,y)， 其中 dydt,y 是列向量，而 t 是标量。(2) tspan：是一个向量，用于指定求解区间 。要求这个向量的分量是有序的，或者单

8、调上tspan(1),ted)升或者单调下降。(3) y0：是初值问题* MERGEFORMAT (1.3)中的初始列向量 。0ur(4) options：用于对解法器缺省的求解参数进行设置。如果不想改变缺省值，可以用空矩阵 来代替。MATLAB 的 ODE 解法器设计了一个结构变量用于设置解法器的各项公共参数，比如精度，步长等。其中已经定义了一些缺省值，用无参数的 odeset 命令可以列出所有的公共参数及其缺省值：AbsTol: positive scalar or vector 1e-6 绝对误差RelTol: positive scalar 1e-3 相对误差NormControl:

9、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 对刚性方程组

10、的解法器可以提供相应的 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 每个解法器可以设置的参数如下表:Parameters ode45 ode23 ode113 ode15s ode23s ode23t ode23tbR

11、elTol, AbsTol, NormControlOutputFcn, OutputSel, Refine, StatsEvents4MaxStep, InitialStepJacobian, JPattern, Vectorized - - -MassMStateDependence MvPatternMassSingular- -InitialSlope - - - - -MaxOrder, BDF - - - - - -(5) 输入参数 p1,p2. 等等是可选的参数,并且由解法器传递给函数 odefun.输出参数的含义是：(1) t：节点列向量。(2) y：在节点处的解的近似值数组，

12、第 i 行对应于 t 的第 i 行那个节点处的近似解。即输出是如下形式：节点 t y1(t) y2(t) ym(t)例题 1、分别用步长 来求初值问题在区间 上的近似解，并与准确解0.5,.,0h0,8进行图形上的比较。2ut* MERGEFORMAT (1.11)2 (0)1tu解：利用如下代码，注意其中作为输入参数的函数的代码的编写方法，以及 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)

13、;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。从图形中可以明显看出，当时间变量较大时比如 所求的近似5t解严重偏离了真解。这说明初值问题* MERGEFORMAT (1.11)在时刻 0 有较好的稳定性，但是在时刻 4 以后可能稳定性就不太好了，为些我们利用*

14、MERGEFORMAT (1.11)构造一个初值问题：* MERGEFORMAT (1.12)241()3tu5它与初值问题* MERGEFORMAT (1.11)有相同的解。将上述代码适当改动来求解初值问题* MERGEFORMAT (1.12)，并画出图形如图 2，发现还是在初始时刻 4 附近有较好的稳定性，所以我们在求解时不能希望在很大的区间上求解。图 1 图 2 例 2、求解二阶 Van der Pol 微分方程所构成的初值问题在区间 上的近似解。其中参数 是一个正数，0,当较大时这个方程是刚性的。* MERGEFORMAT (1.13)2(1) 20,uut解：这是一个高阶方程，要用

15、变量代换化成向量形式的一阶方程组形式* MERGEFORMAT (1.3)，令 则vu化成方程组* MERGEFORMAT (1.14)2(1)0 20,uvut再令 ， 则* MERGEFORMAT (1.14)的一阶向量形式为* (),()Tutvtr 2(,),()TftuvrMERGEFORMAT (1.3)。下列代码中编写了向量值函数 ，由还有一个可选参数 故要注意解法器的调用(,)ftur 形式。输出矩阵 u 的每一行代表解向量在每个节点处的近似函数值。对于本题来说输出矩阵 u 的第二列相当于方程* MERGEFORMAT (1.13)解函数 导数的近似值。可以对 进行试验解法器()ut 0.5,12,01ode45 的求解结果及图像。function zzz %ode13.mmu=10