數(shù)值分析常微分方程邊值問題

數(shù)值分析常微分方程邊值問題第1頁，課件共16頁，創(chuàng)作于2023年2月初值問題歐拉公式: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第2頁，課件共16頁，創(chuàng)作于2023年2月二階Range-Kutta公式一般形式y(tǒng)n+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第3頁，課件共16頁，創(chuàng)作于2023年2月三階Range-Kutta公式一般形式y(tǒng)n+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公式一般形式y(tǒng)n+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第4頁，課件共16頁，創(chuàng)作于2023年2月例1數(shù)值實(shí)驗(yàn):幾種不同求數(shù)值解公式的誤差比較

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第5頁，課件共16頁，創(chuàng)作于2023年2月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第6頁，課件共16頁，創(chuàng)作于2023年2月例2.捕食者與被捕食者模型0≤x≤15

y1(0)=20,y2(0)=20記7/16第7頁，課件共16頁，創(chuàng)作于2023年2月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第8頁，課件共16頁，創(chuàng)作于2023年2月例1.單擺的數(shù)學(xué)模型其中,a=g/L初值條件:(0)=0.4，

’(0)=0第一步:轉(zhuǎn)化為一階方程組令: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第9頁，課件共16頁，創(chuàng)作于2023年2月求解邊值問題的數(shù)值方法算例解:取正整數(shù)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第10頁，課件共16頁，創(chuàng)作于2023年2月精確解symsxydsolve('D2y+y=-x','y(0)=0','y(1)=0','x')符號(hào)系統(tǒng)計(jì)算結(jié)果為：ans=-x+1/sin(1)*sin(x)

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

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

系數(shù)矩陣:11/16第11頁，課件共16頁，創(chuàng)作于2023年2月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計(jì)算結(jié)果：

er=0.0012第12頁，課件共16頁，創(chuàng)作于2023年2月追趕法解三對(duì)角方程組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計(jì)算結(jié)果：error=5.4566e-005O--------數(shù)值解;—--------精確解13/16第13頁，課件共16頁，創(chuàng)作于2023年2月邊值問題差分方法:u(a)=

,u(b)=

一般二階方程:指數(shù)變換

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

ElsevierSDOS(電子期刊)Appliedmathematicsandcomputation14/16第14頁，課件共16頁，創(chuàng)作于2023年2月線性多步法(n=0,1,···)其中,xn+i=x0+(n+i)h,fn+i=f(xn+i,yn+i)局部載斷誤差A(yù)damas顯格式:yn+2=yn+1+h(3fn+1-fn)/2yn+3=yn+2+h(23fn+