数值分析第八章New

解一阶常微分方程欧拉法Range-Kutta公式一阶微分方程组与二阶方程常微分方程边值问题线性多步法简介《数值分析》第八章一阶常微分方程初值问题:数值方法——取离散点:

x0<x1<x2<···<xn

······

其中,

y=y(x)是未知函数,f(x,y)是已知函数.初值y0是已知数据。求未知函数y(x1),y(x2),····,y(xn),······的近似值y1,y2,y3,·····,yn·······称为常微分方程的数值解。这里是的数值逼近.例1.常微分方程与向量场平面区域:0≤x≤1.5,0≤y≤1.5,斜率方向余弦取步长h,记点:在第n个点处一阶向前差商:欧拉公式初值问题:例2.Logistic模型

解析解:欧拉公式取点

(xn

,yn),(n=0,1,2,···

)欧拉公式解的几何解释:取

x=xn+1得:yn+1=yn+hf(xn,yn)点斜式直线方程:y=yn+(x–xn)f(xn,yn),取点

(xn

,yn)(xn+1,yn+1)y’=f(x,y)梯形公式:

左矩形公式用数值积分方法离散常微分方程预-校方法又称为修正的Euler法,算法如下k1=f(xn,yn),

k2=f(xn+1,yn+hk1),由梯形公式推出的预-校方法:设

yn=y(xn),称Rn+1=y(xn+1)-yn+1为局部截断误差.即由泰勒公式Euler公式:yn+1=yn+hf(xn,yn)的局部截断误差y(xn+1)–yn+1=y(xn)–yn+O(h2)=O(h2)Euler公式的局部截断误差记为:

O(h2)称Euler公式具有1阶精度。若局部截断误差为:O(hp

+1)

则称显式单步法具有p阶精度。例3证明修正的Euler法具有2阶精度证由预测公式由Taylor级数设局部截断误差:故修正的Euler法具有2阶精度。三阶Range-Kutta公式一般形式yn+1=yn+h[k1+4k2+k3]/6k1=f(xn,yn),k2=f(xn+0.5h,yn+0.5hk1)k3=f(xn+h,yn–hk1+2hk2)四阶Range-Kutta公式一般形式yn+1=yn+h[k1+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)例4数值实验:几种不同求数值解公式的误差比较n10203040h0.20.10.06670.05RK46.862e-0053.747e-0067.071e-0072.186e-007RK30.00121.529e-0044.517e-0051.906e-005RK20.01230.00260.00115.9612e-004Euler0.10590.05210.03420.0256MATLAB求解常微分方程初值问题命令：

(1)定义一阶微分方程的右端函数；(2)用MATLAB命令ode23()求数值解。使用格式：[T,Y]=ode23('F',Tspan,y0)其中，Tspan=[t0，tN]是常微分方程的求解区域，y0是解的初值实验例题1蛇形曲线的常微分方程初值问题

MATLAB数值求解命令F=inline('1./(1+x.^2)-2*y.^2')；ode23(F,[0,6],0)输出结果为图形

[T,y]=ode23(f,[0,6],0)将得到自变量和函数的离散数据

MATLAB解常微分方程初值问题命令数值求解命令:[x,y]=ode23('f',[a,b],y0)f=inline('y-x.*y.^2');[x,y]=ode23(f,[0,2],1)符号求解命令:dsolve('eqn1',...)symsxydsolve('Dy=y-x*y^2','y(0)=1','x')ans=1/(x-1+2*exp(-x))解析解:两个未知函数的一阶常微分方程组常微分方程组的向量形式记欧拉公式:修改的欧拉公式:(n=0,1,·······,N)经典龙格-库塔公式:捕食者与被捕食者问题

海岛上有狐狸和野兔,当野兔数量增多时，狐狸捕食野兔导致狐群数量增长;大量兔子被捕食使狐群进入饥饿状态其数量下降;狐群数量下降导致兔子被捕食机会减少,兔群数量回升。微分方程模型如下计算x(t)，y(t)当t∈[0，20]时的数据。绘图并分析捕食者和被捕食者的数量变化规律。x(0)=100y(0)=20

平面向量场:——向量场中过点:(100,20)

的轨线MATLAB命令求解: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)functionz=fox(t,y)z(1,:)=y(1)-0.01*y(1).*y(2);z(2,:)=-y(2)+0.02*y(1).*y(2);----y1----y2y1—y2

相位图定义方程右端函数“蝴蝶效应”来源于洛伦兹一次讲演。模型如下求微分方程数值解,

绘出解函数曲线取=8/3，=10，=28。x(0)=0，y(0)=0，z(0)=0.01。t∈[0，80]，微分方程右端函数:记向量[y1，y2，y3]=[x，y，z]，创建函数文件functionz=flo(t,y)z(1,:)=-8*y(1)/3+y(2).*y(3);z(2,:)=-10*(y(2)-y(3));z(3,:)=-y(1).*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))分量x

的误差分量y

的误差分量z

的误差二阶常微分方程问题(简谐振动)(衰减振动)(受迫振动)n

阶勒让德方程n阶贝塞尔方程令一阶常微分方程组:初值条件:常微分方程组例3.单摆的数学模型其中,a=g/L初值条件:(0)=0.4，

’(0)=0第一步:转化为一阶方程组令:y1=,y2=

’

初值条件:y1(0)=0.4,y2(0)=0第二步:求解方程组functionf=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]);[t,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];fork=1:nxk=x(1:k);yk=y(1: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单摆的动态模拟程序人造卫星的轨道模型万有引力定律

地球引力参数:

GM=3.986005×105（km3/s2）

常微分方程例4求解边值问题的数值方法算例解:取正整数n,令h=1/(n+1),xj=jh,(j=0,1,···,n+1

).

将常微分方程离散化

整理,得:

–yj-1+(2–h2)yj–yj+1=xjh2

(j=1,2,···,n)y0=0,yn=01.打靶法;2.高斯消元法–yj-1+(2–h2)yj–yj+1=xjh2

(j=1,2,···,n)三对角方程组AY=F—

y(xn);o

yn线性多步法(n=0,1,···)其中,xn+i=x0+(n+i)h,fn+i=f(xn+i,yn+i)局部载断误差Adamas显格式:yn+2=yn+1+h(3fn+1-fn)/2yn+3=yn+2+h(23fn+2-16fn+1+5fn)/12y’=f(x,y)在区间[xn,xn+1]上插值f(x)≈[(xn+1－x)fn+(x－xn)fn+1]/h二阶Adamas显格式:yn+2=yn+1+h(3fn+1-fn)/2思考题与练习题1.例举几种求