1、Matlab在常微分方程在常微分方程求解中的应用求解中的应用实验目的实验目的 （1 1）学会用学会用MatlabMatlab软件求解微分方程的初值问题软件求解微分方程的初值问题 （2 2）了解微分方程数值解思想，掌握基本的微分方程数值解方法）了解微分方程数值解思想，掌握基本的微分方程数值解方法 （3 3）学会根据实际问题建立简单微分方程数学模型）学会根据实际问题建立简单微分方程数学模型 （4 4）了解计算机数据仿真、数据模拟的基本方法）了解计算机数据仿真、数据模拟的基本方法 n17世纪：初等解法世纪：初等解法n18世纪：初等解法和无穷级数方法世纪：初等解法和无穷级数方法n19世纪：解的存在性、
2、奇点理论、定性世纪：解的存在性、奇点理论、定性理论、稳定性理论理论、稳定性理论l 包含一个自变量和它的未知函数以及未知函数的导数的包含一个自变量和它的未知函数以及未知函数的导数的等式等式l 形成和发展与力学、天文学、物理学及其他自然科学技形成和发展与力学、天文学、物理学及其他自然科学技术的发展互相促进和推动术的发展互相促进和推动 常微分方程常微分方程yxf,ybxaD,:1,yx2,yx2121,yyLyxfyxf xy定理 设函数设函数在区域在区域上连续，且在区域上连续，且在区域D内满足李普希兹内满足李普希兹(Lipschitz)条件，即存在条件，即存在正数正数L，使得对于，使得对于R内任意
3、两点内任意两点与与,恒有恒有则初值问题则初值问题(1)的解的解存在并且唯一。存在并且唯一。常微分方程的解析解常微分方程的解析解 求微分方程（组）的解析解命令:dsolve(方程方程1,方程方程2,方程方程n,初始条件初始条件,自变量自变量)结 果：u=tg(t-c)记号：在表达微分方程时，用字母D表示求微分，D2、D3等表示求高阶微分。D后所跟字母为因变量，自变量可以指定或由系统规则选定为缺省。例如：微分方程 可以表示为D2y=0.022dxyd 解解 输入命令:y=dsolve(D2y+4*Dy+29*y=0,y(0)=0,Dy(0)=15,x)结 果 为:y=3e-2xsin（5x）例例
4、3 求微分方程组的通解.zyxdtdzzyxdtdyzyxdtdx244354332解解 输入命令：x,y,z=dsolve(Dx=2*x-3*y+3*z,Dy=4*x-5*y+3*z,Dz=4*x-4*y+2*z,t)；x=simple(x)%将将x化简化简 y=simple(y)z=simple(z)结 果 为：x=(c1-c2+c3+c2e-3t-c3e-3t)e2t y=(c1e-4t+c2e-4t+c2e-3t-c3e-3t+c1-c2+c3)e2t z=(-c1e-4t+c2e-4t+c1-c2+c3)e2t 虽然说解析解是最精确的，但是实际问虽然说解析解是最精确的，但是实际问题中
5、常要求研究常微分方程的数值解题中常要求研究常微分方程的数值解常微分方程的数值解常微分方程的数值解 常微分方程中只有一些典型方程能求出初等解(用初等函数表示的解)。另外，有些初值问题虽然有初等解，但由于形式太复杂不便于应用。因此，有必要探讨常微分方程初值问题的数值解法。以下主要介绍一阶常微分方程初值问题几种经典数值解方法：欧拉法及改进的欧拉法。其它方法：龙格-库塔法、阿达姆斯方法；一阶微分方程组与高阶方程初值问题的数值解法；二阶常微分方程值问题的差分方法等。求解常微分方程初值问题的数值解的整体思路：00)(),(yxyyxfdxdy （1）)(xynxxxx210,210nyyyy ny1iii
6、xxh寻求准确解在一系列离散节点：上的近似值 称为问题的数值解，数值解所满足的离散方程统称为差称为步长，实用中常取定步长。分格式，建立数值解法的一些途径建立数值解法的一些途径001i)y(xy)f(x,y ,1,2,1,0 ,xynihxi解微分方程：可用以下离散化方法求设1、用差商代替导数、用差商代替导数 若步长h较小，则有hxyhxyxy)()()(故有公式：1-n,0,1,2,i )(),(001xyyyxhfyyiiii此即欧拉法欧拉法Euler法的几何意义:2、使用数值积分、使用数值积分对方程y=f(x,y),两边由xi到xi+1积分，并利用梯形公式，有：)(,()(,(2)(,()
7、()(11111iiiiiixxiixyxfxyxfxxdttytfxyxyii实际应用时，与欧拉公式结合使用：,2,1,0 ),(),(2),()(11)1(1)0(1kyxfyxfhyyyxhfyykiiiiikiiiii的计算。然后继续下一步，取时，当满足，对于已给的精确度）（y y 2i111i)(1)1(1kikikiyyy此即改进的欧拉法改进的欧拉法故有公式：)(),(),(200111xyyyxfyxfhyyiiiiii使得中值定理可知，存在可微的条件下，由微分在函数,10)(xy)()()(1hxyhxyxynnn上的平均斜率在区间为函数称,)()()(11nnnnxxxyhx
8、yxyhxyxyxxyxfxyxfxxxxyxnnnnnnnnnnnnn)()(,),(),(,)(111-111上的平均斜率区间的平均值代替了曲线在处的预测斜率和在处的斜率线在改进的欧拉公式：用曲上的平均斜率在区间了函数的斜率（导数值）代替欧拉公式：用曲线在的精度更高的计算公式计算由此可构造出由上的平均斜率，在区间的加权平均值作为曲线后将这些点处的斜率值值，然内多预测几个点的斜率间的公式；如果设法在区一种计算算法，就可以得到上的平均斜率提供一种间启发：只要对曲线在区11111,nnnnnnnnnyyxxxxyxx3、使用泰勒公式、使用泰勒公式 以此方法为基础，有龙格龙格-库塔法库塔法、线性多
9、步法线性多步法等方法。4、数值公式的精度、数值公式的精度 当一个数值公式的截断误差可表示为O（hk+1）时（k为正整数，h为步长），称它是一个k阶公式阶公式。k越大，则数值公式的精度越高。欧拉法是一阶公式，改进的欧拉法是二阶公式。龙格-库塔法有二阶公式和四阶公式。线性多步法有四阶阿达姆斯外插公式和内插公式。ODE 指令列表 nMATLAB 用於求解常微分方程式的指令:指令方法应用ODE类别ode45Explicit Runge-Kutta(4,5)pair of Dormand-PrinceNonstiff ODEode23Explicit Runge-Kutta(2,3)pair of Bo
10、gacki and ShampineNonstiff ODEode113Variable order Adams-Bashforth-Moulton PECE solverNonstiff ODEode15sNumerical differentiation formulas(NDFS)Stiff ODEode23sModified Rosenbrock formula of order 2Stiff ODEode23t Trapezoidal rule with a“free”interpolantStiff ODEode23tbImplicit Runge-Kutta formula wi
11、th a backward differentiation formula of order twoStiff ODE 适用于 Nonstiff 系統 n速率（即微分值）差异相常大 n使用一般的ode45、ode23或ode113來求解,可能会使得积分的步长(Step Sizes)变得很小,以便降低积分误差至可容忍范围以內,会导致计算时间过长 n专门对付Stiff系统的指令，例如ode15s、ode23s、ode23t及ode23tb 适用于 Nonstiff 系統 一般的常微分方程式都是 Nonstiff 系統 直接采用 ode45、ode23 或 ode113 來求解 提示 n使用 Sim
12、ulink 來求解常微分方程式Simulink是和MATLAB共同使用的一套软件 可使用拖拉的方式來建立动力系统 可直接产生C语言代码或进行动画演示 功能非常强大 n使用ODE指令时，必须先将要求解的ODE表示 成一个函数 输入为t(时间)及y(状态函数:State Variables）输出则为dy（状态变量的微分值）n若ODE函数的文件为odeFile.m,则调用ODE指令 的格式如下：t,y=solver(odeFile,t0,t1,y0)t0,t1 是积分的时间区间 y0代表起始条件（Initial Conditions）solver是前表所列的各种ODE指令 t是输出的时间向量 y是对
13、应的状态变量向量 t,y=solver（f,ts,y0,options）ode45 ode23 ode113ode15sode23s由待解方程写成的m-文件名ts=t0,tf，t0、tf为自变量的初值和终值函数的初值ode23：组合的2/3阶龙格-库塔-芬尔格算法ode45：运用组合的4/5阶龙格-库塔-芬尔格算法自变量值函数值用于设定误差限(缺省时设定相对误差10-3,绝对误差10-6),命令为：options=odeset（reltol,rt,abstol,at）,rt，at：分别为设定的相对误差和绝对误差.1、在解n个未知函数的方程组时，y0和y均为n维向量，m-文件中的待解方程组应以y
14、的分量形式写成.2、使用Matlab软件求数值解时，高阶微分方程必须等价地变换成一阶微分方程组.注意注意:Van der Pol微分方程n1928年荷兰的范德波耳(Van der Pol)为描述LC回路的电子管振荡器建立了著名的vanderPol方程.它在自激振荡理论中有着重要的意义,一直作为数学物理方程中的一个基本方程.这是一个具有可变非线性阻尼的微分方程,代表了一类极为典型的非线性问题.和其他非线性微分方程在数学上无法精确求解一样,人们一直在努力寻找求解这类方程近似解析解的方法,并乐于用Van der Pol方程来检验求解方法的有效性.Van der Pol微分方程其方程形式为：n令 0)
15、1(2yyyy121yyyy将Van der Pol微分方程化成标准形式 1221221)1(yyyyyy),(tyFyy写成向量的形式：为一个向量，代表状态变量n设 =1，ODE文件(vdp1.m)显示如下：type vdp1.m type vdp1.m function dy=vdp1(t,y)mu=1;dy=y(2);mu*(1-y(1)2)*y(2)-y(1);建立了vdp1.m后，即可选用前述ODE指令来求解 n在 =1时，van der Pol方程并非 Stiff系统，所以使用ode45来显示结果 odeBasic01.mode45(vdp1,0 25,3 3);n 0 25代表积
16、分的时间区间，3 3 则代表起始条件n因为没有输出变量，所以上述程序执行结束后，MATLAB 只画出状态变量对时间的图形0510152025-3-2-101234 odeGetData01.mt,y=ode45(vdp1,0 25,3 3);plot(t,y(:,1),t,y(:,2),:);xlabel(Time t);ylabel(Solution y(t)and y(t);legend(y(t),y(t);0510152025-3-2-101234Time tSolution y(t)and y(t)y(t)y(t)n也可以画出 及 在相位平面（Phase Plane）的运动情况 ode
17、PhastPlot01.m)(ty)(tyt,y=ode45(vdp1,0 25,3 3);plot(y(:,1),y(:,2),-o);xlabel(y(t);ylabel(y(t);-3-2-101234-3-2-10123y(t)y(t)n当 值越来越大时，Van der Pol方程就渐变为一个Stiff系统，此时若要求解该类问题，就必须用对应于Stiff系统的指令 n将 值改成 1000，ODE 文件改写成（vdp2.m）：type vdp2.m function dy=vdp2(t,y)mu=1000;dy=y(2);mu*(1-y(1)2)*y(2)-y(1);n选用专门用来解决
18、Stiff 系统问题的 ODE 指令，例如 ode15s，來求解此系统并作图显示 ode15s01.mt,y=ode15s(vdp2,0 3000,2 1);subplot(2,1,1);plot(t,y(:,1),-o);xlabel(Time t);ylabel(y(t);subplot(2,1,2);plot(t,y(:,2),-o);xlabel(Time t);ylabel(y(t);n由上图可知，的变化相当剧烈（超过 ），这也就是 Stiff 系统的特点 050010001500200025003000-4-2024Time ty(t)05001000150020002500300
19、0-2000-1000010002000Time ty(t)(ty1000n绘制二维平面相位图：ode15s02.mn若要产生在某些特定时间点的状态变量值，则调用 ODE 指令的格式可改成：t,y=solver(odeFile,t0,t1,tn,y0)其中t0,t1,tn即是特定时间点所形成的向量 t,y=ode15s(vdp2,0 3000,2 1);subplot(1,1,1);plot(y(:,1),y(:,2),-o);xlabel(y(t);ylabel(y(t)-2.5-2-1.5-1-0.500.511.522.5-1500-1000-500050010001500y(t)y(t
20、)ODE 指令的选项nODE 指令可以接受第四个输入变量，代表积分过程用到的各种选项（Options），此种ODE指令的格式为：t,y=solver(odeFile,t0,tn,y0,options)其中 options 是由 odeset 指令來控制的结构变量 结构变量即包含了积分过程用到的各种选项 odeset的一般格式如下：options=odeset(name1,value1,name2,value2,)nname1的值为value1，name2的值为value2 n也可以只改变一个现有的options中某栏位的值，格式如下：newOptions=odeset(options,name
21、,value);n若要获取某栏位的值，可用 odeget，格式如下：value=odeget(otpions,name);n当使用 odeset 指令時，若无输入任何变量名，则 odeset 会显示所有变量名及其值，并以大括号代表预设值 odeset AbsTol:positive scalar or vector 1e-6 RelTol:positive scalar 1e-3 NormControl:on|off NonNegative:vector of integers OutputFcn:function_handle OutputSel:vector of integers Ref
22、ine:positive integer Stats:on|off InitialStep:positive scalar MaxStep:positive scalar BDF:on|off MaxOrder:1|2|3|4|5 Jacobian:matrix|function_handle JPattern:sparse matrix Vectorized:on|off Mass:matrix|function_handle MStateDependence:none|weak|strong MvPattern:sparse matrix MassSingular:yes|no|maybe
23、 InitialSlope:vector Events:function_handle 由 odeset 产生的 ODE 选项 310610类别栏位名称资料形态预设值说明誤差容忍度之相关栏位RelTol正数量相对误差容忍度AbsTo1正数量或向量绝对误差容忍度积分输出之相关栏位OutPutFcn字串odeplot输出函式（若 ODE 指令无输出变量，則在数值积分执行完毕后，MATLAB会呼叫此输出函数）OutputSel索引向量全部ODE 指令之输出变量的索引值，以決定那些输出变量之元素将被送到输出函式Refine正整数1或4(for ode45)Refine=2 可产生两倍数量的输出点，Re
24、fine=3 可产生三倍数量的输出点，依此类推Statson 或 offoffStats=on 产生计算过程的各种统计资料常用到的档位说明 nf在积分误差容忍度方面，每一次积分所产生的局部误差e(i)，必须满足下列方程式：max(RelTol*,AbsTol(i)其中i代表第i个状态 降低RelTol及AbsTol来求得更精确的积分结果)(ie)(iyodeRelTol01.msubplot(2,1,1);ode45(vdp1,0 25,3 3);title(RelTol=0.01);options=odeset(RelTol,0.00001);subplot(2,1,2);ode45(vdp
25、1,0 25,3 3,options);title(RelTol=0.0001);n第一个图所使用的相对误差值是0.01（预设值），第二个图所使用的相对误差值是0.00001，因此我们得到较细密的点，但所花的计算时间也会比较长 0510152025-4-2024RelTol=0.010510152025-4-2024RelTol=0.0001积分输出说明n积分输出的相关处理 选用一个 OutputFcn OutputFcn 的预设值是“odeplot”，其功能为画出所有的状态变量的图像 其它可用的函式 nodephas2：画出 2-D 的相位平面（Phase Plane）nodephas3：画出 3-D 的相位平面nodeprint：随时在指令视窗显示计算结果options=odeset(OutputFcn,odephas2);%使用 odephas2 进行绘图ode45(vdp1,0 25,3 3,options);-3-2-101234-3-2-10123
