《MATLAB教程【】微分方程》由会员分享，可在线阅读，更多相关《MATLAB教程【】微分方程（16页珍藏版）》请在金锄头文库上搜索。
1、1.4.8 解常微分方程 (ordinary differential equation, ODE),一、引言 微分方程求解的数值算法有多种，常用的有Euler（欧拉法）、Runge Kutta(龙格-库塔法）。 Euler法称一步法，用于一阶微分方程。,当给定步长时：,所以 y1 = y0 + hf (x0,y0) y2 = y1 + hf (x1,y1) yi+1 = yi+ hf (xi,yi),Runge Kutta法,龙格-库塔法：实际上取两点斜率的平均斜率来计算的，其精度高于欧拉算法 。,二、解题方法 1.编写odefile 文件格式：function ydot=odefile(t
2、,y) ydot=待求的函数,(1)建立ode函数文件 % m-function, g1.m function dy=g1(x,y) dy=3*x.2;,(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);,x 的积分区间,y的初始条件,例3 ：解二阶微分方程,解：1.先将方程写成一阶微分方程 令y(1)=x, y(2)=dx/dt,2.建立函数文件yjs.m并存盘 function ydot = y
3、js( 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.m function xdot=wf(t,x) xdot=x(2); x(2)*(1-x(1)2)-x(1);,建立m文件 解微分方程,2.给定区间、初始值;求解微分方程 t,x=ode23(wf, 0,2
4、0,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.m m=1; b=0,0.1,0.1; p=0,0,1; %设置参数 figure axis(0 9 0 4); %设
5、置坐标轴的范围 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);,3.确定求解的条件和要求 T,Y=ode23(F,tspan,y0,
6、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 绝对误差。默认1e6
7、 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
8、: 1 | 2 | 3 | 4 | 5 ODE15s解刚性方程的最高阶数 MaxStep: positive scalar 步长上界,气象学家Lorenz提出一篇论文，名叫一只蝴蝶拍一下翅膀会不会在Taxas州引起龙卷风？论述某系统如果初期条件差一点点，结果会很不稳定，他把这种现象戏称做蝴蝶效应。Lorenz为何要写这篇论文呢？ 这故事发生在1961年的某个冬天，他如往常一般在办公室操作气象电脑。平时，他只需要将温度、湿度、压力等气象数据输入，电脑就会依据三个内建的微分方程式，计算出下一刻可能的气象数据，因此模拟出气象变化图。 这一天，Lorenz想更进一步了解某段纪录的後续变化，他把某时刻的
9、气象数据重新输入电脑，让电脑计算出更多的後续结果。当时，电脑处理数据资料的数度不快，在结果出来之前，足够他喝杯咖啡并和友人闲聊一阵。在一小时後，结果出来了，不过令他目瞪口呆。结果和原资讯两相比较，初期数据还差不多，越到後期，数据差异就越大了，就像是不同的两笔资讯。而问题并不出在电脑，问题是他输入的数据差了0.000127，而这些微的差异却造成天壤之别。所以长期的准确预测天气是不可能的。,例：求解lorlenz方程,例：求解lorlenz方程,%program lor.m axis(10 50 -50 50 -50 50); %设置坐标轴范围 view(3) %设置观察三维图形的视角 hold
10、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);,作业： 1.求微分方程在1,3区间内的数值解，并将结果与解析解（y=x2+(1/x2)进行比较。,2.求微分方程在区间0,2的解。,3.上机演示例1和例2。,再 见,
