MATLAB教程【6】微分方程课件

1、1.4.8 解常微分方程 (ordinary differential equation, ODE) 一、引言微分方程求解的数值算法有多种，常用的有Euler（欧拉法）、Runge Kutta(龙格-库塔法）。Euler法称一步法，用于一阶微分方程。0 x0 x1 x2 xn xn+1当给定步长时：所以 y1 = y0 + hf (x0,y0) y2 = y1 + hf (x1,y1) yi+1 = yi+ hf (xi,yi)第1页，共16页。Runge Kutta法 龙格-库塔法：实际上取两点斜率的平均斜率来计算的，其精度高于欧拉算法 。 下面是一个经典 的四阶龙格库塔公式：第2页，共16

2、页。二、解题方法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，YT是输出的时间列矢量，矩阵Y每个列矢量是解的一个 分量第3页，共16页。(1)建立ode函数文件% m-function, g1.

3、m function dy=g1(x,y) dy=3*x.2;例1：解一阶微分方程在区间2,4的数值解 dy/dx=3x2 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=2xcos2y y(0)=pi/4x 的积分区间y的初始条件第4页，共16页。例3 ：解二阶微分方程解：1.先将方程写成一阶微分方程 令y(1)=x, y(

4、2)=dx/dt,2.建立函数文件yjs.m并存盘 function ydot = yjs( t,y ) ydot=y(2); 4 ;3.调用函数文件解方程 T,Y=ode23(yjs,0:0.1:10,2,1);时间t 的积分区间初始条件第5页，共16页。为方便令x1=x，x2=x分别对x1,x2求一阶导数，整理后写成一阶微分方程组形式 x1=x2 x2=x2(1-x12)-x1例：x+(x2-1)x+x=0 (t=0,20; x0=0,x0=0.25)1.建立m文件wf.mfunction xdot=wf(t,x)xdot=x(2); x(2)*(1-x(1)2)-x(1);建立m文件解微

5、分方程2.给定区间、初始值;求解微分方程t,x=ode23(wf, 0,20,0,0.25)plot(t,x), figure(2),plot(x(:,1),x(:,2)第6页，共16页。例：研究有空气阻力时抛体运动的特征。比较下面三种情况下的抛体的轨道：没有空气阻力；空气阻力与速度一次方成正比；以及空气阻力与速度二次方成正比。质点运动微分方程为空气阻力的三种情况分别对应方程中参数值为b=0,0.1,0.1,p=0,0,1令y(1)=x, y(2)=dx/dt, y(3)=y , y(4)=dy/dt,将方程写成一阶微分方程组第7页，共16页。%program ddqxn.mm=1; b=0,

6、0.1,0.1; p=0,0,1; %设置参数figureaxis(0 9 0 4); %设置坐标轴的范围hold onfor 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).(

7、p/2);第8页，共16页。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的指定值第9页，共16页。 odeset Ab

8、sTol: positive scalar or vector 1e-6 绝对误差。默认1e6 BDF: on | off 在使用ODE15s时是否使用反向差分公式 RelTol: positive scalar 1e-3 相对误差。默认1e3 Events: function 设置事件 InitialStep: positive scalar 解方程时使用的初始步长 Jacobian: matrix | function 设定是否计算雅各比矩阵 JPattern: sparse matrix 若返回稀疏雅各比矩阵，为on Mass: matrix | function 返回质量矩阵 Mass

9、Singular: yes | no | maybe 质量矩阵是否为奇异矩阵 MaxOrder: 1 | 2 | 3 | 4 | 5 ODE15s解刚性方程的最高阶数 MaxStep: positive scalar 步长上界 第10页，共16页。NormControl: on | off 解向量范数的误差控制OutputSel: vector of integers 选择输出解向量哪个分量 Refine: positive integer 细化输出因子 Stats: on | off 显示计算成本统计结果 Vectorized: on | off 设定向量形式的ODE函数文件 OutputF

10、cn: function 输出方式，其选项有：odeplot 按时间顺序画出全部变量的解；odephas2 二维相空间中前两个变量的图形；odephas3 三维相空间中三个变量的图形；Odeprint 打印输出第11页，共16页。 气象学家Lorenz提出一篇论文，名叫一只蝴蝶拍一下翅膀会不会在Taxas州引起龙卷风？论述某系统如果初期条件差一点点，结果会很不稳定，他把这种现象戏称做蝴蝶效应。Lorenz为何要写这篇论文呢？ 这故事发生在1961年的某个冬天，他如往常一般在办公室操作气象电脑。平时，他只需要将温度、湿度、压力等气象数据输入，电脑就会依据三个内建的微分方程式，计算出下一刻可能的气

11、象数据，因此模拟出气象变化图。 这一天，Lorenz想更进一步了解某段纪录的後续变化，他把某时刻的气象数据重新输入电脑，让电脑计算出更多的後续结果。当时，电脑处理数据资料的数度不快，在结果出来之前，足够他喝杯咖啡并和友人闲聊一阵。在一小时後，结果出来了，不过令他目瞪口呆。结果和原资讯两相比较，初期数据还差不多，越到後期，数据差异就越大了，就像是不同的两笔资讯。而问题并不出在电脑，问题是他输入的数据差了0.000127，而这些微的差异却造成天壤之别。所以长期的准确预测天气是不可能的。例：求解lorlenz方程第12页，共16页。例：求解lorlenz方程设y(1)=x, y(2)=y, y(3)=z第13页，共16页。%program lor.maxis(10 50 -50 50 -50 50); %设置坐标轴范围view(3) %设置观察三维图形的视角hold ontitle(Lonrenz Attractor) %图的标题options =odeset(OutputFcn,odephas3) ;%设置输出方式为三维空间三变量图形t,y=ode23(lorfun,0,20,0,0,eps,options);%odefile lorfun.mfunction ydot=l