MATLAB教程【6】微分方程.pptVIP

1.4.8 解常微分方程 (ordinary differential equation, ODE)  一、引言 微分方程求解的数值算法有多种，常用的有Euler（欧拉法）、Runge Kutta(龙格-库塔法）。 Euler法称一步法，用于一阶微分方程。 Runge Kutta法 为方便令x1=x，x2=x’分别对x1,x2求一 阶导数，整理后写成一阶微分方程组 形式 x1’=x2 x2’=x2(1-x12)-x1 * * 0 x0 x1 x2 xn xn+1 当给定步长时： 所以 y1 = y0 + h·f (x0,y0) y2 = y1 + h·f (x1,y1) … yi+1 = yi+ h·f (xi,yi) 龙格-库塔法：实际上取两点斜率的平均斜率来计算的，其精度高于欧拉算法 。 下面是一个经典 的四阶龙格－库塔公式： 二、解题方法 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/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/4 x 的积分区间 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 的积 分区间 初始条件 例：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,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.m m=1; b=[0,0.1,0.1]; p=[0,0,1]; %设置参数 figure axis([0 9 0 4]); %设置坐标轴的范围 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)) %画轨