数值分析常微分方程边值问题

数值分析常微分方程边值问题第一页，共十六页，编辑于2023年，星期三初值问题欧拉公式:yn+1=yn+hf(xn,yn)修改的欧拉公式(Heun方法):k1=f(xn,yn),k2=f(xn+1,yn+hk1)(n=0,1,2,·······,N)(n=0,1,2,·······,N)2/16第二页，共十六页，编辑于2023年，星期三二阶Range-Kutta公式一般形式yn+1=yn+h[c1k1+c2k2]k1=f(xn,yn)k2=f(xn+λ2h,yn+μ21hk1)(n=0,1,····,N)c1=0.5,c2=0.5,λ2=1,

μ21=1c1=0,c2=1,λ2=0.5,

μ21=0.5

yn+1=yn+0.5h[k1+k2]k1=f(xn,yn)k2=f(xn+h,yn+hk1)

yn+1=yn+h

k2k1=f(xn,yn)k2=f(xn+.5h,yn+.5hk1)3/16第三页，共十六页，编辑于2023年，星期三三阶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/16第四页，共十六页，编辑于2023年，星期三例1数值实验:几种不同求数值解公式的误差比较

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.02565/16第五页，共十六页，编辑于2023年，星期三n=input('inputn:=');f=inline('y-x.*y.^2');h=2/n;x=h:h:2;y0=1;k1=f(0,y0);k2=f(.5*h,y0+0.5*h*k1);k3=f(.5*h,y0+.5*h*k2);k4=f(h,y0+h*k3);y(1)=y0+h*(k1+2*k2+2*k3+k4)/6;fork=1:n-1xk=x(k);yk=y(k);k1=f(xk,yk);k2=f(xk+h/2,yk+0.5*h*k1);k3=f(xk+h/2,yk+0.5*h*k2);k4=f(xk+h,yk+h*k3);y(k+1)=yk+h*(k1+2*k2+2*k3+k4)/6;endt=-1:h/2:5;yy=1./(t-1+2*exp(-t));plot(0,1,'ro',x,y,'ro',t,yy)y1=1./(x-1+2*exp(-x));error=max(abs(y-y1))inputn:=20error=3.7475e-0066/16第六页，共十六页，编辑于2023年，星期三例2.捕食者与被捕食者模型0≤x≤15

y1(0)=20,y2(0)=20记7/16第七页，共十六页，编辑于2023年，星期三x(1)=20;y(1)=20;forn=1:length(t)-1tn=t(n);xn=x(n);yn=y(n);k11=f1(xn,yn);k21=f2(xn,yn);k12=f1(xn+.5*h*k11,yn+.5*h*k21);k22=f2(xn+.5*h*k11,yn+.5*h*k21);k13=f1(xn+.5*h*k12,yn+.5*h*k22);k23=f2(xn+.5*h*k12,yn+.5*h*k22);k14=f1(xn+h*k13,yn+h*k23);k24=f2(xn+h*k13,yn+h*k23);x(n+1)=xn+h*(k11+2*k12+2*k13+k14)/6;y(n+1)=yn+h*(k21+2*k22+2*k23+k24)/6;endf1=inline('x-0.01*x*y');f2=inline('-y+0.02*x*y');h=.15;t=0:h:15;----y1----y2

y1—y2

相位图8/16第八页，共十六页，编辑于2023年，星期三例1.单摆的数学模型其中,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]);9/16第九页，共十六页，编辑于2023年，星期三求解边值问题的数值方法算例解:取正整数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.打靶法;3.追赶法10/16第十页，共十六页，编辑于2023年，星期三精确解symsxydsolve('D2y+y=-x','y(0)=0','y(1)=0','x')符号系统计算结果为：ans=-x+1/sin(1)*sin(x)

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

(j=1,2,···,n)方程组AY=f

系数矩阵:11/16第十一页，共十六页，编辑于2023年，星期三n=input('inputn:=');y=inline('sin(x)/sin(1)-x');h=1/(n+1);x=0:h:1;hh=h*h;r=2-hh;u(1)=0;u(n+2)=0;er=1;k=0;whileer>0.0001er=0;k=k+1;forj=2:n+1s=(hh*x(j)+u(j-1)+u(j+1))/r;er=max(er,abs(s-u(j)));u(j)=s;endendv=y(x);plot(x,v,x,u,'o')er=max(abs(u-v))—----y(xn);O----yn迭代格式:yj

=[xjh2+yj-1+yj+1]/(2–h2)12/16输入:n=10计算结果：

er=0.0012第十二页，共十六页，编辑于2023年，星期三追赶法解三对角方程组n=input('inputn:=');h=1/(n+1);hh=h*h;a=-1;b=2-hh;x=0:h:1;f=hh*x;f(1)=[];f(n+1)=[];f(1)=f(1)/b;d=b;fori=2:nc(i-1)=a/d;d=b-a*c(i-1);f(i)=(f(i)-a*f(i-1))/d;endfori=n-1:-1:1f(i)=f(i)-c(i)*f(i+1);endy=sin(x)/sin(1)-x;yk=[0f0];plot(x,y,x,yk,'o')error=max(abs(yk-y))输入:n=10计算结果：error=5.4566e-005O--------数值解;—--------精确解13/16第十三页，共十六页，编辑于2023年，星期三边值问题差分方法:u(a)=

,u(b)=

一般二阶方程:指数变换

1.高精度差分格式;2.组合外推方法;3.大型稀疏方程组解算子参考:

ElsevierSDOS(电子期刊)Appliedmathematicsandcomputation14/16第十四页，共十六页，编辑于2023年，星期三线性多步法(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+