《Matlab在常微分方程求解中应用课件.ppt》由会员分享，可在线阅读，更多相关《Matlab在常微分方程求解中应用课件.ppt（50页珍藏版）》请在三一办公上搜索。
1、Matlab在常微分方程求解中的应用,实验目的（1）学会用Matlab软件求解微分方程的初值问题（2）了解微分方程数值解思想，掌握基本的微分方程数值解方法（3）学会根据实际问题建立简单微分方程数学模型（4）了解计算机数据仿真、数据模拟的基本方法,17世纪：初等解法18世纪：初等解法和无穷级数方法19世纪：解的存在性、奇点理论、定性理论、稳定性理论,包含一个自变量和它的未知函数以及未知函数的导数的等式,形成和发展与力学、天文学、物理学及其他自然科学技术的发展互相促进和推动,常微分方程,定理 设函数,在区域,上连续，且在区域D内满足李普希兹(Lipschitz)条件，即存在正数L，使得对于R内任意
2、两点,与,恒有,则初值问题(1)的解,存在并且唯一。,常微分方程的解析解,结 果：u=tg(t-c),记号：在表达微分方程时，用字母D表示求微分，D2、D3等表示求高阶微分。D后所跟字母为因变量，自变量可以指定或由系统规则选定为缺省。例如：微分方程 可以表示为D2y=0.,解 输入命令:y=dsolve(D2y+4*Dy+29*y=0,y(0)=0,Dy(0)=15,x),结 果 为:y=3e-2xsin（5x）,解 输入命令：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=simp
3、le(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,虽然说解析解是最精确的，但是实际问题中常要求研究常微分方程的数值解,常微分方程的数值解,常微分方程中只有一些典型方程能求出初等解(用初等函数表示的解)。另外，有些初值问题虽然有初等解，但由于形式太复杂不便于应用。因此，有必要探讨常微分方程初值问题的数值解法。以下主要介绍一阶常微分方程初值问题几种经典数值解方法：欧拉法及改进的欧拉法。其它方法：
4、龙格-库塔法、阿达姆斯方法；一阶微分方程组与高阶方程初值问题的数值解法；二阶常微分方程值问题的差分方法等。,求解常微分方程初值问题的数值解的整体思路：,（1）,寻求准确解,在一系列离散节点：,上的近似值,称为问题的数值解，数值解所满足的离散方程统称为差,称为步长，实用中常取定步长。,分格式，,建立数值解法的一些途径,1、用差商代替导数,若步长h较小，则有,故有公式：,此即欧拉法,Euler法的几何意义:,找到了积分的一条近似折线,2、使用数值积分,对方程y=f(x,y),两边由xi到xi+1积分，并利用梯形公式，有：,实际应用时，与欧拉公式结合使用：,此即改进的欧拉法,故有公式：,3、使用泰勒
5、公式,以此方法为基础，有龙格-库塔法、线性多步法等方法。,4、数值公式的精度,当一个数值公式的截断误差可表示为O（hk+1）时（k为正整数，h为步长），称它是一个k阶公式。,k越大，则数值公式的精度越高。,欧拉法是一阶公式，改进的欧拉法是二阶公式。龙格-库塔法有二阶公式和四阶公式。线性多步法有四阶阿达姆斯外插公式和内插公式。,ODE 指令列表,MATLAB 用於求解常微分方程式的指令:,适用于 Nonstiff 系統 速率（即微分值）差异相常大 使用一般的ode45、ode23或ode113來求解,可能会使得积分的步长(Step Sizes)变得很小,以便降低积分误差至可容忍范围以內,会导致计
6、算时间过长 专门对付Stiff系统的指令，例如ode15s、ode23s、ode23t及ode23tb,指令项目繁多，最主要可分两大类：,适用于 Nonstiff 系統 一般的常微分方程式都是 Nonstiff 系統 直接采用 ode45、ode23 或 ode113 來求解,提示,使用 Simulink 來求解常微分方程式Simulink是和MATLAB共同使用的一套软件 可使用拖拉的方式來建立动力系统 可直接产生C语言代码或进行动画演示 功能非常强大,ODE 指令基本用法,使用ODE指令时，必须先将要求解的ODE表示 成一个函数 输入为t(时间)及y(状态函数:State Variable
7、s）输出则为dy（状态变量的微分值）,ODE 指令基本用法,若ODE函数的文件为odeFile.m,则调用ODE指令 的格式如下：t,y=solver(odeFile,t0,t1,y0)t0,t1 是积分的时间区间 y0代表起始条件（Initial Conditions）solver是前表所列的各种ODE指令 t是输出的时间向量 y是对应的状态变量向量,t,y=solver（f,ts,y0,options）,1、在解n个未知函数的方程组时，y0和y均为n维向量，m-文件中的待解方程组应以y的分量形式写成.,2、使用Matlab软件求数值解时，高阶微分方程必须等价地变换成一阶微分方程组.,注意:
8、,Van der Pol微分方程,1928年荷兰的范德波耳(Van der Pol)为描述LC回路的电子管振荡器建立了著名的vanderPol方程.它在自激振荡理论中有着重要的意义,一直作为数学物理方程中的一个基本方程.这是一个具有可变非线性阻尼的微分方程,代表了一类极为典型的非线性问题.和其他非线性微分方程在数学上无法精确求解一样,人们一直在努力寻找求解这类方程近似解析解的方法,并乐于用Van der Pol方程来检验求解方法的有效性.,Van der Pol微分方程其方程形式为：令,将Van der Pol微分方程化成标准形式,写成向量的形式：,为一个向量，代表状态变量,设=1，ODE文件
9、(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指令来求解,在=1时，van der Pol方程并非 Stiff系统，所以使用ode45来显示结果 odeBasic01.m,ode45(vdp1,0 25,3 3);,0 25代表积分的时间区间，3 3 则代表起始条件因为没有输出变量，所以上述程序执行结束后，MATLAB 只画出状态变量对时间的图形,odeGetData01.m,t,y=ode45(vdp1,0 25,3 3);plot(t
10、,y(:,1),t,y(:,2),:);xlabel(Time t);ylabel(Solution y(t)and y(t);legend(y(t),y(t);,也可以画出 及 在相位平面（Phase Plane）的运动情况 odePhastPlot01.m,t,y=ode45(vdp1,0 25,3 3);plot(y(:,1),y(:,2),-o);xlabel(y(t);ylabel(y(t);,当 值越来越大时，Van der Pol方程就渐变为一个Stiff系统，此时若要求解该类问题，就必须用对应于Stiff系统的指令,将 值改成 1000，ODE 文件改写成（vdp2.m）：ty
11、pe vdp2.m function dy=vdp2(t,y)mu=1000;dy=y(2);mu*(1-y(1)2)*y(2)-y(1);,选用专门用来解决 Stiff 系统问题的 ODE 指令，例如 ode15s，來求解此系统并作图显示 ode15s01.m,t,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);,由上图可知，的变化相当剧烈（超过），这
12、也就是 Stiff 系统的特点,绘制二维平面相位图：ode15s02.m若要产生在某些特定时间点的状态变量值，则调用 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),ODE 指令的选项,ODE 指令可以接受第四个输入变量，代表积分过程用到的各种选项（Options），此种ODE指令的格式为：t,y=solver(odeFile,t
13、0,tn,y0,options)其中 options 是由 odeset 指令來控制的结构变量 结构变量即包含了积分过程用到的各种选项 odeset的一般格式如下：options=odeset(name1,value1,name2,value2,)name1的值为value1，name2的值为value2,也可以只改变一个现有的options中某栏位的值，格式如下：newOptions=odeset(options,name,value);若要获取某栏位的值，可用 odeget，格式如下：value=odeget(otpions,name);当使用 odeset 指令時，若无输入任何变量名，则
14、 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 Refine:positive integer Stats:on|off InitialStep:positive scalar MaxStep:positive scalar BD
15、F: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 InitialSlope:vector Events:function_handle,由 odeset 产生的 ODE 选项,常用到的档位说明,f在积分误差容忍度方面，每一次
16、积分所产生的局部误差e(i)，必须满足下列方程式：max(RelTol*,AbsTol(i)其中i代表第i个状态 降低RelTol及AbsTol来求得更精确的积分结果,odeRelTol01.m,subplot(2,1,1);ode45(vdp1,0 25,3 3);title(RelTol=0.01);options=odeset(RelTol,0.00001);subplot(2,1,2);ode45(vdp1,0 25,3 3,options);title(RelTol=0.0001);,第一个图所使用的相对误差值是0.01（预设值），第二个图所使用的相对误差值是0.00001，因此我们得到较细密的点，但所花的计算时间也会比较长,积分输出说明,积分输出的相关处理 选用一个 OutputFcn OutputFcn 的预设值是“odeplot”，其功能为画出所有的状态变量的图像 其它可用的函式 odephas2：画出 2-D 的相位平面（Phase Plane）odephas3：画出 3-D 的相位平面odeprint：随时在指令视窗显示计算结果,options=odeset(OutputFcn,odephas2);%使用 odephas2 进行绘图ode45(vdp1,0 25,3 3,options);,
