實(shí)驗(yàn)四 常微分方程初值問題數(shù)值解法

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

2、驟1、用編好的Euler法、改進(jìn)Euler法計(jì)算書本P167 的例1、P171例題3。（1）取，求解初值問題（2）取，求解初值問題2、用RungeKutta算法計(jì)算P178例題、P285實(shí)驗(yàn)任務(wù)（2）（1）取，求解初值問題（2）求初值問題的解在處的近似值，并與問題的解析解相比較。3、用Matlab繪圖函數(shù)plot(x,y繪制P285實(shí)驗(yàn)任務(wù)（2）的精確解和近似解的圖形。4、使用matlab中的ode45求解P285實(shí)驗(yàn)任務(wù)（2），并繪圖。四、實(shí)驗(yàn)結(jié)果1、Euler算法程序、改進(jìn)Euler算法程序；euler法求解初值問題function x,y=euler_f(ydot_fun, x0, y0

3、, h, N% Euler（向前）公式，其中% ydot_fun - 一階微分方程的函數(shù)% x0, y0 - 初始條件% h - 區(qū)間步長% N - 區(qū)間的個(gè)數(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;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 - 一階微分方程的

4、函數(shù)% x0, y0 - 初始條件% h - 區(qū)間步長% N - 區(qū)間的個(gè)數(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+feval(ydot_fun, x(n+1, ybar;end2、用Euler算法程序、改進(jìn)Euler算法求解P167例題1和 P171例題3的運(yùn)行結(jié)果；用歐拉法計(jì)算：p167例1：&

5、gt;>ydot_fun=inline('y-2*x./y','x','y'>> x,y=euler_f(ydot_fun,0,1,0.1,10x = 0 0.100000000000000 0.200000000000000 0.300000000000000 0.400000000000000 0.500000000000000 0.600000000000000 0.700000000000000 0.800000000000000 0.900000000000000 1.000000000000000y =0.00000

6、0000000000 1.100000000000000 1.191818181818182 1.277437833714722 1.358212599560289 1.435132918657796 1.508966253566332 1.580338237655217 1.649783431047711 1.717779347860087 1.784770832497982用改進(jìn)歐拉公式計(jì)算：>>ydot_fun=inline('y-2*x./y','x','y'>> x,y=euler_r(ydot_fun,0,1,

7、0.1,10x =0 0.100000000000000 0.200000000000000 0.300000000000000 0.400000000000000 0.500000000000000 0.600000000000000 0.700000000000000 0.800000000000000 0.900000000000000 1.000000000000000y = 1.000000000000000 1.095909090909091 1.184096569242997 1.266201360875776 1.343360151483999 1.41640192853690

8、9 1.485955602415669 1.552514091326146 1.616474782752058 1.678166363675186 1.737867401035414P171例題3h=0.1；用歐拉法計(jì)算：ydot_fun=inline('-y+x+1','x','y'>> x,y=euler_f(ydot_fun,0,1,0.1,5x =0 0.100000000000000 0.200000000000000 0.300000000000000 0.400000000000000 0.500000000000000

9、y =1.000000000000000 1.000000000000000 1.010000000000000 1.029000000000000 1.056100000000000 1.090490000000000用改進(jìn)歐拉法計(jì)算：ydot_fun=inline('-y+x+1','x','y'>> x,y=euler_r(ydot_fun,0,1,0.1,5x=0 0.100000000000000 0.200000000000000 0.300000000000000 0.400000000000000 0.50000000

10、0000000y=1.000000000000000 1.005000000000000 1.019025000000000 1.041217625000000 1.070801950625000 1.1070757653156253、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ū)間的個(gè)數(shù)% x - Xn 構(gòu)成的向量% y - Yn 構(gòu)成的向量x=zeros(1,N

11、+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;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實(shí)驗(yàn)任務(wù)（2），計(jì)算結(jié)果如下（其中表示數(shù)值解，表示解析解，結(jié)果

12、保留八位有效數(shù)字）：解：P178>> ydot_fun=inline('y2','x','y'>> x,y=Runge_Kutta4(ydot_fun,0,1,0.05,10;>>sprintf('%.8f %.8f n',x,y 由以上程序可以得出的近似解；由常微分方法算出其的y(x=1/(1-x;ydot_fun=inline('1/(1-x','x','y'>> x,y=Runge_Kutta4(ydot_fun,0,1,0.05

13、,10;>>sprintf('%.8f %.8f n',x,y 由以上程序可以得出的解析解，列表如下：0.050.10.150.20.251.05263156 1.111111071.176470511.249999871.333333121.052631581.11111111 1.176470591.250000001.333333330.30.350.40.450.51.428571091.538461001.666665801.818180401.99999761、1.42857143 1.53846154 1.66666667 1.818181822.00

14、000000P285>>ydot_fun=inline('(-y+x2+4*x-1./2','x','y'>> x,y=Runge_Kutta4(ydot_fun,0,0,0.05,10; sprintf('%.8f %.8f n',x,y 由以上程序可以得出的近似值；>>ydot_fun=inline('exp（-x./2+x.2-1','x','y'>> x,y=Runge_Kutta4(ydot_fun,0,0,0.05,10;

15、sprintf('%.8f %.8f n',x,y 由以上程序可以得出解析解,列表如下：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.02880079-0.04929202-0.03804298-0.02126925 0.001016220.0

16、28800785、P285實(shí)驗(yàn)任務(wù)（2）精確解與近似解的圖形比較>>xx=0:0.05:0.5;>> z=exp(-xx./2+xx.2-1;>> hold on>> plot(x,y,'o'plot(xx,z,'*'>>hold off其圖形如圖一所示： 圖一 精確解與龍格庫塔計(jì)算結(jié)果比較>>ydot_fun=inline('(-y+x2+4*x-1./2','x','y'>> x,y=ode45(ydot_fun, 0,0.5,

17、 0;>> x'y'>> plot(x,y,'*'其圖形如圖二所示：圖二 用ode45求解得到的近似解散點(diǎn)圖6、用matlab中的ode45求解P285實(shí)驗(yàn)任務(wù)（2）解：>>ydot_fun=inline('(-y+x2+4*x-1./2','x','y'>> x,y=ode45(ydot_fun, 0,0.5, 0;>> x'y'得：ans = Columns 1 through 60 0.000100475457260 0.000200

18、950914521 0.000301426371781 0.000401901829042 0.000904279115343 0 -0.000050226371419 -0.000100430028501 -0.000150610971371 -0.000200769200158 -0.0004512196372670.001406656401645 0.001909033687947 0.002411410974249 0.004923297405759 0.007435183837268 0.009947070268778-0.000701102241287 -0.00095041702

19、8058 -0.001199164013417 -0.002434382472993 -0.003655408270323 -0.0048622433804000.012458956700288 0.024958956700288 0.037458956700288 0.049958956700288 0.062458956700288 0.074958956700288-0.006054889775763 -0.011778983079828 -0.017151998201403 -0.022174175468426 -0.026845753781789 -0.031166970574532

20、0.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.212458

21、956700288 0.224958956700288-0.051624995148617 -0.053151689328665 -0.054330092850809 -0.055160428675467 -0.055642918443663 -0.0557777824361190.237458956700288 0.249958956700288 0.262458956700288 0.274958956700288 0.287458956700288 0.299958956700287-0.055565239445031 -0.055005506925134 -0.054098801045

22、891 -0.052845336650564 -0.051245327128055 -0.0492989845633520.312458956700288 0.324958956700288 0.337458956700287 0.349958956700287 0.362458956700287 0.374958956700288-0.047006519789452 -0.044368142346404 -0.041384060353229 -0.038054480657744 -0.034379608888238 -0.0303596494124860.387458956700288 0.399958956700287 0.412458956700287 0.424958956700287 0.43745895670