2022年常微分方程MATLAB解法 .pdf

《2022年常微分方程MATLAB解法 .pdf》由会员分享，可在线阅读，更多相关《2022年常微分方程MATLAB解法 .pdf（9页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、精品资料欢迎下载实验四求微分方程的解一、问题背景与实验目的实际应用问题通过数学建模所归纳而得到的方程，绝大多数都是微分方程，真正能得到代数方程的机会很少另一方面，能够求解的微分方程也是十分有限的，特别是高阶方程和偏微分方程（组） 这就要求我们必须研究微分方程（组）的解法，既要研究微分方程 （组）的解析解法（精确解），更要研究微分方程 （组）的数值解法（近似解） 对微分方程（组）的解析解法(精确解 )，Matlab 有专门的函数可以用，本实验将作一定的介绍本实验将主要研究微分方程(组)的数值解法（近似解），重点介绍Euler 折线法二、相关函数（命令）及简介1dsolve(equ1,equ2,

2、)：Matlab 求微分方程的解析解 equ1、equ2、为方程（或条件）写方程（或条件）时用Dy 表示 y 关于自变量的一阶导数，用用 D2y 表示 y 关于自变量的二阶导数，依此类推2simplify(s) ：对表达式s 使用 maple 的化简规则进行化简例如：syms x simplify(sin(x)2 + cos(x)2) ans=1 3r,how=simple(s)：由于 Matlab 提供了多种化简规则， simple 命令就是对表达式s 用各种规则进行化简， 然后用r 返回最简形式， how 返回形成这种形式所用的规则例如：syms x r,how=simple(cos(x)

3、2-sin(x)2) r = cos(2*x) how = combine 4T,Y = solver( odefun,tspan,y 0) 求微分方程的数值解说明：(1) 其中的 solver 为命令 ode45、ode23、ode113 、ode15s 、ode23s 、ode23t、ode23tb 之一(2) odefun 是显式常微分方程：00)(),(ytyytfdtdy(3) 在积分区间tspan=,0ftt上，从0t到ft，用初始条件0y求解精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 1 页，共 9 页精品资料欢迎下载(4) 要

4、获 得 问 题 在 其 他 指 定 时 间 点,210,ttt上 的 解 ， 则 令tspan= ,210ftttt（要求是单调的）(5) 因为没有一种算法可以有效地解决所有的ODE 问题，为此， Matlab 提供了多种求解器Solver，对于不同的 ODE 问题，采用不同的Solver求解器Solver ODE 类型特点说明ode45 非刚性单步算法； 4、5 阶 Runge-Kutta方程；累计截断误差达3)( x大部分场合的首选算法ode23 非刚性单步算法； 2、3 阶 Runge-Kutta方程；累计截断误差达3)( x使用于精度较低的情形ode113 非刚性多步法； Adams

5、算法；高低精度均可到631010计算时间比ode45 短ode23t 适度刚性采用梯形算法适度刚性情形ode15s 刚性多步法； Gears反向数值微分；精度中等若 ode45 失效时，可尝试使用ode23s 刚性单步法； 2 阶 Rosebrock 算法；低精度当精度较低时，计算时间比 ode15s 短ode23tb 刚性梯形算法；低精度当精度较低时，计算时间比 ode15s 短(6) 要特别的是： ode23、ode45 是极其常用的用来求解非刚性的标准形式的一阶常微分方程 (组)的初值问题的解的Matlab 的常用程序，其中：ode23 采用龙格 -库塔 2 阶算法，用 3 阶公式作误差

6、估计来调节步长，具有低等的精度ode45 则采用龙格 -库塔 4 阶算法，用 5 阶公式作误差估计来调节步长，具有中等的精度5ezplot(x,y,tmin,tmax) ：符号函数的作图命令x,y 为关于参数t 的符号函数， tmin,tmax 为 t 的取值范围6inline()：建立一个内联函数格式：inline(expr, var1, var2,) ，注意括号里的表达式要加引号例：Q = dblquad(inline(y*sin(x), pi, 2*pi, 0, pi) 三、实验内容精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 2 页，共

7、9 页精品资料欢迎下载1. 几个可以直接用Matlab 求微分方程精确解的例子：例 1：求解微分方程22xxexydxdy，并加以验证求解本问题的 Matlab 程序为：syms x y %line1 y=dsolve(Dy+2*x*y=x*exp(-x2),x) %line2 diff(y,x)+2*x*y-x*exp(-x2) %line3 simplify(diff(y,x)+2*x*y-x*exp(-x2) %line4 说明：(1) 行 line1 是用命令定义 x,y 为符号变量这里可以不写， 但为确保正确性，建议写上；(2) 行 line2 是用命令求出的微分方程的解：1/2*e

8、xp(-x2)*x2+exp(-x2)*C1 (3) 行 line3 使用所求得的解这里是将解代入原微分方程，结果应该为0，但这里给出：-x3*exp(-x2)-2*x*exp(-x2)*C1+2*x*(1/2*exp(-x2)*x2+exp(-x2)*C1) (4) 行 line4 用 simplify() 函数对上式进行化简， 结果为 0， 表明)(xyy的确是微分方程的解例 2：求微分方程0 xeyxy在初始条件ey2) 1(下的特解，并画出解函数的图形求解本问题的Matlab 程序为：syms x y y=dsolve(x*Dy+y-exp(x)=0,y(1)=2*exp(1),x)

9、ezplot(y) 微分方程的特解为： y=1/x*exp(x)+1/x* exp (1) (Matlab 格式)，即xeeyx，解函数的图形如图1：-6-4-20246-30-20-1001020304050 x1/x exp(x)+1/x exp(1)图 1 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 3 页，共 9 页精品资料欢迎下载例 3：求微分方程组035yxdtdyeyxdtdxt在初始条件0|, 1|00ttyx下的特解，并画出解函数的图形求解本问题的Matlab 程序为：syms x y t x,y=dsolve(Dx+5*x+

10、y=exp(t),Dy-x-3*y=0,x(0)=1,y(0)=0,t) simple(x); simple(y); ezplot(x,y,0,1.3);axis auto 微分方程的特解 (式子特别长 )以及解函数的图形均略2. 用 ode23、 ode45等求解非刚性的标准形式的一阶常微分方程(组)的初值问题的数值解 (近似解 )例 4：求解微分方程初值问题1)0(2222yxxydxdy的数值解，求解范围为区间0, 0.5fun=inline(-2*y+2*x2+2*x,x,y); x,y=ode23(fun,0,0.5,1); x; y; plot(x,y,o-) x ans = 0.

11、0000 0.0400 0.0900 0.1400 0.1900 0.2400 0.2900 0.3400 0.3900 0.4400 0.4900 0.5000 y ans = 1.0000 0.9247 0.8434 0.7754 0.7199 0.6764 0.6440 0.6222 0.6105 0.6084 0.6154 0.6179 图形结果为图2精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 4 页，共 9 页精品资料欢迎下载00.050.10.150.20.250.30.350.40.450.50.60.650.70.750.80.

12、850.90.951图 2 例 5：求解描述振荡器的经典的Ver der Pol 微分方程.7,0)0( ,1)0(,0)1 (222yyydtdyydtyd分析：令,121dtdxxyx则.)1(,1221221xxxdtdxxdtdx先编写函数文件 verderpol.m：function xprime = verderpol(t,x) global mu; xprime = x(2);mu*(1-x(1)2)*x(2)-x(1); 再编写命令文件 vdp1.m：global mu; mu = 7; y0=1;0 t,x = ode45(verderpol,0,40,y0); x1=x(:

13、,1);x2=x(:,2); plot(t,x1) 图形结果为图 3精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 5 页，共 9 页精品资料欢迎下载0510152025303540-2.5-2-1.5-1-0.500.511.522.5图 3 3. 用 Euler 折线法求解前面讲到过， 能够求解的微分方程也是十分有限的下面介绍用Euler 折线法求微分方程的数值解 (近似解 )的方法Euler 折线法求解的基本思想是将微分方程初值问题00)(),(yxyyxfdxdy化成一个代数方程， 即差分方程，主要步骤是用差商hxyhxy)()(替代微商d

14、xdy,于是：)(),(,()()(00 xyyxyxfhxyhxykkkk记)(,1kkkkxyyhxx，从而)(1hxyykk，则有1,2, 1 ,0).,(,),(1100nkyxhfyyhxxxyykkkkkk例 6：用 Euler 折线法求解微分方程初值问题1)0(,22yyxydxdy的数值解 (步长 h 取 0.4)，求解范围为区间 0,2解：本问题的差分方程为精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 6 页，共 9 页精品资料欢迎下载1, 2, 1 ,0).2),(),(, 4.0, 1,021100nkyxyyxfyxhfy

15、yhxxhyxkkkkkk（其中：相应的 Matlab 程序见附录1数据结果为：0 1.0000 0.4000 1.4000 0.8000 2.1233 1.2000 3.1145 1.6000 4.4593 2.0000 6.3074 图形结果见图 4：00.20.40.60.811.21.41.61.821234567图 4 特别说明： 本问题可进一步利用四阶Runge-Kutta 法求解， 读者可将两个结果在一个图中显示， 并和精确值比较， 看看哪个更“精确”？（相应的Matlab 程序参见附录2） 四、自己动手1. 求微分方程0sin2)1(2xxyyx的通解2. 求微分方程xeyyy

16、xsin52 的通解3. 求微分方程组00yxdtdyyxdtdx精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 7 页，共 9 页精品资料欢迎下载在初始条件0|, 1|00ttyx下的特解，并画出解函数( )yf x的图形4. 分别用 ode23、ode45 求上述第3 题中的微分方程初值问题的数值解(近似解)，求解区间为0,2t利用画图来比较两种求解器之间的差异5. 用 Euler 折线法求解微分方程初值问题1)0(,1232yyxyy的数值解 (步长 h 取 0.1)，求解范围为区间 0,26. 用四阶 Runge-Kutta 法求解微分方程

17、初值问题1)0(,cosyxeyyx的数值解 (步长 h 取 0.1)，求解范围为区间 0,3四阶Runge-Kutta 法的迭代公式为(Euler 折线法实为一阶Runge-Kutta 法)：1, 2, 1 , 0),()2,2()2,2(),()22(6,),(342312143211100nkhLyhxfLLhyhxfLLhyhxfLyxfLLLLLhyyhxxxyykkkkkkkkkkkk相应的 Matlab 程序参见附录2试用该方法求解第5 题中的初值问题7. 用 ode45 方法求上述第6 题的常微分方程初值问题的数值解(近似解 )，从而利用画图来比较两者间的差异五、附录附录 1：

18、(fulu1.m) clear f=sym(y+2*x/y2); a=0; b=2; h=0.4; n=(b-a)/h+1; x=0; y=1; 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 8 页，共 9 页精品资料欢迎下载szj=x,y; for i=1:n-1 y=y+h*subs(f,x,y,x,y); x=x+h; szj=szj;x,y; end szj plot(szj(:,1),szj(:,2) 附录 2：(fulu2.m) clear f=sym(y-exp(x)*cos(x); a=0; b=3; h=0.1; n=(b-a)/h+1; x=0; y=1; szj=x,y; for i=1:n-1 l1=subs(f,x,y,x,y); l2=subs(f,x,y,x+h/2,y+l1*h/2); l3=subs(f,x,y,x+h/2,y+l2*h/2); l4=subs(f,x,y,x+h,y+l3*h); y=y+h*(l1+2*l2+2*l3+l4)/6; x=x+h; szj=szj;x,y; end szj plot(szj(:,1),szj(:,2) 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 9 页，共 9 页