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