最新微分方程(2)PPT课件.ppt

《最新微分方程(2)PPT课件.ppt》由会员分享，可在线阅读，更多相关《最新微分方程(2)PPT课件.ppt（46页珍藏版）》请在得力文库 - 分享文档赚钱的网站上搜索。

1、微分方程(2)微分方程求解方法MATLAB软件的实现微分方程的定义数学建模实例实验目的与内容MATLABMATLAB软件实现软件实现解析解dsolve(eqn1,eqn2,c1,var1,)微分方程组初值条件变量组注意：y Dy，y D2y 自变量名可以省略，默认变量名t。例输入：y=dsolve(Dy=1+y2)y1=dsolve(Dy=1+y2,y(0)=1,x)输出：y=tan(t-C1)（通解，一簇曲线）y1=tan(x+1/4*pi)（特解，一条曲线）例 常系数的二阶微分方程y=dsolve(D2y-2*Dy-3*y=0,x)y=dsolve(D2y-2*Dy-3*y=0,y(0)=

2、1,Dy(0)=0,x)输入 ：x=dsolve(D2x-(1-x2)*Dx+x=0,x(0)=3,Dx(0)=0)上述两例的计算结果怎样？由此得出什么结论？例 非常系数的二阶微分方程例无解析表达式！x=dsolve(Dx)2+x2=1,x(0)=0)例 非线性微分方程x=sin(t)-sin(t)若欲求解的某个数值解，如何求解？t=pi/2;eval(x)输入：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)例输出：（li3.m）数值解1、欧拉法2、龙格库塔法数值求解思想：(变

3、量离散化)引入自变量点列xn yn，在x0 x1x2xn上求y(xn)的近似值yn.通常取等步长 h，即xn=x0+nh，或 xn=xn-1+h，（n=1,2,)。研究常微分方程的数值解法是十分必要的研究常微分方程的数值解法是十分必要的。1)向前欧拉公式:（y=f(x,y)）y(xn+1)y(xn)+h f(xn,y(xn)(迭代式)yn+1 yn+h f(xn,yn)(近似式)特点：f(x,y)取值于区间xn,xn+1的左端点.1、欧拉方法在小区间xn,xn+1上用差商代替微商(近似),yn+1 yn+h f(xn+1,yn+1)特点：f(x,y)取值于区间xn,xn+1的右端点.非线性方程

4、,称隐式公式。yn+1=yn+h f(xn,yn)2)向后欧拉公式方法：迭代（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);(向前）end例 1 观察向前欧拉、向后欧拉算法计算情况。与精确解进行比较。误差有多大？解：1）解析解：y=x+e-x y=dsolve(Dy=-y+x+1,y(0)=1,x)2）向前欧拉法：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

5、+1+h)/(1+h)y=f(x,y)=-y+x+1;x1(1)=0;y1(1)=1;y2(1)=1;h=0.1;(died.m)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向前欧拉解，向前欧拉解，y2向后欧拉解）向后欧拉解）x=0:0.1:1;y=x+exp(-x)（解析解）（解析解）plot(x,y,x1,y1,k:,x1,y2,r-)x精确解向前欧拉向后欧拉01110.11.004811.00910.21.01871.01

6、1.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.28741.324111.36791.34871.3855（1）步长）步长h=0.1的数值解比较表的数值解比较表结果（2）步长）步长h=0.01的数值解比较表的数值解比较表x精确解向前欧拉向后欧拉01110.11.00481.00441.00530.21.01871.01791.01950.31.04081.039

7、71.04190.41.07031.06901.07170.51.10651.10501.10800.61.14881.14721.15040.71.19661.19481.19830.81.24931.24751.25110.91.30661.30471.308411.36791.36601.3697显然迭代步长显然迭代步长h 的选取对精度有影响。的选取对精度有影响。图形显示 有什么方法可以使精度提高？向后欧拉法向后欧拉法 对方程y=f(x,y),两边由xi到xi+1积分，并利用梯形公式，有：使用数值积分使用数值积分即梯形法：梯形公式梯形公式 改进的欧拉公式改进的欧拉公式yn+hf(xn,y

8、n)以例1为例，用改进欧拉公式编程计算，再与精确解的比较。yn+1=yn+(h/2)*(-yn+xn+1)+(-yn+1+xn+1+1)=yn+(h/2)*(-yn+xn+1)-(yn+h*(-yn+xn+1)+xn+1+1 =yn+(h/2)*(1-h)*xn+xn+1+2-h+(h-2)*yn died1.mx精确解向前欧拉向后欧拉改进欧拉011110.11.004811.00911.0050.21.01871.011.02641.0190.31.04081.0291.05131.04120.41.07031.05611.08301.07080.51.10651.09051.12091.1

9、0710.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的数值解比较表的数值解比较表结果使用泰勒公式使用泰勒公式 以此方法为基础，有龙格龙格-库塔法库塔法、线性多步法线性多步法等方法。数值公式的精度数值公式的精度 当一个数值公式的截断误差可表示为O（hk+1）时（k为正整数，h为步长），称它是一个k阶公式阶公式。k越大，则数值公式的精度越高。欧拉法是一阶公式，改进

10、的欧拉法是二阶公式。龙格-库塔法有二阶公式和四阶公式。线性多步法有四阶阿达姆斯外插公式和内插公式。t，x=solver（f,ts,x0,options）ode23 ode45 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：分别为设定的相对误差和绝对误差.Matlab软件

11、计算数值解1）首先建立M-文件 （weif.m）function f=weif(x,y)f=-y+x+1;2）求解：x，y=ode23(weif,0,1,1)3)作图形:plot(x,y,r);4)与精确解进行比较 hold on ezplot(x+exp(-x),0,1)例1 y=-y+x+1，y(0)=1标准形式：y=f(x,y)1、在解n个未知函数的方程组时，x0和x均为n维向量，m-函数文件中的待解方程组应以x的分量形式写成.2、使用Matlab软件求数值解时，高阶微分方程必须等价地变换成一阶微分方程组.注意注意:注意1：1、建立、建立M文件函数文件函数 function xdot=f

12、un(t,x,y)xdot=f1(t,x(t),y(t)；f2(t,x(t),y(t)；2、数值计算、数值计算（执行以下命令）（执行以下命令）t,x,y=ode23(t,x,y=ode23(funfun,t,t0 0,t,tf f,x,x0 0,y,y0 0)注意：执行命令不能写在注意：执行命令不能写在MM函数文件中。函数文件中。xd(1)=f1(t,x(t),y(t);xd(2)=f2(t,x(t),y(t);xdot=xd;%列向量列向量例如：例如：令令注意2：function xdot=fun1(t,x,y)（fun1.m)xdot=f(t,x(t),y(t)；x(t)；t,x,y=od

13、e23(t,x,y=ode23(fun1fun1,t,t0 0,t,tf f,x,x0 0,y,y0 0)M-M-文件函数如文件函数如何写呢？何写呢？注意：注意：y(t)y(t)是是原方程的解。原方程的解。x(t)x(t)只是中间只是中间变量。变量。如果方程形式是：如果方程形式是：z=f(t,z,z)?例例2 Van der pol 方程方程：令令 y1=x(t),y2=x(t);该方程是否有解析解？（1）编写M文件(文件名为 vdpol.m):function yp=vdpol(t,y);yp=y(2);(1-y(1)2)*y(2)-y(1);（2）编写程序如下:（vdj.m）t,y=ode

14、23(vdpol,0,20,3,0);y1=y(:,1);%原方程的解 y2=y(:,2);plot(t,y1,t,y2,-)%y1(t),y2(t)曲线图 pause,plot(y1,y2),grid,%相轨迹图，即y2(y1)曲线 蓝色曲线 y（1）；（原方程解）红色曲线 y（2）；计算结果例3 考虑Lorenz模型：其中参数=8/3，=10，=28解：1）编写M函数文件(lorenz.m);2)数值求解并画三维空间的相平面轨线；(ltest.m)1、lorenz.mfunction xdot=lorenz(t,x)xdot=-8/3,0,x(2);0,-10,10;-x(2),28,-1

15、*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计算结果如下图图中，x1的图形为实线(蓝），x2的图形为“*”线绿），x3的图形为“+”线（红）.取t0，tf=0，10。若自变量区间取0，20、0，40，计算结果如下：曲线呈震荡发散状三维图形的混沌状ltest.m在(x1,x2,x3)空间中,观察到如下现象：1、该曲线包含两个“圆盘”，每一个都是由螺线形轨道构成。某些轨道几乎是垂直地离开圆盘中

16、一个而进入另一个。2、随着t的增加，x(t)先绕一个圆盘几圈，然后“跳”到另一个圆盘中。绕第二个圆盘几圈，又跳回原来的圆盘。并以这样的方式继续下去，在每个圆盘上绕的圈数是随机的。思考：该空间曲线与初始点x0的选择有关吗？请观察以下实验：1）x0=0 0.1 0.1;t0,tf=0,30;解向量y2）x00=0.01 0.11 0.11;t0,tf=0,30;解向量x y x=(y1-x1,y2-x2,y3-x3)1)用MATLAB软件掌握微分方程解的三种表示方法;2)掌握求微分方程数值解的欧拉方法,了解龙格-库塔方法的思想;3)通过实例学习怎样建立微分方程模型；实验目的与内容实验目的实验内容1、求微分方程的解析解,并画出它们的图形，y=y+2x,y(0)=1,0 x1;y+ycos(x)=0,y(0)=1,y(0)=0;2、用向前欧拉公式和改进的欧拉公式求方程 y=y-2x/y,y(0)=1的数值解(0 x1,h=0.1)要求编写程序。3、Rossler微分方程组：当固定参数b=2,c=4时，试讨论随参数a由小到大变化（如a(0,0.65)而方程解的变化情况，并且画出空间曲线图形，观察空间曲线是否形成混沌状？结束语结束语谢谢大家聆听！谢谢大家聆听！46