MATLAB教程微分方程课件

单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,#,Runge Kutta,法,龙格,-,库塔法：实际上取两点斜率的平均斜率来计算的，其精度高于欧拉算法。,下面是一个经典 的四阶龙格库塔公式：,Runge Kutta法 龙格-库塔法：实际,1,二、解题方法,1.,编写,odefile,文件格式：,function ydot=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,，,Y,T,是输出的时间列矢量，矩阵,Y,每个列矢量是解的一个 分量,二、解题方法2.选择和调用解常微分方程的指令(常用的有od,2,(1),建立,ode,函数文件,%m-function,g1.m,function dy=g1(x,y),dy=3*x.2;,例,1,：解一阶微分方程在区间,2,4,的数值解,d,y,/d,x,=3,x,2,y,(2)=0.5,(2),调用函数文件解微分方程,x,num_y=ode23(g1,2,4,0.5);,%m-function,g3.m,function dy=g3(x,y),dy=2*x*cos(y)2;,x,num_y=ode23(g3,0,2,pi/4);,例,2,：解一阶微分方程在区,间,0,，,2,的数值解,dy/dx=2xcos,2,y,y(0)=pi/4,x,的积分区间,y,的初始条件,(1)建立ode函数文件例1：解一阶微分方程在区间2,4,3,例,3,：解二阶微分方程,解：,1.,先将方程写成一阶微分方程,令,y,(1)=,x,y,(2)=d,x,/d,t,2.,建立函数文件,yjs.m,并存盘,function ydot=yjs(t,y),ydot=y(2);,4;,3.,调用函数文件解方程,T,Y=ode23(yjs,0:0.1:10,2,1);,时间,t,的积,分区间,初始条件,例3：解二阶微分方程解：1.先将方程写成一阶微分方程2.建,4,为方便令,x,1,=x,，,x,2,=x,分别对,x,1,x,2,求一,阶导数，整理后写成一阶微分方程组,形式,x,1,=x,2,x,2,=x,2,(1-x,1,2,)-x,1,例：,x+(x,2,-1)x+x=0(,t=0,20;x,0,=0,x,0,=0.25),1.,建立,m,文件,wf.m,function xdot=wf(t,x),xdot=x(2);,x(2)*,(1-x(1)2)-x(1);,建立,m,文件,解微分方程,2.,给定区间、初始值,;,求解微分方程,t,x=ode23(wf,0,20,0,0.25),plot(t,x),figure(2),plot(x(:,1),x(:,2),为方便令x1=x，x2=x分别对x1,x2求一例：x+,5,例：研究有空气阻力时抛体运动的特征。比较下面三种情况下的抛体的轨道：没有空气阻力；空气阻力与速度一次方成正比；以及空气阻力与速度二次方成正比。,质点运动微分方程为,空气阻力的三种情况分别对应方程中参数值为,b=0,0.1,0.1,p=0,0,1,令,y,(1)=,x,y,(2)=d,x,/d,t,y,(3)=,y,y,(4)=d,y,/d,t,将方程写成一阶微分方程组,例：研究有空气阻力时抛体运动的特征。比较下面三种情况下的抛体,6,%program ddqxn.m,m=1;b=0,0.1,0.1;p=0,0,1;,%,设置参数,figure,axis(0 9 0 4);,%,设置坐标轴的范围,hold on,for i=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,如下：,function ydot=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);,%program ddqxn.m函数文件ddqxnfun,7,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,的指定值,3.确定求解的条件和要求,8,odeset,AbsTol:,positive scalar or vector 1e-6,绝对误差。默认,1e,6,BDF:on|off,在使用,ODE15s,时是否使用反向差分公式,RelTol,:positive scalar 1e-3,相对误差。默认,1e,3,Events:function,设置事件,InitialStep:positive scalar,解方程时使用的初始步长,Jacobian:matrix|function,设定是否计算雅各比矩阵,JPattern:sparse matrix,若返回稀疏雅各比矩阵，为,on,Mass:matrix|function,返回质量矩阵,MassSingular:yes|no|maybe,质量矩阵是否为奇异矩阵,MaxOrder:1|2|3|4|5,ODE15s,解刚性方程的最高阶数,MaxStep:positive scalar,步长上界,odeset,9,NormControl:on|off,解向量范数的误差控制,OutputSel:vector of integers,选择输出解向量哪个分量,Refine:positive integer,细化输出因子,Stats:on|off,显示计算成本统计结果,Vectorized:on|off,设定向量形式的,ODE,函数文件,OutputFcn:function,输出方式，其选项有：,odeplot,按时间顺序画出全部变量的解；,odephas2,二维相空间中前两个变量的图形；,odephas3,三维相空间中三个变量的图形；,Odeprint,打印输出,NormControl:on|off,10,气象学家,Lorenz,提出一篇论文，名叫一只蝴蝶拍一下翅膀会不会在,Taxas,州引起龙卷风？论述某系统如果初期条件差一点点，结果会很不稳定，他把这种现象戏称做蝴蝶效应。,Lorenz,为何要写这篇论文呢？,这故事发生在,1961,年的某个冬天，他如往常一般在办公室操作气象电脑。平时，他只需要将温度、湿度、压力等气象数据输入，电脑就会依据三个内建的微分方程式，计算出下一刻可能的气象数据，因此模拟出气象变化图。,这一天，,Lorenz,想更进一步了解某段纪录的後续变化，他把某时刻的气象数据重新输入电脑，让电脑计算出更多的後续结果。当时，电脑处理数据资料的数度不快，在结果出来之前，足够他喝杯咖啡并和友人闲聊一阵。在一小时後，结果出来了，不过令他目瞪口呆。结果和原资讯两相比较，初期数据还差不多，越到後期，数据差异就越大了，就像是不同的两笔资讯。而问题并不出在电脑，问题是他输入的数据差了,0.000127,，而这些微的差异却造成天壤之别。所以长期的准确预测天气是不可能的。,例：求解,lorlenz,方程,气象学家Lorenz提出一篇论文，名叫一只蝴,11,例：求解,lorlenz,方程,设,y,(1)=,x,y,(2)=,y,y,(3)=,z,例：求解lorlenz方程设y(1)=x,y(2)=y,12,%program lor.m,axis(10 50-50 50-50 50);,%,设置坐标轴范围,view(3),%,设置观察三维图形的视角,hold on,title(Lonrenz Attractor),%,图的标题,options=odeset(OutputFcn,odephas3);,%,设置输出方式为三维空间三变量图形,t,y=ode23(lorfun,0,20,0,0,eps,options);,%odefile lorfun.m,function ydot=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);,%program lor.m%odefile lorfu,13,作业：,1.,求微分方程在,1,3,区间内的数值解，并将结果与解析解（,y=x,2,+(1/x,2,),进行比较。,2.,求微分方程在区间,0,2,的解。,3.,上机演示例,1,和例,2,。,作业：2.求微分方程在区间0,2的解。3.上机演示例1和,14,再 见,再 见,15,