1.4.8 解常微分方程解常微分方程 (ordinary differential equation, ODE) 一、引言一、引言v微分方程求解的数值算法有多种，常用的有微分方程求解的数值算法有多种，常用的有Euler（欧拉法）、（欧拉法）、Runge Kutta(龙格龙格-库塔法）。库塔法）。vEuler法称一步法，用于一阶微分方程。法称一步法，用于一阶微分方程。0 x0 x1 x2 xn xn+1v当给定步长时：当给定步长时：所以所以 y1 = y0 + hf (x0,y0) y2 = y1 + hf (x1,y1) yi+1 = yi+ hf (xi,yi)MATLAB教程6微分方程vRunge Kutta法法 龙格龙格-库塔法：实际上取库塔法：实际上取两点斜率的平均斜率来计算两点斜率的平均斜率来计算的，其精度高于欧拉算法的，其精度高于欧拉算法 。 下面是一个经典下面是一个经典 的四阶龙格库塔公式：的四阶龙格库塔公式：MATLAB教程6微分方程二、解题方法二、解题方法1.编写编写odefile文件格式：文件格式：function ydot=ode) 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每个列矢量是解的一个每个列矢量是解的一个 分量分量MATLAB教程6微分方程(1)建立建立ode函数文件函数文件% m-function, g1.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的初始条件的初始条件MATLAB教程6微分方程例例3 ：解二阶微分方程：解二阶微分方程解：解：1.先将方程写成一阶微分方程先将方程写成一阶微分方程 令令y(1)=x, y(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 的积的积分区间分区间初始条件初始条件MATLAB教程6微分方程为方便令为方便令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);1.建立建立m文件文件2.解微分方程解微分方程2.给定区间、初始值;求解微分方程t,x=ode23(wf, 0,20,0,0.25)plot(t,x), figure(2),plot(x(:,1),x(:,2)MATLAB教程6微分方程例：研究有空气阻力时抛体运动的特征。比较下面三种情况例：研究有空气阻力时抛体运动的特征。比较下面三种情况下的抛体的轨道：没有空气阻力；空气阻力与速度一次方成下的抛体的轨道：没有空气阻力；空气阻力与速度一次方成正比；以及空气阻力与速度二次方成正比。正比；以及空气阻力与速度二次方成正比。质点运动微分方程为质点运动微分方程为空气阻力的三种情况分别对应方程中参数值为空气阻力的三种情况分别对应方程中参数值为b=0,0.1,0.1,p=0,0,1令令y(1)=x, y(2)=dx/dt, y(3)=y , y(4)=dy/dt,将方程写成一阶微分将方程写成一阶微分方程组方程组MATLAB教程6微分方程%program ddqxn.mm=1; b=0,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).(p/2);MATLAB教程6微分方程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的指定值的指定值MATLAB教程6微分方程 odeset AbsTol: 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 返回质量矩阵返回质量矩阵 MassSingular: yes | no | maybe 质量矩阵是否为奇异矩阵质量矩阵是否为奇异矩阵 MaxOrder: 1 | 2 | 3 | 4 | 5 ODE15s解刚性方程的最高阶数解刚性方程的最高阶数 MaxStep: positive scalar 步长上界步长上界 MATLAB教程6微分方程NormControl: on | off 解向量范数的误差控制解向量范数的误差控制OutputSel: vector of integers 选择输出解向量哪个分量选择输出解向量哪个分量 Refine: positive integer 细化输出因子细化输出因子 Stats: on | off 显示计算成本统计结果显示计算成本统计结果 Vectorized: on | off 设定向量形式的设定向量形式的ODE函数文件函数文件 OutputFcn: function 输出方式，其选项有：输出方式，其选项有：odeplot 按时间顺序画出全部变量的解；按时间顺序画出全部变量的解；odephas2 二维相空间中前两个变量的图形；二维相空间中前两个变量的图形；odephas3 三维相空间中三个变量的图形；三维相空间中三个变量的图形；Odeprint 打印输出打印输出MATLAB教程6微分方程 气象学家气象学家Lorenz提出一篇论文，名叫一只蝴蝶拍一下翅膀会不会在提出一篇论文，名叫一只蝴蝶拍一下翅膀会不会在Taxas州引起龙卷风？论述某系统如果初期条件差一点点，结果会很不稳州引起龙卷风？论述某系统如果初期条件差一点点，结果会很不稳定，他把这种现象戏称做蝴蝶效应。定，他把这种现象戏称做蝴蝶效应。Lorenz为何要写这篇论文呢？为何要写这篇论文呢？ 这故事发生在这故事发生在1961年的某个冬天，他如往常一般在办公室操作气象电年的某个冬天，他如往常一般在办公室操作气象电脑。平时，他只需要将温度、湿度、压力等气象数据输入，电脑就会依据脑。平时，他只需要将温度、湿度、压力等气象数据输入，电脑就会依据三个内建的微分方程式，计算出下一刻可能的气象数据，因此模拟出气象三个内建的微分方程式，计算出下一刻可能的气象数据，因此模拟出气象变化图。变化图。 这一天，这一天，Lorenz想更进一步了解某段纪录的後续变化，他把某时刻的气想更进一步了解某段纪录的後续变化，他把某时刻的气象数据重新输入电脑，让电脑计算出更多的後续结果。当时，电脑处理数象数据重新输入电脑，让电脑计算出更多的後续结果。当时，电脑处理数据资料的数度不快，在结果出来之前，足够他喝杯咖啡并和友人闲聊一阵。据资料的数度不快，在结果出来之前，足够他喝杯咖啡并和友人闲聊一阵。在一小时後，结果出来了，不过令他目瞪口呆。结果和原资讯两相比较，在一小时後，结果出来了，不过令他目瞪口呆。结果和原资讯两相比较，初期数据还差不多，越到後期，数据差异就越大了，就像是不同的两笔资初期数据还差不多，越到後期，数据差异就越大了，就像是不同的两笔资讯。而问题并不出在电脑，问题是他输入的数据差了讯。而问题并不出在电脑，问题是他输入的数据差了0.000127，而这些微，而这些微的差异却造成天壤之别。所以长期的准确预测天气是不可能的。的差异却造成天壤之别。所以长期的准确预测天气是不可能的。例：求解lorlenz方程MATLAB教程6微分方程例：求解lorlenz方程设设y(1)=x, y(2)=y, y(3)=zMATLAB教程6微分方程%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);%odefunction 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);MATLAB教程6微分方程作业：1.求微分方程在求微分方程在1,3区间内的区间内的数值解，并将结果与解析解数值解，并将结果与解析解（y=x2+(1/x2)进行比较。进行比较。2.求微分方程在区间求微分方程在区间0,2的解。的解。3.上机演示例上机演示例1和例和例2。MATLAB教程6微分方程MATLAB教程6微分方程