数值分析常微分方程边值问题幻灯片.ppt

《数值分析常微分方程边值问题幻灯片.ppt》由会员分享，可在线阅读，更多相关《数值分析常微分方程边值问题幻灯片.ppt（16页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

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

2、.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第3页，共16页，编辑于2022年，星期六三阶三阶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公式一般形式公式一般形式yn+1=yn+hk1+2k2+2k3+k4/6k1=f(xn,yn),k2=f(xn+0.5h,y

3、n+0.5hk1)k3=f(xn+0.5h,yn+0.5hk2),k4=f(xn+h,yn+hk3)4/16第4页，共16页，编辑于2022年，星期六例例1数值实验数值实验:几种不同求数值解公式的误差比较几种不同求数值解公式的误差比较 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-005 1.906e-005RK2 0.0123 0.0026 0.0011 5.9612e-004Euler 0.1059 0.0521 0

4、.0342 0.02565/16第5页，共16页，编辑于2022年，星期六n=input(input n:=);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;for k=1:n-1 xk=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

5、*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)input n:=20error=3.7475e-0066/16第6页，共16页，编辑于2022年，星期六例例2.捕食者与被捕食者模型捕食者与被捕食者模型0 x 15 y1(0)=20,y2(0)=20记记7/16第7页，共16页，编辑于2022年，星期六x(1)=20;y(1)=20;for n=1:length(t)-1

6、tn=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

7、)/6;endf1=inline(x-0.01*x*y);f2=inline(-y+0.02*x*y);h=.15;t=0:h:15;-y1 -y2 y1y2 相位图相位图8/16第8页，共16页，编辑于2022年，星期六例例1.单摆的数学模型单摆的数学模型其中其中,a=g/L初值条件初值条件:(0)=0.4，(0)=0 第一步第一步:转化为一阶方程组转化为一阶方程组令令: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.2o

8、de23(dan,0,2,0.4,0);9/16第9页，共16页，编辑于2022年，星期六求解求解边值问题边值问题的数值方法算例的数值方法算例解解:取正整数取正整数n,令令h=1/(n+1),xj=jh,(j=0,1,n+1).将常微分方程离散化将常微分方程离散化 整理整理,得得:yj-1+(2 h2)yj yj+1=xj h2 (j=1,2,n)y0=0,yn=01.高斯高斯-赛德尔迭代赛德尔迭代;2.打靶法打靶法;3.追赶法追赶法10/16第10页，共16页，编辑于2022年，星期六精确解精确解syms x ydsolve(D2y+y=-x,y(0)=0,y(1)=0,x)符号系统计算结果

9、为：符号系统计算结果为：ans=-x+1/sin(1)*sin(x)yj-1+(2 h2)yj yj+1=xj h2 (j=1,2,n)方程组方程组AY=f 系数矩阵系数矩阵:11/16第11页，共16页，编辑于2022年，星期六n=input(input n:=);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;while er0.0001 er=0;k=k+1;for j=2:n+1 s=(hh*x(j)+u(j-1)+u(j+1)/r;er=max(er,abs(s-u(

10、j);u(j)=s;endendv=y(x);plot(x,v,x,u,o)er=max(abs(u-v)-y(xn);O -yn迭代格式迭代格式:yj=xj h2+yj-1+yj+1/(2 h2)12/16输入输入:n=10计算结果：计算结果：er=0.0012第12页，共16页，编辑于2022年，星期六追赶法解三对角方程组追赶法解三对角方程组n=input(input n:=);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;for i=2:n c(i-1)=a/d;d=b-a*c(i-1)

11、;f(i)=(f(i)-a*f(i-1)/d;endfor i=n-1:-1:1 f(i)=f(i)-c(i)*f(i+1);endy=sin(x)/sin(1)-x;yk=0 f 0;plot(x,y,x,yk,o)error=max(abs(yk-y)输入输入:n=10计算结果：计算结果：error=5.4566e-005O -数值解数值解;-精确解精确解13/16第13页，共16页，编辑于2022年，星期六边值问题差分方法边值问题差分方法:u(a)=,u(b)=一般二阶方程一般二阶方程:指数变换指数变换 1.高精度差分格式高精度差分格式;2.组合外推方法组合外推方法;3.大型稀疏方程组解

12、算子大型稀疏方程组解算子参考参考:Elsevier SDOS(电子期刊电子期刊)Applied mathematics and computation 14/16第14页，共16页，编辑于2022年，星期六线性多步法线性多步法(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)/2 yn+3=yn+2+h(23fn+2-16fn+1+5fn)/1215/16第15页，共16页，编辑于2022年，星期六y=f(x,y)在区间在区间xn,xn+1上插值上插值f(x)(xn+1x)fn+(xxn)fn+1/h二阶二阶Adamas显格式显格式:yn+2=yn+1+h(3fn+1-fn)/216/16第16页，共16页，编辑于2022年，星期六