1、1/36解一阶常微分方程欧拉法解一阶常微分方程欧拉法Range-Kutta公式公式一阶微分方程组与一阶微分方程组与二阶方程二阶方程常微分方程边值问题常微分方程边值问题线性多步法简介线性多步法简介数值分析 第八章第八章2/36一阶常微分方程初值问题一阶常微分方程初值问题:000)(),(yxyxxyxfy数值方法数值方法取离散点取离散点:x0 x1 x2 xn 其中其中,y=y(x)是未知函数是未知函数,f(x,y)是已知函数是已知函数.初值初值 y0 是已知数据。是已知数据。求未知函数求未知函数 y(x1),y(x2),y(xn),的近似值的近似值y1,y2,y3,yn)(nxyny称为常微分
2、方程的数值解。称为常微分方程的数值解。这里这里 是是 的数值逼近的数值逼近.3/36例例1.常微分方程与向量场常微分方程与向量场平面区域平面区域:0 x 1.5,0 y 1.5,)(2),(yxyxf ),(yxfy 0)0(yy ijijjikyxf tan),(斜率斜率方向余弦方向余弦1,11sin,cos22ijijijijijkkk 4/36 000)(),(yxyxxyxfy取步长取步长 h,记点记点:),(1nnnnyxfhyy 在第在第 n 个点处个点处),2,1,0(,0 nnhxxn)()(1)(010 xyxyhxy 一阶向前差商一阶向前差商:),(0001yxfhyy )
3、,(0001yxhfyy ),(1nnnnyxhfyy 欧拉公式欧拉公式),2,1,0(n初值问题初值问题:5/36例例2.Logistic模型模型)1(yyy )exp(911)(xxy 解析解解析解:1.0)0(y)1(1nnnnyhyyy 欧拉公式欧拉公式),2,1,0(n6/3601234500.20.40.60.81取点取点(xn,yn),(n=0,1,2,)欧拉公式解的几何解释欧拉公式解的几何解释:取取 x=xn+1 得得:yn+1=yn+h f(xn,yn)点斜式直线方程点斜式直线方程:y=yn+(x xn)f(xn,yn),取点取点(xn,yn)(xn+1,yn+1)7/36y
4、=f(x,y)1)(,()()(1nnxxnndxxyxfxyxy),(),(2)(,(111 nnnnxxyxfyxfhdxxyxfnn梯形公式梯形公式:),(1nnnnyxhfyy 左矩形公式左矩形公式),()(,(1nnxxyxhfdxxyxfnn 用数值积分方法离散常微分方程用数值积分方法离散常微分方程 11)(,()(nnnnxxxxdxxyxfdxxy8/36),(nnnnyxhfyy1),(),(2111 nnnnnnyxfyxfhyy),(),(2111 nnnnnnyxfyxfhyy预预-校方法又称为修正的校方法又称为修正的Euler法法,算法如下算法如下 k1=f(xn,y
5、n),k2=f(xn+1,yn+h k1),2211kkhyynn 由梯形公式推出的预由梯形公式推出的预-校方法校方法:9/36设设 yn=y(xn),称称 Rn+1=y(xn+1)-yn+1为为局部截断误差局部截断误差.)(2)()()()()(2111 yxxxyxxxyxynnnnnnn )(2),()()(21 yhyxhfxyxynnnn 即即由泰勒公式由泰勒公式Euler公式公式:yn+1=yn+hf(xn,yn)的局部截断误差的局部截断误差y(xn+1)yn+1=y(xn)yn+O(h2)=O(h2)Euler公式的局部截断误差记为公式的局部截断误差记为:O(h2)称称Euler
6、公式具有公式具有1阶精度。阶精度。10/36若局部截断误差为若局部截断误差为:O(h p+1)则称显式单步法具有则称显式单步法具有 p 阶精度阶精度。例例 3 证明修正的证明修正的Euler法具有法具有2阶精度阶精度),(),(2111 nnnnnnyxfyxfhyy),(1nnnnyxhfyy 证证 由预测公式由预测公式),(,(),(11nnnnnnyxhfyhxfyxf )(),(),(),(),(2hOyxfyxhfyxhfyxfnnynnnnxnn ),()(nnnyxfxy ),(),(),()(nnnnynnxnyxfyxfyxfxy 11/36)()(2)()()(321hOx
7、yhxyhxyxynnnn 由由Taylor 级数级数)()()(),(211hOxyhxyyxfnnnn ),(),(2111 nnnnnnyxfyxfhyy)()()()(231hOxyhxyxyhyynnnnn )()(311hOyxynn nnyxy)(设设局部截断误差局部截断误差:故修正的故修正的Euler法法具有具有2阶精度。阶精度。12/36三阶三阶Range-Kutta公式一般形式公式一般形式yn+1=yn+hk1+4k2+k3/6k1=f(xn,yn),k2=f(xn+0.5h,yn+0.5hk1)k3=f(xn+h,yn hk1+2hk2)四阶四阶Range-Kutta公式
8、一般形式公式一般形式yn+1=yn+hk1+2k2+2k3+k4/6k1=f(xn,yn),k2=f(xn+0.5h,yn+0.5hk1)k3=f(xn+0.5h,yn+0.5hk2),k4=f(xn+h,yn+hk3)13/36 1)0(20,2yxxyydxdyxexxy 211)(例例4数值实验数值实验:几种不同求数值解公式的误差比较几种不同求数值解公式的误差比较 n 10 20 30 40 h 0.2 0.1 0.0667 0.05RK4 6.862e-005 3.747e-006 7.071e-007 2.186e-007RK3 0.0012 1.529e-004 4.517e-00
9、5 1.906e-005RK2 0.0123 0.0026 0.0011 5.9612e-004Euler 0.1059 0.0521 0.0342 0.025614/36MATLAB求解常微分方程初值问题命令求解常微分方程初值问题命令：(1)定义一阶微分方程的右端函数定义一阶微分方程的右端函数；(2)用用MATLAB命令命令ode23()求数值解。求数值解。使用格式使用格式：T,Y=ode23(F,Tspan,y0)其中其中，Tspan=t0，tN是常微分方程的求解区域是常微分方程的求解区域，y0是解的初值是解的初值15/36实验例题实验例题1 蛇形曲线的常微分方程初值问题蛇形曲线的常微分方
10、程初值问题 22211yxy 0)0(yMATLAB数值求解命令数值求解命令F=inline(1./(1+x.2)-2*y.2)；ode23(F,0,6,0)输出结果为图形输出结果为图形 012345600.10.20.30.40.5T,y=ode23(f,0,6,0)将得到将得到自变量和函数的离散数据自变量和函数的离散数据 16/36MATLAB解常微分方程初值问题命令解常微分方程初值问题命令 1)0(20,2yxxyydxdy数值求解命令数值求解命令:x,y=ode23(f,a,b,y0)f=inline(y-x.*y.2);x,y=ode23(f,0,2,1)符号求解命令符号求解命令:d
11、solve(eqn1,.)syms x ydsolve(Dy=y-x*y2,y(0)=1,x)ans=1/(x-1+2*exp(-x)xexxy 211)(解析解解析解:17/36 ),(),(21yxtfdtdyyxtfdtdx两个未知函数的一阶常微分方程组两个未知函数的一阶常微分方程组00)(xtx 00)(yty 常微分方程组的向量形式常微分方程组的向量形式),(YtFdtdY 00)(YtY 记记 yxY ),(),(),(21YtfYtfYtF18/36欧拉公式欧拉公式:TyxY,000),(1nnnnYthFYY ),2,1,0(n ),(),(2111nnnnnnnnnnyxtf
12、yxtfhyxyx ),(),(2111nnnnnnnnnnyxthfyyyxthfxx19/366/2243211KKKKhYYnn ),()5.0,5.0()5.0,5.0(),(3423121hKYhtFKhKYhtFKhKYhtfKYtFKnnnnnnnn 修改的欧拉公式修改的欧拉公式:2/211KKhYYnn )5.0,5.0(),(121hKYhtFKYtFKnnnn (n=0,1,N)经典龙格经典龙格-库塔公式库塔公式:20/36捕食者与被捕食者问题捕食者与被捕食者问题 海岛上有狐狸和野兔海岛上有狐狸和野兔,当野兔数量增多时，狐狸捕食当野兔数量增多时，狐狸捕食野兔导致狐群数量增长
13、野兔导致狐群数量增长;大量兔子被捕食使狐群进入大量兔子被捕食使狐群进入饥饿状态其数量下降饥饿状态其数量下降;狐群数量下降导致兔子被捕食狐群数量下降导致兔子被捕食机会减少机会减少,兔群数量回升。微分方程模型如下兔群数量回升。微分方程模型如下 xyydtdyxyxdtdx02.001.0计算计算 x(t)，y(t)当当t0，20时的数据。绘图并分时的数据。绘图并分析捕食者和被捕食者的数量变化规律。析捕食者和被捕食者的数量变化规律。x(0)=100y(0)=20 21/36 xyyyxyxx02.001.0 xyyxyxyxtfyxtf02.001.0),(),(21 xyyxyxppyx02.00
14、1.0平面向量场平面向量场:向量场中向量场中过点过点:(100,20)的轨线的轨线22/36MATLAB命令求解命令求解:Y0=100,20;t,Y=ode23(fox,0,20,Y0);x=Y(:,1);y=Y(:,2);figure(1),plot(t,x,b,t,y,r)figure(2),plot(x,y)function z=fox(t,y)z(1,:)=y(1)-0.01*y(1).*y(2);z(2,:)=-y(2)+0.02*y(1).*y(2);-y1 -y2 y1y2 相位图相位图定义方程右端函数定义方程右端函数05101520010020030040023/36“蝴蝶效应
15、蝴蝶效应”来源于洛伦兹一次讲演来源于洛伦兹一次讲演。模型如下模型如下 zyxydtdzzydtdyyzxdtdx )(求微分方程数值解求微分方程数值解,绘出解函数曲线绘出解函数曲线 取取 =8/3，=10，=28。x(0)=0，y(0)=0，z(0)=0.01。t0，80，zyxyzyxyxzyxtF2810103/8),(微分方程右端函数微分方程右端函数:24/36记向量记向量 y1，y2，y3=x，y，z，创建函数文件创建函数文件function z=flo(t,y)z(1,:)=-8*y(1)/3+y(2).*y(3);z(2,:)=-10*(y(2)-y(3);z(3,:)=-y(1)
16、.*y(2)+28*y(2)-y(3);用用MATLAB命令求解并绘出命令求解并绘出Y-X平面的投影图平面的投影图 P0=0;0;0.01;T,P=ode23(flo,0,80,P0);figure(1),plot(P(:,2),P(:,1)figure(2),comet3(P(:,1),P(:,2),P(:,3)25/36分量分量 x 的误差的误差分量分量 y 的误差的误差分量分量 z 的误差的误差26/36二阶常微分方程问题二阶常微分方程问题02 xkx )(2tfxkxpx 02 xkxpx(简谐振动简谐振动)(衰减振动衰减振动)(受迫振动受迫振动)n 阶勒让德方程阶勒让德方程0)(22
17、2 ynxyxyxn 阶贝塞尔方程阶贝塞尔方程0)1(2)1(2 ynnyxyx27/36 0000)(,)(),(mxyyxyyyxfy令令 00)(),(YxYYxFdxdY ),(21221yyxfyyy一阶常微分方程组一阶常微分方程组:初值条件初值条件:002001)(,)(mxyyxy ),()(),(212yyxfxyYxF 000myY常微分方程组常微分方程组)()(1xyxy)()(2xyxy )()()(21xyxyxY28/36例例3.单摆的数学模型单摆的数学模型 sina 其中其中,a=g/L初值条件初值条件:(0)=0.4，(0)=0 第一步第一步:转化为一阶方程组转化
18、为一阶方程组 1221sin yayyy令令:y1=,y2=初值条件初值条件:y1(0)=0.4,y2(0)=0第二步第二步:求解方程组求解方程组function f=dan(x,y)f(1,:)=y(2);f(2,:)=-9.8*sin(y(1)/3.2;L=3.2ode23(dan,0,2,0.4,0);29/36t,thata=ode23(dan,0,2.755,0.6,0);R=3.2;n=length(t);alpha=thata(:,1);x=R*sin(alpha);y=R*cos(alpha);X=0,0;Y=0,-3.5;for k=1:n xk=x(1:k);yk=y(1:
19、k);Xk=x(k);Yk=y(k);plot(xk,-yk,.-r,Xk,-Yk,o,0,Xk,0,-Yk),axis(-2.5,2.5,-3.5,0)pause(.5)end-2-1012-3.5-3-2.5-2-1.5-1-0.50单摆的动态模拟程序单摆的动态模拟程序30/36人造卫星的轨道模型人造卫星的轨道模型万有引力定律万有引力定律 22|yxMmGF 2/32222)(yxGMxdtxd 2/32222)(yxGMydtyd 222222,yxyyxxyxMmGF地球引力参数地球引力参数:GM=3.986005105（km3/s2）常微分方程常微分方程31/36例例4 求解求解边值
20、问题边值问题的数值方法算例的数值方法算例 0)1(,0)0(10,0yyxxyy解解:取正整数取正整数n,令令h=1/(n+1),xj=jh,(j=0,1,n+1).将常微分方程离散化将常微分方程离散化 02211 jjjjjxyhyyy整理整理,得得:yj-1+(2 h2)yj yj+1=xj h2 (j=1,2,n)y0=0,yn=01.打靶法打靶法;2.高斯消元法高斯消元法32/36yj-1+(2 h2)yj yj+1=xj h2 (j=1,2,n)nnhhhhA 22222112112112三对角方程组三对角方程组 AY=F y(xn);o yn00.20.40.60.8100.020
21、.040.060.0833/36线性多步法线性多步法kiinikiiniknfhyy010 (n=0,1,)其中其中,xn+i=x0+(n+i)h,fn+i=f(xn+i,yn+i)局部载断误差局部载断误差kiinikiiniknknxyhxyxyT010)()()(Adamas显格式显格式:yn+2=yn+1+h(3fn+1-fn)/2 yn+3=yn+2+h(23fn+2-16fn+1+5fn)/1234/36y=f(x,y)21)(,()()(12nnxxnndxxyxfxyxy 2121)(,()(nnnnxxxxdxxyxfdxxy在区间在区间xn,xn+1上插值上插值f(x)(xn
22、+1x)fn+(xxn)fn+1/h 21)()(111nnxxnnnndxfxxfxxh2/323211122nnnnffhfhfhh 二阶二阶Adamas显格式显格式:yn+2=yn+1+h(3fn+1-fn)/235/36思考题与练习题思考题与练习题1.例举几种求解常微分方程初值问题的数值法例举几种求解常微分方程初值问题的数值法,简述它们各自特点？简述它们各自特点？2.求常微分方程初值问题数值法的局部截断误差求常微分方程初值问题数值法的局部截断误差 与精度有何关系？与精度有何关系？3.如何利用一阶常微分方程组的数值方法求解高如何利用一阶常微分方程组的数值方法求解高 阶常微分方程初值问题阶
23、常微分方程初值问题?4.4.试列举出三个数学物理或工程中常见的二阶常试列举出三个数学物理或工程中常见的二阶常 微分方程微分方程(或方程组或方程组),),简述已有的解析方法并简述已有的解析方法并 与数值法比较与数值法比较.36/36 xdttexI022cos1)(5.将勒让德第二类椭圆积分将勒让德第二类椭圆积分转换为常微分方程初值问题转换为常微分方程初值问题,并比较常微分方程数值并比较常微分方程数值解法与数值积分方法相异与相同之处解法与数值积分方法相异与相同之处.000)(),(yxyxxxfy 000)(),(yxyxxyfy4.用四阶用四阶Range-Kutta公式求解下列问题有何特公式求解下列问题有何特 殊之处殊之处
