1、實 驗 報 告課程名稱 數值分析 實驗項目 常微分方程問題初值問題數值解法 一. 實驗目的1、 理解如何在計算機上實現用Euler法、改進Euler法、RungeKutta算法求一階常微分方程初值問題 的數值解。2、 利用圖形直觀分析近似解和準確解之間的誤差。3、 學會Matlab提供的ode45函數求解微分方程初值問題。二、實驗要求（1） 按照題目要求完成實驗內容；（2） 寫出相應的Matlab 程序；（3） 給出實驗結果(可以用表格展示實驗結果)；（4） 分析和討論實驗結果并提出可能的優化實驗。（5） 寫出實驗報告。三、實驗步驟1、用編好的Euler法、改進Euler法計算書本P167 的

2、例1、P171例題3。（1）取，求解初值問題（2）取，求解初值問題2、用RungeKutta算法計算P178例題、P285實驗任務（2）（1）取，求解初值問題（2）求初值問題的解在處的近似值，并與問題的解析解相比較。3、用Matlab繪圖函數plot(x,y)繪制P285實驗任務（2）的精確解和近似解的圖形。4、使用matlab中的ode45求解P285實驗任務（2），并繪圖。四、實驗結果1、Euler算法程序、改進Euler算法程序；Euler算法程序：function x,y=euler_f(ydot_fun, x0, y0, h, N)% Euler（向前）公式，其中% ydot_fun

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

4、件% h - 區間步長% N - 區間的個數% x - Xn 構成的向量% y - Yn 構成的向量x=zeros(1,N+1); y=zeros(1,N+1); x(1)=x0; y(1)=y0;for n=1:N x(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算法程序、改進Euler算法求解P167例題1的運行結果；（1.）Euler算法程序：x = Colu

5、mns 1 through 8 0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 Columns 9 through 11 0.8000 0.9000 1.0000y = Columns 1 through 8 1.0000 1.1000 1.1918 1.2774 1.3582 1.4351 1.5090 1.5803 Columns 9 through 11 1.6498 1.7178 1.7848（2）改進Euler算法：x = Columns 1 through 7 0 0.1000 0.2000 0.3000 0.4000 0.50

6、00 0.6000 Columns 8 through 11 0.7000 0.8000 0.9000 1.0000y = Columns 1 through 7 1.0000 1.0959 1.1841 1.2662 1.3434 1.4164 1.4860 Columns 8 through 11 1.5525 1.6165 1.6782 1.73793、RungeKutta算法程序；function x,y=Runge_Kutta4(ydot_fun, x0, y0, h, N)% 標準四階Runge_Kutta公式，其中% ydot_fun - 一階微分方程的函數% x0, y0 -

7、初始條件% h - 區間步長% N - 區間的個數% x - Xn 構成的向量% y - Yn 構成的向量x=zeros(1,N+1); y=zeros(1,N+1); x(1)=x0; y(1)=y0;for n=1:N x(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); k4=h*feval(ydot_fun, x(n)+h, y(n)+k3); y(n+1)=

8、y(n)+1/6*(k1+2*k2+2*k3+k4);end4、 用RungeKutta算法求解P178例題、P285實驗任務（2），計算結果如下（其中表示數值解，表示解析解，結果保留八位有效數字）：求解P178例題：x = 0 0.1000 0.2000 0.3000 0.4000 0.5000y = 1.0000 1.1111 1.2500 1.4286 1.6667 2.00000.050.10.150.20.25-0.0222 -0.0388-0.0498-0.0552-0.0550-0.0222 -0.0388-0.0498-0.0552-0.05500.30.350.40.450.

9、5-0.0493 -0.0380-0.02130.00100.0288-0.0493 -0.0380-0.02130.00100.02885、P285實驗任務（2）精確解與近似解的圖形比較6、用matlab中的ode45求解P285實驗任務（2）ans = Columns 1 through 7 0 0.0001 0.0002 0.0003 0.0004 0.0009 0.0014 0 -0.0001 -0.0001 -0.0002 -0.0002 -0.0005 -0.0007 Columns 8 through 14 0.0019 0.0024 0.0049 0.0074 0.0099 0

10、.0125 0.0250 -0.0010 -0.0012 -0.0024 -0.0037 -0.0049 -0.0061 -0.0118 Columns 15 through 21 0.0375 0.0500 0.0625 0.0750 0.0875 0.1000 0.1125 -0.0172 -0.0222 -0.0268 -0.0312 -0.0351 -0.0388 -0.0420 Columns 22 through 28 0.1250 0.1375 0.1500 0.1625 0.1750 0.1875 0.2000 -0.0450 -0.0475 -0.0497 -0.0516 -

11、0.0532 -0.0543 -0.0552 Columns 29 through 35 0.2125 0.2250 0.2375 0.2500 0.2625 0.2750 0.2875 -0.0556 -0.0558 -0.0556 -0.0550 -0.0541 -0.0528 -0.0512 Columns 36 through 42 0.3000 0.3125 0.3250 0.3375 0.3500 0.3625 0.3750 -0.0493 -0.0470 -0.0444 -0.0414 -0.0381 -0.0344 -0.0304 Columns 43 through 49 0

12、.3875 0.4000 0.4125 0.4250 0.4375 0.4500 0.4625 -0.0260 -0.0213 -0.0162 -0.0108 -0.0051 0.0010 0.0074 Columns 50 through 53 0.4718 0.4812 0.4906 0.50000.0125 0.0177 0.0232 0.0288五、討論分析誤差一開始變大，然后變小，最后又變大六、改進實驗建議可以通過提高保留的位數，使RungeKutta算法結果更精確，也可以使數值解和解析解的比較更加精確實驗四 常微分方程初值問題數值解法（這里只提供解答格式請同學自己按照上機實驗報告格

13、式填寫）1、餓uler法求解初值問題function x,y=euler_f(ydot_fun, x0, y0, h, N)% Euler（向前）公式，其中% ydot_fun - 一階微分方程的函數% x0, y0 - 初始條件% h - 區間步長% N - 區間的個數% x - Xn 構成的向量% y - Yn 構成的向量x=zeros(1,N+1); y=zeros(1,N+1); x(1)=x0; y(1)=y0;for n=1:N x(n+1)=x(n)+h; y(n+1)=y(n)+h*feval(ydot_fun, x(n), y(n);end用歐拉法計算p167例1>&g

14、t;ydot_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.400000000000000 0.5000000000000000.600000000000000 0.700000000000000 0.800000000000000 0.900000000000000 1.00000000

15、0000000y =.000000000000000 1.100000000000000 1.191818181818182 1.277437833714722 1.358212599560289 1.4351329186577961.508966253566332 1.580338237655217 1.649783431047711 1.717779347860087 1.7847708324979822、改進Euler公式function x,y=euler_r(ydot_fun, x0, y0, h, N)% 改進Euler公式，其中% ydot_fun - 一階微分方程的函數% x0

16、, y0 - 初始條件% h - 區間步長% N - 區間的個數% x - Xn 構成的向量% y - Yn 構成的向量x=zeros(1,N+1); y=zeros(1,N+1); x(1)=x0; y(1)=y0;for n=1:N x(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);end用改進歐拉公式計算P167例題1>>ydot_fun=inline('y-2*

17、x./y','x','y');>> x,y=euler_r(ydot_fun,0,1,0.1,10)x =Columns 1 through 60 0.100000000000000 0.200000000000000 0.300000000000000 0.400000000000000 0.500000000000000Columns 7 through 110.600000000000000 0.700000000000000 0.800000000000000 0.900000000000000 1.000000000000000y

18、=Columns 1 through 61.000000000000000 1.095909090909091 1.184096569242997 1.266201360875776 1.343360151483999 1.416401928536909Columns 7 through 111.485955602415669 1.552514091326146 1.616474782752058 1.678166363675186 1.737867401035414用歐拉法、改進歐拉法計算P171例題3ydot_fun=inline('-y+x+1','x',

19、'y');>> x,y=euler_f(ydot_fun,0,1,0.1,5)x =0 0.100000000000000 0.200000000000000 0.300000000000000 0.400000000000000 0.500000000000000y =1.000000000000000 1.000000000000000 1.010000000000000 1.029000000000000 1.056100000000000 1.090490000000000用改進歐拉法計算P171例題3ydot_fun=inline('-y+x+1&

20、#39;,'x','y');>> x,y=euler_r(ydot_fun,0,1,0.1,5)x=0 0.100000000000000 0.200000000000000 0.300000000000000 0.400000000000000 0.500000000000000y=1.000000000000000 1.005000000000000 1.019025000000000 1.041217625000000 1.070801950625000 1.1070757653156253、標準四階Runge_Kutta公式function

21、x,y=Runge_Kutta4(ydot_fun, x0, y0, h, N)% 標準四階Runge_Kutta公式，其中% ydot_fun - 一階微分方程的函數% x0, y0 - 初始條件% h - 區間步長% N - 區間的個數% x - Xn 構成的向量% y - Yn 構成的向量x=zeros(1,N+1); y=zeros(1,N+1); x(1)=x0; y(1)=y0;for n=1:N x(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)

22、; 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);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.300000000

23、000000 0.400000000000000 0.500000000000000y =1.000000000000000 1.111110490052194 1.249997992047015 1.428566186301445 1.666653257250323 1.999963258950669用龍格庫塔方法計算P285實驗任務（2）>>ydot_fun=inline('(-y+x2+4*x-1)./2','x','y');>> x,y=Runge_Kutta4(ydot_fun,0,0,0.05,10)x =Co

24、lumns 1 through 60 0.050000000000000 0.100000000000000 0.150000000000000 0.200000000000000 0.250000000000000Columns 7 through 110.300000000000000 0.350000000000000 0.400000000000000 0.450000000000000 0.500000000000000y = Columns 1 through 60 -0.022190087076823 -0.038770573733692 -0.049756511058554 -

25、0.055162578526671 -0.055003093175772Columns 7 through 11-0.049292018554665 -0.038042973450910 -0.021269240403008 0.001016225997586 0.028800791009418>>xx=0:0.05:0.5;>> z=exp(-xx./2)+xx.2-1;>> hold on>> plot(x,y,'o');plot(xx,z,'*')>>hold off其圖形如圖一所示圖一 精確解與

26、龍格庫塔計算結果比較5、 用Matlab提供的ode45求解P285實驗任務（2）>>ydot_fun=inline('(-y+x2+4*x-1)./2','x','y');>> x,y=ode45(ydot_fun, 0,0.5, 0);>> x'y'>>plot(x,y,*)ans = Columns 1 through 60 0.000100475457260 0.000200950914521 0.000301426371781 0.000401901829042 0.000

27、904279115343 0 -0.000050226371419 -0.000100430028501 -0.000150610971371 -0.000200769200158 -0.000451219637267 0.001406656401645 0.001909033687947 0.002411410974249 0.004923297405759 0.007435183837268 0.009947070268778 -0.000701102241287 -0.000950417028058 -0.001199164013417 -0.002434382472993 -0.003

28、655408270323 -0.0048622433804000.012458956700288 0.024958956700288 0.037458956700288 0.049958956700288 0.062458956700288 0.074958956700288 -0.006054889775763 -0.011778983079828 -0.017151998201403 -0.022174175468426 -0.026845753781789 -0.0311669705745320.087458956700288 0.099958956700288 0.112458956700288 0.124958956700288 0.137458956700288 0.149958956700288-0.035138061683373 -0.038759261501758 -0.042030803032876 -0.044952917848809 -0.047525835962155 -0.0497497859783920.162458956700288 0.174958956700288 0.187458956700288 0.199958956700288 0.212458956700288 0.224958956700288-0.051