《MATLAB教程【6】微分方程ppt课件》由会员分享,可在线阅读,更多相关《MATLAB教程【6】微分方程ppt课件(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。,再 见,