微分方程的matlab求解数学建模学习教案.pptx

《微分方程的matlab求解数学建模学习教案.pptx》由会员分享，可在线阅读，更多相关《微分方程的matlab求解数学建模学习教案.pptx（44页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、会计学1微分方程微分方程(wi fn fn chn)的的matlab求解数学建模求解数学建模第一页，共44页。特殊(tsh)情形：2、一阶微分方程(wi fn fn chn)组的一般形式：初始条件：y(x0)=y0微分方程(wi fn fn chn)模型第1页/共44页第二页，共44页。图形(txng)解 tyo简单的微分方程。复杂、大型的微分方程。返 回解析(ji x)解 y=f(t)数值(shz)解(ti,yi)欧拉方法改进欧拉方法 梯形法龙格-库塔法微分方程求解方法简介第2页/共44页第三页，共44页。微分方程(wi fn fn chn)数值解1、欧拉法2、龙格库塔法数值求解思想：(变量

2、(binling)离散化)引入自变量(binling)点列xn yn，在x0 x1x2xn上求y(xn)的近似值yn.通常取等步长 h，即xn=x0+nh，或 xn=xn-1+h，（n=1,2,)。第3页/共44页第四页，共44页。1)向前欧拉公式:（y=f(x,y)）y(xn+1)y(xn)+h f(xn,y(xn)(迭代式)yn+1 yn+h f(xn,yn)(近似式)特点(tdin)：f(x,y)取值于区间xn,xn+1的左端点.在小区间xn,xn+1上用差商代替微商(近似),1、欧拉方法(fngf)第4页/共44页第五页，共44页。yn+1 yn+h f(xn+1,yn+1)特点：f(

3、x,y)取值于区间xn,xn+1的右端点(dun din).非线性方程,称隐式公式。yn+1=yn+h f(xn,yn)2)向后欧拉公式(gngsh)方法(fngf)：迭代（y=f(x,y)）x=;y=;x(1)=x0;y(1)=y0;for n=1:k x(n+1)=x(n)+n*h;y(n+1)=y(n)+h*f(x(n),y(n);(向前）end1、欧拉方法第5页/共44页第六页，共44页。例 1 观察(gunch)向前欧拉、向后欧拉算法计算情况。与精确解进行比较。误差有多大？解：1）解析(ji x)解：y=x+e-x1、欧拉方法(fngf)第6页/共44页第七页，共44页。2）向前(x

4、in qin)欧拉法：yn+1=yn+h(-yn+xn+1)=(1-h)yn+h xn+h 3）向后欧拉法：yn+1=yn+h(-yn+1+xn+1+1)转化 yn+1=(yn+h xn+1+h)/(1+h)y=f(x,y)=-y+x+1;1、欧拉方法(fngf)第7页/共44页第八页，共44页。x1(1)=0;y1(1)=1;y2(1)=1;h=0.1;()for k=1:10 x1(k+1)=x1(k)+h;y1(k+1)=(1-h)*y1(k)+h*x1(k)+h;y2(k+1)=(y2(k)+h*x1(k+1)+h)/(1+h);endx1,y1,y2,（y1向前向前(xin qin)

5、欧拉解，欧拉解，y2向后欧拉解）向后欧拉解）x=0:0.1:1;y=x+exp(-x)（解析解）（解析解）plot(x,y,x1,y1,k:,x1,y2,r-)1、欧拉方法(fngf)第8页/共44页第九页，共44页。x精确解向前欧拉向后欧拉01110.11.004811.00910.21.01871.011.02640.31.04081.0291.05130.41.07031.05611.08300.51.10651.09051.12090.61.14881.13141.16450.71.19661.17831.21320.81.24931.23051.26650.91.30661.2874

6、1.324111.36791.34871.3855（1）步长的数值）步长的数值(shz)解比较表解比较表计算结果第9页/共44页第十页，共44页。（2）步长的数值）步长的数值(shz)解比较表解比较表x精确解向前欧拉向后欧拉01110.11.00481.00441.00530.21.01871.01791.01950.31.04081.03971.04190.41.07031.06901.07170.51.10651.10501.10800.61.14881.14721.15040.71.19661.19481.19830.81.24931.24751.25110.91.30661.30471

7、.308411.36791.36601.3697结论结论(jiln)：显然迭代步长：显然迭代步长h 的选取对精度有影响。的选取对精度有影响。第10页/共44页第十一页，共44页。图形(txng)显示 有什么方法可以使精度提高？向后欧拉法向后欧拉法返 回第11页/共44页第十二页，共44页。梯形(txng)公式 改进(gijn)欧拉公式yn+hf(xn,yn)返 回第12页/共44页第十三页，共44页。x精确解向前欧拉向后欧拉改进欧拉011110.11.004811.00911.0050.21.01871.011.02641.0190.31.04081.0291.05131.04120.41.0

8、7031.05611.08301.07080.51.10651.09051.12091.10710.61.14881.13141.16451.14940.71.19661.17831.21321.19720.81.24931.23051.26651.25000.91.30661.28741.32411.307211.36791.34871.38551.3685步长 h=0.1 的数值(shz)解比较表使用(shyng)改进欧拉公式的例第13页/共44页第十四页，共44页。2、龙格-库塔法 龙格龙格-库塔法是利用泰勒展式将库塔法是利用泰勒展式将y(x+h)在在x处展开，并取其前面若干项来近似处展

9、开，并取其前面若干项来近似(jn s)y(x+h)而得到公式而得到公式y(x+h)y(x)+h j(x,y(x),h)如果如果y(xn)yn，则，则y(xn+1)的近似的近似(jn s)值为：值为：yn+1=yn+h j(xn,yn,h),n=0,1,若若 y(x+h)-y(x)+h j(x,y(x),h)=O(h p+1)，则称以上迭代公式为则称以上迭代公式为p阶公式，阶公式，p的大小反映了截断误差的高低，高阶高精度。要得到一个的大小反映了截断误差的高低，高阶高精度。要得到一个p阶公式，关键在于如何选取阶公式，关键在于如何选取j(x,y(x),h)使之满足阶的要求。使之满足阶的要求。返 回第

10、14页/共44页第十五页，共44页。微分方程(wi fn fn chn)图解法欲将微分方程解的全局信息(xnx)形象化、直观化。对 于 一 阶 微 分 方 程(wi fn fn chn)dy/dx=f(x,y)，如果给出平面上任意一点(x,y)，就能够确定出解y=f(x)在该点(x,y)处的斜率f(x,y)。从图象上看，给出平面上的一系列点，通过每一点(x0,y0)，可以画出一条通过点(x0,y0)、斜率为f(x0,y0)的短直线。这样的短直线布满整个坐标平面，形成的图形就称为斜率场或方向场。第15页/共44页第十六页，共44页。微分方程(wi fn fn chn)图解法第16页/共44页第十

11、七页，共44页。相平面(pngmin)轨迹表示微分方程的解 微分方程(wi fn fn chn)图解法 利用微分方程的数值(shz)解法，可以得到其数值(shz)解：(x(t)，y(t)在t取离散值时的取值列向量X,Y；然后分别独立地作出函数x(t)和y(t)的曲线，如图4.2，其初值条件为(5,5)。第17页/共44页第十八页，共44页。微分方程(wi fn fn chn)图解法第18页/共44页第十九页，共44页。微分方程(wi fn fn chn)图解法 如果撇开自变量的取值T，直接利用X,Y的分量作为坐标，就可以在xoy平面上画出解的轨迹(guj)，称为相平面轨迹(guj)图。第19页

12、/共44页第二十页，共44页。微分方程(wi fn fn chn)图解法第20页/共44页第二十一页，共44页。MATLAB软件(run jin)求解解析(ji x)解y=dsolve(eqn1,eqn2,c1,x)微分方程组初值条件指明变量注意：y Dy，y D2y 自变量名可以省略，默认变量名t。第21页/共44页第二十二页，共44页。例输入(shr)：y=dsolve(Dy=1+y2)y1=dsolve(Dy=1+y2,y(0)=1,x)输出(shch)：y=tan(t-C1)（通解）y1=tan(x+1/4*pi)（特解）MATLAB软件(run jin)求解第22页/共44页第二十三

13、页，共44页。例 常系数(xsh)的二阶微分方程y=dsolve(D2y-2*Dy-3*y=0,x)y=dsolve(D2y-2*Dy-3*y=0,y(0)=1,Dy(0)=0,x)输入(shr)：y=C1*exp(-x)+C2*exp(3*x)y=3/4*exp(-x)+1/4*exp(3*x)结果：第23页/共44页第二十四页，共44页。x=dsolve(D2x-(1-x2)*Dx+x=0,x(0)=3,Dx(0)=0)例 非常系数的二阶微分方程无解析(ji x)表达式！第24页/共44页第二十五页，共44页。x=dsolve(Dx)2+x2=1,x(0)=0)例 非线性微分方程(wi f

14、n fn chn)x=sin(t)-sin(t)若欲求解的某个(mu)数值解，如何求解？t=pi/2;eval(x)MATLAB软件(run jin)求解第25页/共44页第二十六页，共44页。输入(shr)：x,y=dsolve(Dx=3*x+4*y,Dy=-4*x+3*y)x,y=dsolve(Dx=3*x+4*y,Dy=-4*x+3*y,x(0)=0,y(0)=1)例输出(shch)：x=1/2*exp(7*t)-1/2*exp(-t)y=1/2*exp(-t)+1/2*exp(7*t)返 回MATLAB软件(run jin)求解第26页/共44页第二十七页，共44页。MATLAB软件(

15、run jin)求解数值(shz)解其中（1）Fun表示由微分方程(组)写成的m文件名；（2）y0表示为函数(hnsh)的初值；（3）options 用于设定误差限(可以默认)。程序为 options=odeset(reltol,rt,abstol,at)这里的rt和at分别为设定的相对误差和绝对误差。(rt=1e-7)t,y=ode45(Fun,t0，tf,y0,options)第27页/共44页第二十八页，共44页。1）首先建立M-文件 （weif.m）function f=weif(x,y)f=-y+x+1;2）求解(qi ji)：x，y=ode23(weif,0,1,1)3)作图形:p

16、lot(x,y,r);4)与精确解进行比较 hold on ezplot(x+exp(-x),0,1)例1 y=-y+x+1，y(0)=1标准(biozhn)形式：y=f(x,y)范例(fnl)第28页/共44页第二十九页，共44页。使用Matlab软件求数值解时，高阶微分方程必须等价(dngji)地变换成一阶微分方程组.注意注意(zh(zh y):y):选择一组状态变量第29页/共44页第三十页，共44页。注意(zh y)1、建立、建立M文件函数文件函数(hnsh)function xdot=fun(t,x,y)xdot=x2(t)；x3(t)；f(t,x1(t),x2(t),xn(t)；2

17、、数值计算（执行以下命令）、数值计算（执行以下命令）t,x1,x2,xn=ode45(fun,t0,tf,x1(0),x2(0),xn(0)第30页/共44页第三十一页，共44页。例例2 Van der pol 方程方程：令令 y1=x(t),y2=x(t);该方程(fngchng)无解析解！范例(fnl)第31页/共44页第三十二页，共44页。（1）编写(binxi)M文件(文件名为 vdpol.m):function yp=vdpol(t,y);yp=y(2);(1-y(1)2)*y(2)-y(1);（2）编写程序如下(rxi):（vdj.m）t,y=ode23(vdpol,0,20,3,

18、0);y1=y(:,1);%原方程的解 y2=y(:,2);plot(t,y1,t,y2,-)%y1(t),y2(t)曲线图 pause,plot(y1,y2),grid,%相轨迹图，即y2(y1)曲线范例(fnl)第32页/共44页第三十三页，共44页。蓝色曲线(qxin)y（1）；（原方程解）红色曲线(qxin)y（2）；计算结果范例(fnl)第33页/共44页第三十四页，共44页。范例(fnl)第34页/共44页第三十五页，共44页。例3 考虑(kol)Lorenz模型：其中(qzhng)参数=8/3，=10，=28解：1）编写(binxi)M函数文件(lorenz.m);2)数值求解并

19、画三维空间的相平面轨线；(ltest.m)范例第35页/共44页第三十六页，共44页。1、lorenz.mfunction xdot=lorenz(t,x)xdot=-8/3,0,x(2);0,-10,10;-x(2),28,-1*x;2、ltest.mx0=0 0 0.1;t,x=ode45(lorenz,0,10,x0);plot(t,x(:,1),-,t,x(:,2),*,t,x(:,3),+)pauseplot3(x(:,1),x(:,2),x(:,3),grid on计算结果如下计算结果如下(rxi)图图范例(fnl)第36页/共44页第三十七页，共44页。图中，x1的图形(txng

20、)为实线(蓝），x2的图形(txng)为“*”线（绿），x3的图形(txng)为“+”线（红）.取t0，tf=0，10。若自变量区间(q jin)取0，20、0，40，计算结果如下：范例(fnl)第37页/共44页第三十八页，共44页。曲线(qxin)呈震荡发散状三维图形(txng)的混沌状ltest.m第38页/共44页第三十九页，共44页。观察(gunch)结果：1、该曲线包含两个“圆盘”，每一个都是由螺线形轨道构成。某些轨道几乎(jh)是垂直地离开圆盘中一个而进入另一个。2、随着t的增加，x(t)先绕一个圆盘几圈，然后“跳”到另一个圆盘中。绕第二个圆盘几圈，又跳回原来(yunli)的圆盘

21、。并以这样的方式继续下去，在每个圆盘上绕的圈数是随机的。思考：该空间曲线与初始点x0的选择有关吗？第39页/共44页第四十页，共44页。1）x0=0 0.1 0.1;t0,tf=0,30;解向量解向量(xingling)y2）x00=0.01 0.11 0.11;t0,tf=0,30;解向量解向量(xingling)x y x=(y1-x1,y2-x2,y3-x3)返 回第40页/共44页第四十一页，共44页。1、Apollo卫星的运动(yndng)轨迹的绘制实验(shyn)内容第41页/共44页第四十二页，共44页。2、用向前欧拉公式(gngsh)和改进的欧拉公式(gngsh)求方程 y=y-2x/y,y(0)=1的数值解(0 x1,h=0.1)要求编写程序。实验(shyn)内容第42页/共44页第四十三页，共44页。3、Rossler微分方程(wi fn fn chn)组：当固定参数b=2,c=4时，试讨论随参数a由小到大变化（如a(0,0.65)而方程解的变化情况，并且画出空间曲线(qxin)图形，观察空间曲线(qxin)是否形成混沌状？4、操练一。返 回实验(shyn)内容第43页/共44页第四十四页，共44页。