MATLAB教程【6】微分方程

1.4.8解常微分方程(ordinarydifferentialequation,ODE)一、引言微分方程求解的数值算法有多种，常用的有Euler（欧拉法）、RungeKutta(龙格-库塔法）。Euler法称一步法，用于一阶微分方程。00yxyyxfdxdy)(),,(0x0x1x2xnxn+1当给定步长时：hyyxxyydxdyf010101(x)所以y1=y0+h·f(x0,y0)y2=y1+h·f(x1,y1)…yi+1=yi+h·f(xi,yi)RungeKutta法),(),()(21121211hKyxfKyxfKkkhyynnnnnn龙格-库塔法：实际上取两点斜率的平均斜率来计算的，其精度高于欧拉算法。)(4);););()(23121432113nnnnnnnnnnhKh,yxfKK2hy2hxfKK2hy2hxfK,yxfKK2K2KK6hyy,(,(下面是一个经典的四阶龙格－库塔公式：二、解题方法1.编写odefile文件格式：functionydot=odefile(t,y)ydot=[待求的函数]2.选择和调用解常微分方程的指令(常用的有ode23、ode45)指令格式：[T,Y]=ode23(‘F’,tspan,y0,options,p1,p2,…)F求解的odefile的文件名tspan单调递增(减)的积分区间y0初始条件矢量options用odeset建立的优化选项，如用默认值则不必输入p1,p2传递给F的参数，T，YT是输出的时间列矢量，矩阵Y每个列矢量是解的一个分量(1)建立ode函数文件%m-function,g1.mfunctiondy=g1(x,y)dy=3*x.^2;例1：解一阶微分方程在区间[2,4]的数值解dy/dx=3x2y(2)=0.5(2)调用函数文件解微分方程[x,num_y]=ode23('g1',[2,4],0.5);%m-function,g3.mfunctiondy=g3(x,y)dy=2*x*cos(y)^2;[x,num_y]=ode23('g3',[0,2],pi/4);例2：解一阶微分方程在区间[0，2]的数值解dy/dx=2xcos2yy(0)=pi/4x的积分区间y的初始条件例3：解二阶微分方程12)4(0022vxaadtxd解：1.先将方程写成一阶微分方程令y(1)=x,y(2)=dx/dt,4)2()2()1(dtdyydtdy2.建立函数文件yjs.m并存盘functionydot=yjs(t,y)ydot=[y(2);4];3.调用函数文件解方程[T,Y]=ode23('yjs',[0:0.1:10],[2,1]);时间t的积分区间初始条件为方便令x1=x，x2=x’分别对x1,x2求一阶导数，整理后写成一阶微分方程组形式x1’=x2x2’=x2(1-x12)-x1例：x’’+(x2-1)x’+x=0(t=[0,20];x0=0,x0’=0.25)1.建立m文件wf.mfunctionxdot=wf(t,x)xdot=[x(2);x(2)*(1-x(1)^2)-x(1)];1.建立m文件2.解微分方程2.给定区间、初始值;求解微分方程[t,x]=ode23('wf',[0,20],[0,0.25])plot(t,x),figure(2),plot(x(:,1),x(:,2))例：研究有空气阻力时抛体运动的特征。比较下面三种情况下的抛体的轨道：没有空气阻力；空气阻力与速度一次方成正比；以及空气阻力与速度二次方成正比。·vgr2p22vbmtddm质点运动微分方程为空气阻力的三种情况分别对应方程中参数值为b=[0,0.1,0.1],p=[0,0,1]令y(1)=x,y(2)=dx/dt,y(3)=y,y(4)=dy/dt,将方程写成一阶微分方程组(4)d(3)d(2)d(1)dyyyytt222222(4))(2)((4)(4)(4))(2)((2)(2)ppyymbygdtdyyymbydtdy%%programddqxn.mm=1;b=[0,0.1,0.1];p=[0,0,1];%设置参数figureaxis([0904]);%设置坐标轴的范围holdonfori=1:3%解微分方程三次[t,y]=ode45('ddqxnfun',[0:0.001:2],[0,5,0,8],[],b(i),p(i),m);comet(y(:,1),y(:,3))%画轨道的彗星轨迹end函数文件ddqxnfun.m如下：functionydot=ddqxnfun(t,y,flag,b,p,m)ydot=[y(2)；-b/m*y(2)*(y(2).^2+y(4).^2).^(p/2);y(4);-9.8-b/m*y(4)*(y(2).^2+y(4).^2).^(p/2)];3.确定求解的条件和要求[T,Y]=ode23(‘F’,tspan,y0,options,p1,p2,…)Odeset建立和修改优化选项语句格式：options=odeset(‘name1’,value1,‘name2’,value2,…)指定各个参数的取值，对不指定取值的参数，取默认值options=odeset(odeopts,‘name1’,value1,…)使用原来的优化选项，但对其中部分参数指定新值options=odeset(oldopts,newopts)合并了两个优化选项oldopts和newopts，重复部分取newopts的指定值odesetAbsTol:[positivescalarorvector{1e-6}]绝对误差。默认1e－6BDF:[on|{off}]在使用ODE15s时是否使用反向差分公式RelTol:[positivescalar{1e-3}]相对误差。默认1e－3Events:[function]设置事件InitialStep:[positivescalar]解方程时使用的初始步长Jacobian:[matrix|function]设定是否计算雅各比矩阵JPattern:[sparsematrix]若返回稀疏雅各比矩阵，为‘on’Mass:[matrix|function]返回质量矩阵MassSingular:[yes|no|{maybe}]质量矩阵是否为奇异矩阵MaxOrder:[1|2|3|4|{5}]ODE15s解刚性方程的最高阶数MaxStep:[positivescalar]步长上界NormControl:[on|{off}]解向量范数的误差控制OutputSel:[vectorofintegers]选择输出解向量哪个分量Refine:[positiveinteger]细化输出因子Stats:[on|{off}]显示计算成本统计结果Vectorized:[on|{off}]设定向量形式的ODE函数文件OutputFcn:[function]输出方式，其选项有：odeplot按时间顺序画出全部变量的解；odephas2二维相空间中前两个变量的图形；odephas3三维相空间中三个变量的图形；Odeprint打印输出气象学家Lorenz提出一篇论文，名叫「一只蝴蝶拍一下翅膀会不会在Taxas州引起龙卷风？」论述某系统如果初期条件差一点点，结果会很不稳定，他把这种现象戏称做「蝴蝶效应」。Lorenz为何要写这篇论文呢？这故事发生在1961年的某个冬天，他如往常一般在办公室操作气象电脑。平时，他只需要将温度、湿度、压力等气象数据输入，电脑就会依据三个内建的微分方程式，计算出下一刻可能的气象数据，因此模拟出气象变化图。这一天，Lorenz想更进一步了解某段纪录的後续变化，他把某时刻的气象数据重新输入电脑，让电脑计算出更多的後续结果。当时，电脑处理数据资料的数度不快，在结果出来之前，足够他喝杯咖啡并和友人闲聊一阵。在一小时後，结果出来了，不过令他目瞪口呆。结果和原资讯两相比较，初期数据还差不多，越到後期，数据差异就越大了，就像是不同的两笔资讯。而问题并不出在电脑，问题是他输入的数据差了0.000127，而这些微的差异却造成天壤之别。所以长期的准确预测天气是不可能的。例：求解lorlenz方程例：求解lorlenz方程zyxydtdzzydtdyyzxdtdx35101038zyxydtdyzydtdyyzxdtdy35(3)1010(2)38(1)设y(1)=x,y(2)=y,y(3)=z%%programlor.maxis([1050-5050-5050]);%设置坐标轴范围view(3)%设置观察三维图形的视角holdontitle('LonrenzAttractor')%图的标题options=odeset('OutputFcn','odephas3');%设置输出方式为三维空间三变量图形[t,y]=ode23('lorfun',[0,20],[0,0,eps],options);%%odefilelorfun.mfunctionydot=lorfun(t,y)ydot=[-8/3*y(1)+y(2)*y(3);-10*y(2)+10*y(3);-y(2)*y(1)+35*y(2)-y(3)];作业：1.求微分方程在[1,3]区间内的数值解，并将结果与解析解（y=x2+(1/x2))进行比较。2(1)42yxxydxdy605005222)()('yyxdxdyxdxyd2.求微分方程在区间[0,2]的解。3.上机演示例1和例2。