常微分方程MATLAB解法_第1页
常微分方程MATLAB解法_第2页
常微分方程MATLAB解法_第3页
常微分方程MATLAB解法_第4页
常微分方程MATLAB解法_第5页
已阅读5页,还剩4页未读 继续免费阅读

下载本文档

版权说明:本文档由用户提供并上传,收益归属内容提供方,若内容存在侵权,请进行举报或认领

文档简介

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

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

3、)2-sin(x)2)r = cos(2*x)how = combine4T,Y = solver(odefun,tspan,y0) 求微分方程的数值解说明:(1) 其中的 solver为命令 ode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb 之一(2) odefun 是显式常微分方程:(3) 在积分区间 tspan=上,从到,用初始条件求解(4) 要获得问题在其他指定时间点上的解,则令 tspan= (要求是单调的)(5) 因为没有一种算法可以有效地解决所有的 ODE 问题,为此,Matlab 提供了多种求解器 Solver,对于不同的ODE 问

4、题,采用不同的Solver求解器SolverODE类型特点说明ode45非刚性 单步算法;4、5阶Runge-Kutta方程;累计截断误差达大部分场合的首选算法ode23非刚性单步算法;2、3阶Runge-Kutta方程;累计截断误差达使用于精度较低的情形ode113非刚性多步法;Adams算法;高低精度均可到计算时间比 ode45 短ode23t适度刚性采用梯形算法适度刚性情形ode15s刚性多步法;Gear's反向数值微分;精度中等若 ode45 失效时,可尝试使用ode23s刚性单步法;2阶 Rosebrock 算法;低精度当精度较低时,计算时间比 ode15s 短ode23tb

5、刚性梯形算法;低精度当精度较低时,计算时间比 ode15s 短(6) 要特别的是:ode23、ode45 是极其常用的用来求解非刚性的标准形式的一阶常微分方程(组)的初值问题的解的 Matlab 的常用程序,其中:ode23 采用龙格-库塔2 阶算法,用3 阶公式作误差估计来调节步长,具有低等的精度ode45 则采用龙格-库塔4 阶算法,用5 阶公式作误差估计来调节步长,具有中等的精度5ezplot(x,y,tmin,tmax):符号函数的作图命令x,y 为关于参数t 的符号函数,tmin,tmax 为 t 的取值范围6inline():建立一个内联函数格式:inline('expr&

6、#39;, 'var1', 'var2',) ,注意括号里的表达式要加引号例:Q = dblquad(inline('y*sin(x)'), pi, 2*pi, 0, pi)三、实验内容1. 几个可以直接用 Matlab 求微分方程精确解的例子:例1:求解微分方程,并加以验证求解本问题的Matlab 程序为:syms x y %line1y=dsolve('Dy+2*x*y=x*exp(-x2)','x') %line2diff(y,x)+2*x*y-x*exp(-x2) %line3simplify(diff(y

7、,x)+2*x*y-x*exp(-x2) %line4说明:(1) 行line1是用命令定义x,y为符号变量这里可以不写,但为确保正确性,建议写上;(2) 行line2是用命令求出的微分方程的解:1/2*exp(-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, 表明的确是微分方程的解例2:求微分方程在初始条件下的特解,

8、并画出解函数的图形求解本问题的 Matlab 程序为:syms x yy=dsolve('x*Dy+y-exp(x)=0','y(1)=2*exp(1)','x')ezplot(y)微分方程的特解为:y=1/x*exp(x)+1/x* exp (1) (Matlab格式),即,解函数的图形如图 1:图1例3:求微分方程组在初始条件下的特解,并画出解函数的图形求解本问题的 Matlab 程序为:syms x y tx,y=dsolve('Dx+5*x+y=exp(t)','Dy-x-3*y=0','x(0)=

9、1','y(0)=0','t')simple(x);simple(y);ezplot(x,y,0,1.3);axis auto微分方程的特解(式子特别长)以及解函数的图形均略2. 用ode23、ode45等求解非刚性的标准形式的一阶常微分方程(组)的初值问题的数值解(近似解)例4:求解微分方程初值问题的数值解,求解范围为区间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-

10、')>> x'ans =0.0000 0.0400 0.0900 0.1400 0.1900 0.24000.2900 0.3400 0.3900 0.4400 0.4900 0.5000>> y'ans =1.0000 0.9247 0.8434 0.7754 0.7199 0.67640.6440 0.6222 0.6105 0.6084 0.6154 0.6179图形结果为图 2图2例 5:求解描述振荡器的经典的 Ver der Pol 微分方程 分析:令则先编写函数文件verderpol.m:function xprime = verde

11、rpol(t,x)global mu;xprime = x(2);mu*(1-x(1)2)*x(2)-x(1);再编写命令文件vdp1.m:global mu;mu = 7;y0=1;0t,x = ode45('verderpol',0,40,y0);x1=x(:,1);x2=x(:,2);plot(t,x1)图形结果为图3图33. 用 Euler 折线法求解前面讲到过,能够求解的微分方程也是十分有限的下面介绍用 Euler 折线法求微分方程的数值解(近似解)的方法Euler 折线法求解的基本思想是将微分方程初值问题化成一个代数方程,即差分方程,主要步骤是用差商替代微商,于是:

12、记,从而,则有例 6:用 Euler 折线法求解微分方程初值问题的数值解(步长h取0.4),求解范围为区间0,2解:本问题的差分方程为相应的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:图4特别说明:本问题可进一步利用四阶 Runge-Kutta 法求解,读者可将两个结果在一个图中显示,并和精确值比较,看看哪个更“精确”?(相应的 Matlab 程序参见附录 2)四、自己动手1. 求微分方程的通解2. 求微分方程的通解3. 求微

13、分方程组 在初始条件下的特解,并画出解函数的图形4. 分别用 ode23、ode45 求上述第 3 题中的微分方程初值问题的数值解(近似解),求解区间为利用画图来比较两种求解器之间的差异5. 用 Euler 折线法求解微分方程初值问题的数值解(步长h取0.1),求解范围为区间0,26. 用四阶 Runge-Kutta 法求解微分方程初值问题的数值解(步长h取0.1),求解范围为区间0,3四阶 Runge-Kutta 法的迭代公式为(Euler 折线法实为一阶 Runge-Kutta 法):相应的 Matlab 程序参见附录 2试用该方法求解第5题中的初值问题7. 用 ode45 方法求上述第

14、6 题的常微分方程初值问题的数值解(近似解),从而利用画图来比较两者间的差异五、附录附录 1:(fulu1.m)clearf=sym('y+2*x/y2');a=0;b=2;h=0.4;n=(b-a)/h+1;x=0;y=1;szj=x,y;for i=1:n-1y=y+h*subs(f,'x','y',x,y);x=x+h;szj=szj;x,y;endszjplot(szj(:,1),szj(:,2)附录 2:(fulu2.m)clearf=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-1l1=subs(f,'x','y'

温馨提示

  • 1. 本站所有资源如无特殊说明,都需要本地电脑安装OFFICE2007和PDF阅读器。图纸软件为CAD,CAXA,PROE,UG,SolidWorks等.压缩文件请下载最新的WinRAR软件解压。
  • 2. 本站的文档不包含任何第三方提供的附件图纸等,如果需要附件,请联系上传者。文件的所有权益归上传用户所有。
  • 3. 本站RAR压缩包中若带图纸,网页内容里面会有图纸预览,若没有图纸预览就没有图纸。
  • 4. 未经权益所有人同意不得将文件中的内容挪作商业或盈利用途。
  • 5. 人人文库网仅提供信息存储空间,仅对用户上传内容的表现方式做保护处理,对用户上传分享的文档内容本身不做任何修改或编辑,并不能对任何下载内容负责。
  • 6. 下载文件中如有侵权或不适当内容,请与我们联系,我们立即纠正。
  • 7. 本站不保证下载资源的准确性、安全性和完整性, 同时也不承担用户因使用这些下载资源对自己和他人造成任何形式的伤害或损失。

评论

0/150

提交评论