常微分方程初值問題數(shù)值解法比較

1、 常微分方程初值問題數(shù)值解法的比較 數(shù)值計(jì)算實(shí)踐課程設(shè)計(jì)報(bào)告課題名稱常微分方程初值問題數(shù)值解法的比較完成時(shí)間2013-1-17姓名班級(jí)學(xué)號(hào)成績一. 實(shí)驗(yàn)?zāi)康募皟?nèi)容1實(shí)驗(yàn)?zāi)康模海?） 了解常微分方程初值問題的理論背景以及初值問題穩(wěn)定性、收斂性的研究；（2） 熟練掌握歐拉法、改進(jìn)歐拉法、龍格-庫塔法以及截?cái)嗾`差分析；（3） 比較歐拉法、改進(jìn)歐拉法及龍格-庫塔法，能夠選擇合適的方法進(jìn)行問題的研究計(jì)算； 2實(shí)驗(yàn)內(nèi)容：求微分方程（歐拉法求解）求微分方程（改進(jìn)歐拉法求解）求微分方程（龍格-庫塔求解）根據(jù)實(shí)驗(yàn)的結(jié)果進(jìn)行分析，了解一般方法的的優(yōu)缺點(diǎn)，穩(wěn)定性，收斂性以及截?cái)嗾`差的分析，針對相應(yīng)問題拿出有效方法得

2、出最優(yōu)的結(jié)果。二相關(guān)背景知識(shí)介紹以及初值問題穩(wěn)定性的研究： 在科學(xué)與工程問題中，常微分方程表述物理量的變化規(guī)律，應(yīng)用非常廣泛，比如，天體運(yùn)動(dòng)的軌跡，機(jī)器人控制，化學(xué)反應(yīng)過程的描述和控制以及電路瞬態(tài)過程分析等。這些問題中要求解隨時(shí)間變化的物理量，即位置函數(shù)表示時(shí)間，而微分方程描述了未知函數(shù)與它的一階或高階導(dǎo)數(shù)之間的關(guān)系。 考慮一階常微分方程的初值問題 ，如果存在實(shí)數(shù)使得則稱f關(guān)于y滿足利普希茨條件，L稱為利普希茨常數(shù)。 對于常微分方程初值問題，考慮初值的擾動(dòng)是問題的解發(fā)生偏差的情形。若時(shí)的偏差被控制在有界范圍內(nèi)，則稱該初值問題是穩(wěn)定的，否則該初值問題不穩(wěn)定的。特別地，若時(shí)的偏差收斂于零，則稱該初

3、值問題是漸進(jìn)穩(wěn)定的。 對于初值問題穩(wěn)定性的研究，易知其準(zhǔn)確解為，假定初值經(jīng)過擾動(dòng)后變?yōu)?，對于擾動(dòng)后的解為因此帶來的擾動(dòng)誤差為，因此考慮時(shí)的值，它取決于。易知，若 ，則原問題是穩(wěn)定的；若 ，則原問題是不穩(wěn)定的；若 ，則原問題是漸進(jìn)穩(wěn)定的。實(shí)際遇到的大多數(shù)常微分方程初值問題都是穩(wěn)定的，因此在后面的討論數(shù)值解法時(shí)這常常是默認(rèn)條件。 1 歐拉法：依據(jù)：積分曲線上一點(diǎn)的切線斜率等于函數(shù)值。方法：推進(jìn)法，初始點(diǎn)出發(fā)，依照方向場在改點(diǎn)的方向推進(jìn)到 向前歐拉法的得到:(1)將在處泰勒展開取h的線性部分，得（2）將初值問題中得導(dǎo)數(shù)用向前差商來代替有 ，因此（3）將兩邊同時(shí)對x的區(qū)間上積分 對右端用左矩形公式得，

4、此方法稱向前歐拉法，也叫顯示歐拉法。（4）對右端用右矩形公式得，也叫隱式歐拉法。誤差分析：1.稱為計(jì)算時(shí)的局部截?cái)嗾`差； 2.如果數(shù)值方法的局部截?cái)嗾`差為，那么稱這種數(shù)值方法的階數(shù)是p，其實(shí)p為非負(fù)整數(shù)。通常情況下，步長h越小，p越高，則局部截?cái)嗾`差越小，計(jì)算精初泰勒展開有則有可見歐拉方法是一階方法，精度不是很高。2 改進(jìn)歐拉方法： 梯形公式：對右端用梯形公式得+顯然梯形公式是隱式公式。 改進(jìn)歐拉公式：先用歐拉公式求的一個(gè)初步的近似值，成為預(yù)測值，預(yù)測值的精度可能達(dá)不到要求，在用梯形公式將他校正一次，記為，這個(gè)結(jié)果成為校正值。 預(yù)測： 校正： 誤差分析：記為改進(jìn)歐拉公式在處的截?cái)嗾`差， 記因此

5、有，表示在出的局部截?cái)嗾`差。由此得，梯形公式的局部截?cái)嗾`差為，因此改進(jìn)歐拉的截?cái)嗾`差為，可見改進(jìn)歐拉的方法是二階方法，改進(jìn)歐拉方法優(yōu)于歐拉方法。3 龍格庫塔法: 根據(jù)拉格朗日微分中值定理，記得到，這樣，只要給出一種計(jì)算的算法，就可以得到相應(yīng)的計(jì)算公式。歐拉公式可以寫為 改進(jìn)歐拉公式可以寫成因此推出一般的推出廣式稱為p階龍格-庫塔方法，簡稱p階R-K方法。因?yàn)?這里的均為常數(shù)。因?yàn)榻o定的系數(shù)不唯一，因此這里的常數(shù)有無窮多個(gè)解，下面是特殊情況下和一般情況下得結(jié)果（一階龍格-庫塔）當(dāng)r=1時(shí)，這就是歐拉法。（二階龍格-庫塔）當(dāng)r=2時(shí)，這就是改進(jìn)歐拉法。三階和四階龍格-庫塔也只是在一般情況下得結(jié)果。

6、三階四階其局部截?cái)嗾`差是： 三程序代碼 歐拉法代碼如下：（1）向前歐拉法：function y=Euler1(fun,x0,y0,xN,N)%fun為一階微分方程，x0，y0為初始條件，xN為取值范圍的一個(gè)端點(diǎn)，N為區(qū)間個(gè)數(shù)x=zeros(1,N+1);x=zeros(1,N+1);x(1)=x0;y(1)=y0; h=(xN-x0)/N; %h為區(qū)間步長for(n=1:N) x(n+1)=x(n)+h; y(n+1)=y(n)+h*feval(fun,x(n),y(n); %根據(jù)向前歐拉公式計(jì)算y值endT=x,y（2）向后歐拉法：function y=Euler2(fun,x0,y0,xN

7、,N)% fun為一階微分方程，x0，y0為初始條件，xN為取值范圍的一個(gè)端點(diǎn)，N為區(qū)間個(gè)數(shù)x=zeros(1,N+1);x=zeros(1,N+1);x(1)=x0;y(1)=y0;h=(xN-x0)/N;% h為區(qū)間步長for(n=1:N) x(n+1)=x(n)+h; z0=y(n)+h*feval(fun,x(n),y(n); for(k=1:3) z1=y(n)+h*feval(fun,x(n+1),z0); if(abs(z1-z0)1e-3) break; end z0=z1; end y(n+1)=z1; %根據(jù)向后歐拉公式計(jì)算y值endT=x,y改進(jìn)歐拉代碼如下：functi

8、on x,y=Gaijineuler(f,x0,y0,xZ,h)%f為一階微分方程的函數(shù)；x0，y0為初始條件；xZ為取值范圍的一個(gè)端點(diǎn)，h為區(qū)間步長n=fix(xZ-x0)/h); %計(jì)算分點(diǎn)數(shù)y(1)=y0;x(1)=x0;for i=1:n x(i+1)=x0+i*h; yp=y(i)+h*feval(f,x(i),y(i); yc=y(i)+h*feval(f,x(i+1),yp); y(i+1)=(yp+yc)/2; %根據(jù)改進(jìn)歐拉公式計(jì)算結(jié)果endx=x;y=y;13.龍格-庫塔代碼如下：（三階龍格-庫塔）function R=Longgekuta3(f,a,b,aZ,h)%a，b

9、為端點(diǎn)，h為步長，aZ為初值n=(b-a)/h; T=zeros(1,n+1); %定義向量 Y=zeros(1,n+1); T=a:h:b; %計(jì)算各個(gè)分點(diǎn)Y(1)=aZ; %初值賦予for j=1:n k1=feval(f,T(j),Y(j); k2=feval(f,T(j)+h/2,Y(j)+k1*h/2); k3=feval(f,T(j)+h,Y(j)-h*k1+k2*2*h); Y(j+1)=Y(j)+(k1+4*k2+k3)*h/6; %根據(jù)公式計(jì)算endR=T Y;（四階龍格-庫塔）function R=Longgekuta4(f,a,b,aZ,h)% a，b為端點(diǎn)，h為步長，a

10、Z為初值n=(b-a)/h; T=zeros(1,n+1); %定義向量Y=zeros(1,n+1); T=a:h:b; %計(jì)算各個(gè)分點(diǎn)Y(1)=aZ; %初值賦予for j=1:n k1=feval(f,T(j),Y(j); k2=feval(f,T(j)+h/2,Y(j)+k1*h/2); k3=feval(f,T(j)+h/2,Y(j)+k2*h/2); k4=feval(f,T(j)+h,Y(j)+k3*h); Y(j+1)=Y(j)+(k1+2*k2+2*k3+k4)*h/6; %根據(jù)公式計(jì)算endR=T Y;四數(shù)值結(jié)果：輸入：定義M文件并保存ffx.m Euler1(ffx,0,1

11、,1,10)結(jié)果：T = 0 1.0000 0.1000 1.1000 0.2000 1.2100 0.3000 1.3310 0.4000 1.4641 0.5000 1.6105 0.6000 1.7716 0.7000 1.9487 0.8000 2.1436 0.9000 2.3579 1.0000 2.5937ans = Columns 1 through 7 1.0000 1.1000 1.2100 1.3310 1.4641 1.6105 1.7716 Columns 8 through 11 1.9487 2.1436 2.3579 2.5937向前歐拉公式結(jié)果：輸入：Eule

12、r2(ffx,0,1,1,10)結(jié)果T = 0 1.0000 0.1000 1.1110 0.2000 1.2344 0.3000 1.3716 0.4000 1.5240 0.5000 1.6933 0.6000 1.8814 0.7000 2.0904 0.8000 2.3227 0.9000 2.5807 1.0000 2.8674ans = Columns 1 through 7 1.0000 1.1110 1.2344 1.3716 1.5240 1.6933 1.8814 Columns 8 through 112.0904 2.3227 2.5807 2.8674改進(jìn)歐拉公式結(jié)果

13、： 輸入：x,y=Gaijineuler(f,0,1,1,0.1)結(jié)果：x = 0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1.0000y = 1.0000 1.0959 1.1841 1.2662 1.3434 1.4164 1.4860 1.5525 1.6165 1.6782 1.7379龍格-庫塔計(jì)算結(jié)果：（三階）輸入：Longgekuta3(ff,0,1,1,0.1) 結(jié)果：ans =Longgekuta3(ff,0,1,1,0.1)ans = 0 1.0000 0.1000 1.1048 0.20

14、00 1.2188 0.3000 1.3411 0.4000 1.4711 0.5000 1.6082 0.6000 1.7520 0.7000 1.9021 0.8000 2.0580 0.9000 2.2195 1.0000 2.3863（四階）輸入：Longgekuta4(ff,0,1,1,0.1) 結(jié)果：ans = 0 1.0000 0.1000 1.1048 0.2000 1.2188 0.3000 1.3411 0.4000 1.4711 0.5000 1.6082 0.6000 1.7520 0.7000 1.9021 0.8000 2.0580 0.9000 2.2195 1.

15、0000 2.3863五計(jì)算結(jié)果分析:方法顯示歐拉簡單精度低隱式歐拉穩(wěn)定性最好精度低，計(jì)算量大梯形公式精度提高計(jì)算量大中點(diǎn)公式精度提高，顯式多一個(gè)初值，可能影響精度1.收斂性：若某算法對于任意固定的 x = xi = x0 + i h，當(dāng) h0 ( 同時(shí) i ) 時(shí)有 yi y( xi )，則稱該算法是收斂的。以下討論的都是單步法（指在計(jì)算時(shí)只用到它前一步的信息）歐拉法，龍格-庫塔法都是單步法的例子；例：就初值問題 考察歐拉顯式格式的收斂性。解：該問題的精確解為歐拉公式為 對于固定的 x = xi = x0 + i h，有 設(shè)單步法具有p階精度，滿足利普希茨條件，其整體的截?cái)嗾`差：判斷單步法的

16、收斂性，歸結(jié)為驗(yàn)證增量函數(shù)能否滿足利普希茨條件。對于歐拉法，由于其增量函數(shù)就是f（x，y），因此當(dāng)f（x，y）關(guān)于y滿足利普希茨條件時(shí)候它就是收斂的。對于改進(jìn)歐拉法和龍格-庫塔方法，找到增量函數(shù)，滿足利普希茨條件的時(shí)候，找到利普希茨常數(shù)，這些方法都是收斂的。 2.穩(wěn)定性： 例：考察初值問題在區(qū)間0,0.5上的解，分別用歐拉顯，隱式格式，改進(jìn)歐拉格式計(jì)算數(shù)值解：節(jié)點(diǎn)歐拉顯式歐拉隱式改進(jìn)歐拉法精確解01.00001.00001.00001.00000.1-2.00002.5*10(-1)2.50004.9787*10(-2)0.24.00006.2500*10(-2)6.25002.4788*10

17、(-3)0.3-8.00001.5625*10(-2)1.5625*10(1)1.2341*10(-4)0.41.6000*10(1)3.9063*10(-3)3.9063*10(1)6.1442*10(-6)0.5-3.2000*10(1)9.7656*10(-1)9.7656*10(1)3.0590*10(-7)若某算法在計(jì)算過程中任一步產(chǎn)生的誤差在以后的計(jì)算中都逐步衰減，則稱該算法是絕對穩(wěn)定的.考察顯示歐拉法：，由此可見，要保證初始誤差以后逐步衰減必須滿足：考察后退歐拉法：由此可見，據(jù)對絕對穩(wěn)定絕對穩(wěn)定絕對穩(wěn)定的區(qū)域?yàn)樘菪喂剑耗Ｐ停^對穩(wěn)定區(qū)域?yàn)橄鄬Ψ€(wěn)定區(qū)間為改進(jìn)歐拉公式：相對穩(wěn)定區(qū)域

18、為穩(wěn)定區(qū)間為四階龍格-庫塔：相對穩(wěn)定區(qū)域：相對穩(wěn)定區(qū)間：考察初值問題在區(qū)間0,1上的解節(jié)點(diǎn)改進(jìn)歐拉法四節(jié)龍格-庫塔法準(zhǔn)確解01110.21.1866671.1832291.1832160.41.3483121.3416671.3416410.61.4937041.4832811.4832400.81.6278611.6125141.61245211.7542051.7321421.732051 變步長的龍格-庫塔法：實(shí)際計(jì)算中，步長h的選擇是一個(gè)十分重要的問題，若步長h取的太小，將增加許多不必要的的計(jì)算量，若步長取的太大，其計(jì)算競速很難保證，因此，與數(shù)值積分一樣，微分方程的數(shù)值解也有一個(gè)選擇步

19、長的問題。由四階龍格-庫塔法的性質(zhì)，其局部截?cái)嗾`差為這里的c與在內(nèi)的值有關(guān)系。然后將步長折半，即取為步長，從跨兩步到，在求一個(gè)近似值，每跨一步的誤差為，因此，比較之后誤差的大約減少。六計(jì)算時(shí)出現(xiàn)的問題，解決方法及體會(huì)： 對于初值問題，在用不同方法計(jì)算離散點(diǎn)的函數(shù)值的時(shí)候，有精確值給定的情況下，不同方法跟精確值的差別很大，也就是不是一種好方法，這與實(shí)際問題的穩(wěn)定性有很大的關(guān)系，只要所求的初值問題，能夠落在使用方法的穩(wěn)定區(qū)間內(nèi)，所得到的結(jié)果和精確值的離散點(diǎn)的函數(shù)值相差不會(huì)很大，并且在穩(wěn)定區(qū)間內(nèi)，改進(jìn)歐拉法優(yōu)于歐拉法，三階之后的龍格-庫塔方法優(yōu)于改進(jìn)歐拉法，離散點(diǎn)的函數(shù)值更加接近精確值。所謂的歐拉法

20、就是一階的龍格-庫塔方法，改進(jìn)歐拉法也是在顯示歐拉法預(yù)測和梯形公式校正的基礎(chǔ)上作出的改進(jìn)，也是二階龍格-庫塔的特殊情況。所謂的三階，四階甚至n階龍格-庫塔方法，也是在待定系數(shù)有無數(shù)解得情況下的一般情況。 一階常微分初值問題的收斂性，只要其增量函數(shù)關(guān)于y能滿足利普希茨條件，這就收斂。而起在計(jì)算過程過有避免不了的誤差，歐拉法精度相對來說比較低，龍格-庫塔的階數(shù)越高，精度就越高，所得的結(jié)果和精確值越接近，截?cái)嗾`差越小。同一個(gè)問題，想要減少截?cái)嗾`差，我們可以采用變步長的的方法，將步長變?yōu)?，這樣得出的結(jié)果，根據(jù)不同的方法，截?cái)嗾`差會(huì)有相應(yīng)的改變。而在每種方法的程序設(shè)計(jì)方面，主要思想還是根據(jù)每種方法的公式

21、，找到步長，分點(diǎn)數(shù)，用循環(huán)語句，得出相應(yīng)的結(jié)果，只是在龍格-庫塔方法里面的三階，四階的公式里面的待定系數(shù)有無數(shù)個(gè)解，我們給出的也是一般情況的三階，四階龍格-庫塔方法，精度，誤差相對于歐拉法，改進(jìn)歐拉法，梯形公式都優(yōu)于他們。數(shù)值計(jì)算實(shí)踐課程設(shè)計(jì)原始數(shù)據(jù)記錄實(shí)驗(yàn)名稱：常微分方程初值問題數(shù)值解法研究 實(shí)驗(yàn)時(shí)間： 2013 年 1 月 14 日 輸入：定義M文件并保存ffx.m Euler1(ffx,0,1,1,10)結(jié)果：T = 0 1.0000 0.1000 1.1000 0.2000 1.2100 0.3000 1.3310 0.4000 1.4641 0.5000 1.6105 0.6000

22、1.7716 0.7000 1.9487 0.8000 2.1436 0.9000 2.3579 1.0000 2.5937ans = Columns 1 through 7 1.0000 1.1000 1.2100 1.3310 1.4641 1.6105 1.7716 Columns 8 through 11 1.9487 2.1436 2.3579 2.5937向前歐拉公式結(jié)果：輸入：Euler2(ffx,0,1,1,10)結(jié)果T = 0 1.0000 0.1000 1.1110 0.2000 1.2344 0.3000 1.3716 0.4000 1.5240 0.5000 1.6933 0.6000 1.8814 0.7000 2.0904 0.8000 2.3227 0.9000 2.5807 1.0000 2.8674ans = Columns 1 through 7 1.0000 1.1110 1.2344 1.3716 1.5240 1.6933 1.8814 Columns 8 through 112.0904 2.3227 2.5807 2.8674改進(jìn)歐拉公式結(jié)果： 輸入：x,y=Gaijineuler(f,0,1,1,0.1)結(jié)果：x = 0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8