試驗四常微分方程初值問題數(shù)值解法講解

1、實 驗 報 告課程名稱 數(shù)值分析實驗項目 常微分方程問題初值問題數(shù)值解法. 實驗?zāi)康?、理解如何在計算機上實現(xiàn)用 Euler 法、改進(jìn) Euler 法、Runge Kutta 算法求一 階常微分方程初值問題y （x） f （x, y）,x a,b y（a） y1的數(shù)值解。2、利用圖形直觀分析近似解和準(zhǔn)確解之間的誤差。3、學(xué)會 Matlab 提供的 ode45 函數(shù)求解微分方程初值問題。二、實驗要求（ 1） 按照題目要求完成實驗內(nèi)容；（ 2） 寫出相應(yīng)的 Matlab 程序；（ 3） 給出實驗結(jié)果 （可以用表格展示實驗結(jié)果 ）；（ 4） 分析和討論實驗結(jié)果并提出可能的優(yōu)化實驗。（ 5） 寫出實驗

2、報告。三、實驗步驟1、用編好的 Euler法、改進(jìn) Euler法計算書本 P167 的例 1、P171例題 3（1）取 h 0.1，求解初值問題y （x） y 2x,x 0,1,yy（0） 1（2）取 h 0.1，求解初值問題y （x） y x 1,x 0,0.5,y（0） 12、用 RungeKutta算法計算 P178例題、 P285實驗任務(wù)（ 2）1)取 h 0.1，求解初值問題y (x) y2,x 0,0.5,y(0) 1(2)求初值問題y (x) 1 ( y x2 4x 1), x 0,0.5,2y(0) 0x 的解 y(x) 在 xi ih(h 0.05) 處的近似值 yi ，并與

3、問題的解析解 y(x) e 2 x2 1 相比較。3、用 Matlab 繪圖函數(shù) plot(x,y) 繪制 P285實驗任務(wù)( 2)的精確解和近似解的圖 形。4、使用 matlab中的 ode45求解 P285實驗任務(wù)( 2)，并繪圖。四、實驗結(jié)果1、Euler 算法程序、改進(jìn) Euler 算法程序；Euler算法程序： function x,y=euler_f(ydot_fun, x0, y0, h, N) % Euler (向前)公式，其中%ydot_fun-%x0, y0 -%h -%N -%x -Xn%y -Yn一階微分方程的函數(shù)初始條件區(qū)間步長區(qū)間的個數(shù)構(gòu)成的向量構(gòu)成的向量x=zer

4、os(1,N+1); y=zeros(1,N+1); x(1)=x0; y(1)=y0; for n=1:Nx(n+1)=x(n)+h;y(n+1)=y(n)+h*feval(ydot_fun, x(n), y(n); end改進(jìn) Euler 算法程序：function x,y=euler_r(ydot_fun, x0, y0, h, N) 改進(jìn) Euler 公式，其中 ydot_fun - x0, y0 -h%一階微分方程的函數(shù)初始條件區(qū)間步長 區(qū)間的個數(shù) Xn 構(gòu)成的向量% y - Yn 構(gòu)成的向量x=zeros(1,N+1); y=zeros(1,N+1); x(1)=x0; y(1)=

5、y0; for n=1:Nx(n+1)=x(n)+h; ybar=y(n)+h*feval(ydot_fun, x(n), y(n); y(n+1)=y(n)+h/2*(feval(ydot_fun, x(n), y(n)+feval(ydot_fun, x(n+1), ybar); end2、用 Euler算法程序、改進(jìn) Euler算法求解 P167例題 1的運行結(jié)果；（1.）Euler 算法程序：Columns 1 through 80 0.1000 0.2000 0.3000 0.4000 0.50000.70000.6000Columns 9 through 110.8000 0.90

6、00 1.00001.50900.50001.4164Columns 1 through 81.0000 1.1000 1.1918 1.2774 1.3582 1.4351 1.5803Columns 9 through 111.6498 1.7178 1.7848（2）改進(jìn) Euler 算法：x =Columns 1 through 70 0.1000 0.2000 0.3000 0.4000 0.6000Columns 8 through 110.7000 0.8000 0.9000 1.0000y =Columns 1 through 71.0000 1.0959 1.1841 1.2

7、662 1.3434 1.4860Columns 8 through 111.5525 1.6165 1.6782 1.73793、RungeKutta 算法程序；function x,y=Runge_Kutta4(ydot_fun, x0, y0, h, N) % 標(biāo)準(zhǔn)四階 Runge_Kutta 公式，其中 %ydot_fun - 一階微分方程的函數(shù)%x0, y0 - 初始條件%h- 區(qū)間步長%N- 區(qū)間的個數(shù)% x - Xn 構(gòu)成的向量% y - Yn 構(gòu)成的向量x=zeros(1,N+1); y=zeros(1,N+1); x(1)=x0; y(1)=y0; for n=1:Nx(n+

8、1)=x(n)+h;k1=h*feval(ydot_fun, x(n), y(n);k2=h*feval(ydot_fun, x(n)+1/2*h, y(n)+1/2*k1);k3=h*feval(ydot_fun, x(n)+1/2*h, y(n)+1/2*k2);k4=h*feval(ydot_fun, x(n)+h, y(n)+k3); y(n+1)=y(n)+1/6*(k1+2*k2+2*k3+k4);end4、用 RungeKutta 算法求解 P178例題、P285實驗任務(wù)（2），計算結(jié)果如下（其 中 yi表示數(shù)值解， y（xi） 表示解析解，結(jié)果保留八位有效數(shù)字） ：求解 P17

9、8 例題：0 0.1000 0.2000 0.3000 0.4000 0.50001.0000 1.1111 1.2500 1.4286 1.6667 2.0000xi0.050.10.150.20.25yi-0.0222-0.0388-0.0498-0.0552-0.0550y(xi)-0.0222-0.0388-0.0498-0.0552-0.0550xi0.30.350.40.450.5yi-0.0493-0.0380-0.02130.00100.0288y(xi)-0.0493-0.0380-0.02130.00100.02885、P285 實驗任務(wù)（ 2）精確解與近似解的圖形比較6、

10、用 matlab 中的 ode45求解 P285 實驗任務(wù)（ 2） ans =Columns 1 through 700.00010.00020.00030.00040.00090.00140-0.0001-0.0001-0.0002-0.0002-0.0005-0.0007Columns 8 through 140.00190.00240.00490.00740.00990.01250.0250-0.0010-0.0012-0.0024-0.0037-0.0049-0.0061-0.0118Columns 15 through 210.0375 0.0500 0.0625 0.0750 0.

11、0875 0.1000 0.1125-0.0172 -0.0222 -0.0268 -0.0312 -0.0351 -0.0388Columns 22 through 280.1250 0.1375 0.1500 0.1625 0.1750 0.1875-0.0450 -0.0475 -0.0497 -0.0516 -0.0532 -0.0543Columns 29 through 350.2125 0.2250 0.2375 0.2500 0.2625 0.2750-0.0556 -0.0558 -0.0556 -0.0550 -0.0541 -0.0528Columns 36 throug

12、h 420.3000 0.3125 0.3250 0.3375 0.3500 0.3625-0.0493 -0.0470 -0.0444 -0.0414 -0.0381 -0.0344Columns 43 through 490.3875 0.4000 0.4125 0.4250 0.4375 0.4500-0.0260 -0.0213 -0.0162 -0.0108 -0.0051 0.0010Columns 50 through 530.4718 0.4812 0.49060.50000.0125 0.0177 0.02320.0288-0.04200.2000-0.05520.2875-

13、0.05120.3750-0.03040.46250.0074五、討論分析 誤差一開始變大，然后變小，最后又變大六、改進(jìn)實驗建議可以通過提高保留的位數(shù)，使 Runge Kutta 算法結(jié)果更精確，也可以使數(shù)值解 和解析解的比較更加精確實驗四 常微分方程初值問題數(shù)值解法這里只提供解答格式請同學(xué)自己按照上機實驗報告格式填寫）1、餓 uler 法求解初值問題function x,y=euler_f(ydot_fun, x0, y0, h, N) % Euler (向前)公式，其中% ydot_fun -% x0, y0 -%h-%N-%x- Xn%y- Yn一階微分方程的函數(shù)初始條件區(qū)間步長區(qū)間的個

14、數(shù)構(gòu)成的向量構(gòu)成的向量x=zeros(1,N+1); y=zeros(1,N+1); x(1)=x0; y(1)=y0;for n=1:Nx(n+1)=x(n)+h;y(n+1)=y(n)+h*feval(ydot_fun, x(n), y(n);end用歐拉法計算 p167例 1ydot_fun=inline(y-2*x./y,x,y); x,y=euler_f(ydot_fun,0,1,0.1,10)x = Columns 1 through 110 0.100000000000000 0.200000000000000 0.300000000000000 0.50000000000000

15、00.600000000000000 0.700000000000000 0.800000000000000 1.000000000000000y =.000000000000000 1.1000000000000001.358212599560289 1.4351329186577961.508966253566332 1.5803382376552171.1918181818181821.6497834310477111.7847708324979820.4000000000000000.9000000000000001.2774378337147221.7177793478600872、

16、改進(jìn) Euler 公式function x,y=euler_r(ydot_fun, x0, y0, h, N)% 改進(jìn) Euler 公式，其中%ydot_fun - 一階微分方程的函數(shù)%x0, y0 - 初始條件%h-區(qū)間步長%N-區(qū)間的個數(shù)%x-Xn 構(gòu)成的向量%y-Yn 構(gòu)成的向量x=zeros(1,N+1); y=zeros(1,N+1); x(1)=x0; y(1)=y0; for n=1:Nx(n+1)=x(n)+h; ybar=y(n)+h*feval(ydot_fun, x(n), y(n); y(n+1)=y(n)+h/2*(feval(ydot_fun, x(n), y(n)

17、+feval(ydot_fun, x(n+1), ybar); end用改進(jìn)歐拉公式計算 P167 例題 1ydot_fun=inline(y-2*x./y,x,y); x,y=euler_r(ydot_fun,0,1,0.1,10)y =Columns 1 through 61.0000000000000001.343360151483999 Columns 7 through 11 1.485955602415669 1.7378674010354141.0959090909090911.4164019285369091.5525140913261461.1840965692429971.

18、6164747827520581.2662013608757761.678166363675186x =Columns 1 through 60 0.100000000000000 0.200000000000000 0.300000000000000 0.4000000000000000.500000000000000Columns 7 through 110.600000000000000 0.700000000000000 0.800000000000000 0.9000000000000001.000000000000000用歐拉法、改進(jìn)歐拉法計算 P171 例題 3ydot_fun=

19、inline(-y+x+1,x,y); x,y=euler_f(ydot_fun,0,1,0.1,5)x =00.1000000000000000.2000000000000000.3000000000000000.400000000000000 0.500000000000000y =1.000000000000000 1.000000000000000 1.010000000000000 1.056100000000000 1.090490000000000 用改進(jìn)歐拉法計算 P171例題 3 ydot_fun=inline(-y+x+1,x,y); x,y=euler_r(ydot_fu

20、n,0,1,0.1,5)x=0 0.100000000000000 0.200000000000000 0.400000000000000 0.5000000000000001.0290000000000000.300000000000000y=1.000000000000000 1.005000000000000 1.0190250000000001.070801950625000 1.1070757653156251.0412176250000003、標(biāo)準(zhǔn)四階 Runge_Kutta 公式function x,y=Runge_Kutta4(ydot_fun, x0, y0, h, N) %

21、 標(biāo)準(zhǔn)四階 Runge_Kutta 公式，其中 %ydot_fun - 一階微分方程的函數(shù)%x0, y0 - 初始條件% h - 區(qū)間步長%N- 區(qū)間的個數(shù)%x- Xn 構(gòu)成的向量% y - Yn 構(gòu)成的向量x=zeros(1,N+1); y=zeros(1,N+1); x(1)=x0; y(1)=y0; for n=1:Nx(n+1)=x(n)+h;k1=h*feval(ydot_fun, x(n), y(n);k2=h*feval(ydot_fun, x(n)+1/2*h, y(n)+1/2*k1);k3=h*feval(ydot_fun, x(n)+1/2*h, y(n)+1/2*k2)

22、; k4=h*feval(ydot_fun, x(n)+h, y(n)+k3); y(n+1)=y(n)+1/6*(k1+2*k2+2*k3+k4); end用四階龍格庫塔計算 P178 例題 ydot_fun=inline(y2,x,y); x,y=Runge_Kutta4(ydot_fun,0,1,0.1,5)x =0 0.100000000000000 0.200000000000000 0.400000000000000 0.5000000000000000.300000000000000y =1.000000000000000 1.111110490052194 1.24999799

23、2047015 1.4285661863014451.666653257250323 1.999963258950669用龍格庫塔方法計算 P285實驗任務(wù)( 2)ydot_fun=inline(-y+x2+4*x-1)./2,x,y); x,y=Runge_Kutta4(ydot_fun,0,0,0.05,10)x =Columns 1 through 6 0 0.050000000000000 0.250000000000000 Columns 7 through 11 0.300000000000000 0.500000000000000y = Columns 1 through 6 0

24、 -0.022190087076823 -0.055003093175772Columns 7 through 110.1000000000000000.1500000000000000.2000000000000000.3500000000000000.4000000000000000.450000000000000-0.038770573733692-0.049756511058554-0.055162578526671-0.049292018554665 -0.038042973450910-0.0212692404030080.0010162259975860.028800791009

25、418xx=0:0.05:0.5; z=exp(-xx./2)+xx.2-1; hold on plot(x,y,o);plot(xx,z,*)hold off 其圖形如圖一所示0.030.020.010-0.01-0.02-0.03-0.04-0.050 0.1 0.2 0.3 0.4 0.5-0.06圖一 精確解與龍格庫塔計算結(jié)果比較5、 用 Matlab 提供的 ode45求解 P285 實驗任務(wù)( 2)ydot_fun=inline(-y+x2+4*x-1)./2,x,y); x,y=ode45(ydot_fun, 0,0.5, 0); x;yplot(x,y, * )ans = Co

26、lumns 1 through 60 0.000100475457260 0.0002009509145210.0003014263717810.0004019018290420.0009042791153430 -0.000050226371419-0.000100430028501-0.000150610971371-0.000200769200158-0.0004512196372670.0014066564016450.0019090336879470.0024114109742490.0049232974057590.0074351838372680.009947070268778-

27、0.000701102241287-0.000950417028058-0.001199164013417-0.002434382472993-0.003655408270323-0.0048622433804000.0124589567002880.0249589567002880.0374589567002880.0499589567002880.0624589567002880.074958956700288-0.006054889775763-0.011778983079828-0.017151998201403-0.022174175468426-0.026845753781789-

28、0.0311669705745320.0874589567002880.0999589567002880.1124589567002880.1249589567002880.1374589567002880.149958956700288-0.035138061683373-0.038759261501758-0.042030803032876-0.044952917848809-0.047525835962155-0.0497497859783920.1624589567002880.1749589567002880.1874589567002880.1999589567002880.2124589567002880.224958956700288-0.05162