實驗四常微分方程初值問題數(shù)值解法重點講義

1、佛山科學(xué)技術(shù)學(xué)院實 驗 報 告課程名稱 數(shù)值分析 實驗項目 常微分方程初值問題數(shù)值解法 專業(yè)班級 12.數(shù)學(xué)與應(yīng)用數(shù)學(xué)（師范） 姓名 葉楚欣 學(xué)號 2012214103 指導(dǎo)教師 黃國順 成 績 日 期 一. 實驗?zāi)康?、 理解如何在計算機上實現(xiàn)用Euler法、改進Euler法、RungeKutta算法求一階常微分方程初值問題 的數(shù)值解。2、 利用圖形直觀分析近似解和準確解之間的誤差。3、 學(xué)會Matlab提供的ode45函數(shù)求解微分方程初值問題。二、實驗要求（1） 按照題目要求完成實驗內(nèi)容；（2） 寫出相應(yīng)的Matlab 程序；（3） 給出實驗結(jié)果(可以用表格展示實驗結(jié)果)；（4） 分析和討

2、論實驗結(jié)果并提出可能的優(yōu)化實驗。（5） 寫出實驗報告。三、實驗步驟1、用編好的Euler法、改進Euler法計算書本P167 的例1、P171例題3。（1）取，求解初值問題（2）取，求解初值問題2、用RungeKutta算法計算P178例題、P285實驗任務(wù)（2）（1）取，求解初值問題（2）求初值問題的解在處的近似值，并與問題的解析解相比較。3、用Matlab繪圖函數(shù)plot(x,y)繪制P285實驗任務(wù)（2）的精確解和近似解的圖形。4、使用matlab中的ode45求解P285實驗任務(wù)（2），并繪圖。四、實驗結(jié)果1、Euler算法程序、改進Euler算法程序；2、用Euler算法程序、改進E

3、uler算法求解P167例題1的運行結(jié)果；3、RungeKutta算法程序；4、用RungeKutta算法求解P178例題、P285實驗任務(wù)（2），計算結(jié)果如下（其中表示數(shù)值解，表示解析解，結(jié)果保留八位有效數(shù)字）：0.050.10.150.20.25-0.02219009-0.03877057-0.04975651-0.05516258-0.05500309-0.02219009-0.03877058-0.04975651-0.05516258-0.055003100.30.350.40.450.5-0.04929202 -0.03804297-0.021269240.001016230.02

4、880079-0.04929202-0.03804298-0.021269250.001016220.028800785、P285實驗任務(wù)（2）精確解與近似解的圖形比較（圖片貼到此處）6、用matlab中的ode45求解P285實驗任務(wù)（2） MATLAB desktop keyboard shortcuts, such as Ctrl+S, are now customizable. In addition, many keyboard shortcuts have changed for improved consistency across the desktop. To customi

5、ze keyboard shortcuts, use Preferences. From there, you can also restore previous default settings by selecting "R2009a Windows Default Set" from the active settings drop-down list. For more information, see Help. Click here if you do not want to see this message again. >> ydot_fun=i

6、nline('y-2*x./y','x','y');>> x,y=euler_f(ydot_fun,0,1,0.1,10)x = Columns 1 through 4 0 0.100000000000000 0.200000000000000 0.300000000000000 Columns 5 through 8 0.400000000000000 0.500000000000000 0.600000000000000 0.700000000000000 Columns 9 through 11 0.80000000000000

7、0 0.900000000000000 1.000000000000000y = Columns 1 through 4 1.000000000000000 1.100000000000000 1.191818181818182 1.277437833714722 Columns 5 through 8 1.358212599560289 1.435132918657796 1.508966253566332 1.580338237655217 Columns 9 through 11 1.649783431047711 1.717779347860087 1.784770832497982&

8、gt;> ydot_fun=inline('y-2*x./y','x','y');>> x,y=euler_r(ydot_fun,0,1,0.1,10)x = Columns 1 through 4 0 0.100000000000000 0.200000000000000 0.300000000000000 Columns 5 through 8 0.400000000000000 0.500000000000000 0.600000000000000 0.700000000000000 Columns 9 through 1

9、1 0.800000000000000 0.900000000000000 1.000000000000000y = Columns 1 through 4 1.000000000000000 1.095909090909091 1.184096569242997 1.266201360875776 Columns 5 through 8 1.343360151483999 1.416401928536909 1.485955602415669 1.552514091326146 Columns 9 through 11 1.616474782752058 1.678166363675186

10、1.737867401035414>> ydot_fun=inline('-y+x+1','x','y');>> x,y=euler_f(ydot_fun,0,1,0.1,5)x = Columns 1 through 4 0 0.100000000000000 0.200000000000000 0.300000000000000 Columns 5 through 6 0.400000000000000 0.500000000000000y = Columns 1 through 4 1.000000000000000

11、 1.000000000000000 1.010000000000000 1.029000000000000 Columns 5 through 6 1.056100000000000 1.090490000000000>> ydot_fun=inline('-y+x+1','x','y');>> x,y=euler_r(ydot_fun,0,1,0.1,5)x = Columns 1 through 4 0 0.100000000000000 0.200000000000000 0.300000000000000 Col

12、umns 5 through 6 0.400000000000000 0.500000000000000y = Columns 1 through 4 1.000000000000000 1.005000000000000 1.019025000000000 1.041217625000000 Columns 5 through 6 1.070801950625000 1.107075765315625>> ydot_fun=inline('y2','x','y');>> x,y=Runge_Kutta4(ydot_fun

13、,0,1,0.1,5)x = Columns 1 through 4 0 0.100000000000000 0.200000000000000 0.300000000000000 Columns 5 through 6 0.400000000000000 0.500000000000000y = Columns 1 through 4 1.000000000000000 1.111110490052194 1.249997992047015 1.428566186301445 Columns 5 through 6 1.666653257250323 1.999963258950669>

14、;> 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 4 0 0.050000000000000 0.100000000000000 0.150000000000000 Columns 5 through 8 0.200000000000000 0.250000000000000 0.300000000000000 0.350000000000000 Columns

15、 9 through 11 0.400000000000000 0.450000000000000 0.500000000000000y = Columns 1 through 4 0 -0.022190087076823 -0.038770573733692 -0.049756511058554 Columns 5 through 8 -0.055162578526671 -0.055003093175772 -0.049292018554665 -0.038042973450910 Columns 9 through 11 -0.021269240403008 0.001016225997

16、586 0.028800791009418ydot_fun=inline('(-y+x2+4*x-1)./2','x','y');x,y=Runge_Kutta4(ydot_fun,0,0,0.05,10); sprintf('%.8f,%.8fn',x,y)ans =0.00000000,0.050000000.10000000,0.150000000.20000000,0.250000000.30000000,0.350000000.40000000,0.450000000.50000000,0.00000000-0.0221

17、9009,-0.03877057-0.04975651,-0.05516258-0.05500309,-0.04929202-0.03804297,-0.021269240.00101623,0.02880079xx=0:0.05:0.5;z=exp(-xx./2)+xx.2-1;sprintf('%.8f,%.8fn',xx,z)ans =0.00000000,0.050000000.10000000,0.150000000.20000000,0.250000000.30000000,0.350000000.40000000,0.450000000.50000000,0.00

18、000000-0.02219009,-0.03877058-0.04975651,-0.05516258-0.05500310,-0.04929202-0.03804298,-0.021269250.00101622,0.02880078>> xx=0:0.05:0.5;z=exp(-xx./2)+xx.2-1;hold on>> plot(x,y,'o');plot(xx,z,'*');>> hold off>> ydot_fun=inline('(-y+x2+4*x-1)./2','x&

19、#39;,'y');>> x,y=ode45(ydot_fun,0,0.5,0);>> x'y'ans = Columns 1 through 4 0 0.000100475457260 0.000200950914521 0.000301426371781 0 -0.000050226371419 -0.000100430028501 -0.000150610971371 Columns 5 through 8 0.000401901829042 0.000904279115343 0.001406656401645 0.0019090

20、33687947 -0.000200769200158 -0.000451219637267 -0.000701102241287 -0.000950417028058 Columns 9 through 12 0.002411410974249 0.004923297405759 0.007435183837268 0.009947070268778 -0.001199164013417 -0.002434382472993 -0.003655408270323 -0.004862243380400 Columns 13 through 16 0.012458956700288 0.0249

21、58956700288 0.037458956700288 0.049958956700288 -0.006054889775763 -0.011778983079828 -0.017151998201403 -0.022174175468426 Columns 17 through 20 0.062458956700288 0.074958956700288 0.087458956700288 0.099958956700288 -0.026845753781789 -0.031166970574532 -0.035138061683373 -0.038759261501758 Column

22、s 21 through 24 0.112458956700288 0.124958956700288 0.137458956700288 0.149958956700288 -0.042030803032876 -0.044952917848809 -0.047525835962155 -0.049749785978392 Columns 25 through 28 0.162458956700288 0.174958956700288 0.187458956700288 0.199958956700288 -0.051624995148617 -0.053151689328665 -0.0

23、54330092850809 -0.055160428675467 Columns 29 through 32 0.212458956700288 0.224958956700288 0.237458956700288 0.249958956700288 -0.055642918443663 -0.055777782436119 -0.055565239445031 -0.055005506925134 Columns 33 through 36 0.262458956700288 0.274958956700288 0.287458956700288 0.299958956700287 -0

24、.054098801045891 -0.052845336650564 -0.051245327128055 -0.049298984563352 Columns 37 through 40 0.312458956700288 0.324958956700288 0.337458956700287 0.349958956700287 -0.047006519789452 -0.044368142346404 -0.041384060353229 -0.038054480657744 Columns 41 through 44 0.362458956700287 0.37495895670028

25、8 0.387458956700288 0.399958956700287 -0.034379608888238 -0.030359649412486 -0.025994805209740 -0.021285278019953 Columns 45 through 48 0.412458956700287 0.424958956700287 0.437458956700287 0.449958956700287 -0.016231268395211 -0.010832975658719 -0.005090597776859 0.000995668492165 Columns 49 throug

26、h 52 0.462458956700287 0.471844217525216 0.481229478350144 0.490614739175072 0.007425627547137 0.012479158933356 0.017726249535970 0.023166817935270 Column 53 0.500000000000000 0.028800783076419>> plot(x,y,'*')>>五、討論分析（通過最后繪制的圖形直觀分析近似解與準確解之間的誤差，自己補充）從圖中可以看到4階的Runge-Kutta畫出來的點基本和原

27、函數(shù)的對應(yīng)點相重合，精度之高，從數(shù)據(jù)來看，起碼有9位有效數(shù)字。六、改進實驗建議（自己補充）可以通過輸出數(shù)據(jù)的精度控制函數(shù)實驗四 常微分方程初值問題數(shù)值解法（這里只提供解答格式請同學(xué)自己按照上機實驗報告格式填寫）1、餓uler法求解初值問題function x,y=euler_f(ydot_fun, x0, y0, h, 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

28、)=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>>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.40000000000000

29、0 0.5000000000000000.600000000000000 0.700000000000000 0.800000000000000 0.900000000000000 1.000000000000000y =.000000000000000 1.100000000000000 1.191818181818182 1.277437833714722 1.358212599560289 1.4351329186577961.508966253566332 1.580338237655217 1.649783431047711 1.717779347860087 1.784770832

30、4979822、改進Euler公式function x,y=euler_r(ydot_fun, x0, y0, h, 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: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(yd

31、ot_fun, x(n), y(n)+feval(ydot_fun, x(n+1), ybar);end用改進歐拉公式計算P167例題1>>ydot_fun=inline('y-2*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

32、through 110.600000000000000 0.700000000000000 0.800000000000000 0.900000000000000 1.000000000000000y =Columns 1 through 61.000000000000000 1.095909090909091 1.184096569242997 1.266201360875776 1.343360151483999 1.416401928536909Columns 7 through 111.485955602415669 1.552514091326146 1.61647478275205

33、8 1.678166363675186 1.737867401035414用歐拉法、改進歐拉法計算P171例題3ydot_fun=inline('-y+x+1','x','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.01000000000

34、0000 1.029000000000000 1.056100000000000 1.090490000000000用改進歐拉法計算P171例題3ydot_fun=inline('-y+x+1','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.00500000000000

35、0 1.019025000000000 1.041217625000000 1.070801950625000 1.1070757653156253、標準四階Runge_Kutta公式function x,y=Runge_Kutta4(ydot_fun, x0, y0, h, 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;fo

36、r 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)=y(n)+1/6*(k1+2*k2+2*k3+k4);end用四階龍格庫塔計算P178例題>> ydot_fun=inline('y2','x','y

37、9;);>> x,y=Runge_Kutta4(ydot_fun,0,1,0.1,5)x =0 0.100000000000000 0.200000000000000 0.300000000000000 0.400000000000000 0.500000000000000y =1.000000000000000 1.111110490052194 1.249997992047015 1.428566186301445 1.666653257250323 1.999963258950669用龍格庫塔方法計算P285實驗任務(wù)（2）>>ydot_fun=inline(

38、9;(-y+x2+4*x-1)./2','x','y');>> x,y=Runge_Kutta4(ydot_fun,0,0,0.05,10)x =Columns 1 through 60 0.050000000000000 0.100000000000000 0.150000000000000 0.200000000000000 0.250000000000000Columns 7 through 110.300000000000000 0.350000000000000 0.400000000000000 0.450000000000000

39、 0.500000000000000y = Columns 1 through 60 -0.022190087076823 -0.038770573733692 -0.049756511058554 -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

40、;>> hold on>> plot(x,y,'o');plot(xx,z,'*')>>hold off其圖形如圖一所示圖一 精確解與龍格庫塔計算結(jié)果比較4、 用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'y'>>plot(x,y,*)ans

41、 = Columns 1 through 60 0.000100475457260 0.000200950914521 0.000301426371781 0.000401901829042 0.000904279115343 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.003655408270323 -0.0048622433804000.012458956700288 0.024958956700288 0.037458956700288 0.049958956700288 0.062458956700288 0.074958956700288 -0.0060548897757