常微分方程算法的應用（楊牧原20096313）常微分方程算法背景科學研究和工程技術中許多問題在數學上往往歸結為微分方程求解問題。為了確定微分方程的解，一般要加上定解條件，根據不同的情況，這些定解條件主要有初始條件和邊界條件。只含初值條件作為定解條件的微分方程求解問題稱為初值問題；例如在天文學中研究星體運動，空間技術中研究物體飛行等，都需要求解常微分方程初值問題。只含邊界條件作為定解條件的微分方程求解問題稱為邊值問題。除特殊情形外，微分方程一般求不出解析解，即使有的能求出解析解，其函數表示也比較復雜，計算量比較大，而且實際問題往往只要求在某一時刻解的函數值。為了解決這個問題，有兩種方法可以逼近原方程的解。第一種方法是：將原微分方程化簡為可以準確求解的微分方程，然后使用化簡后的方程解近似原方程的解。第二種方法是：將求原微分方程的解析解，轉化成求原微分方程的數值解，這是實際中最常用的方法。本課程設計所指常微分方程算法，即求微分方程解初值問題解的數值方法。對于微分方程解初值問題解的數值方法，常用的有歐拉方法、改進的歐拉法、龍格庫塔方法、單步法、線性多步法等。二、數值解的一般概念常微分方程初值問題的數值解是指通過一定的近似方法得出準確解在一列離散點上的近似值。數值解的特征是步進式，即在點的近似值是由等若干點處的近似值的信息給出的遞推公式。若依賴于前面步的值，則稱為步法；稱為單步法。利用在的精確解借助某種算法計算出，則稱為該方法的局部截斷誤差。如果一個算法的局部截斷誤差是，則稱該方法是階的；而利用數值解得到的與微分方程的精確解之差稱為整體截斷誤差，即是該數值方法的誤差。對于固定的，取，用某種算法得到，如有=0，則稱該方法是收斂的。注意，因是固定的，隨著，數值解的步數。在實際計算時由于舍入誤差不可避免，實際得到數值解是，穩定性即研究是否隨著計算步驟的增加而增加。通常所提的穩定性是通過模型方程來討論的。若當某一步有舍入誤差時，在以后的計算中誤差不會逐步擴大，則稱這種穩定性為絕對穩定性。三、常用的數值方法（1）Euler法：Euler法的局部截斷誤差為，整體截斷誤差為，即一階收斂。對于模型問題，當時，Euler法是數值穩定的。隱式Euler法的誤差與Euler法相同，但是無條件穩定：即對任意步長，隱式Euler法都是穩定的。梯形法的誤差比Euler法高一階，也是無條件穩定的。改進Euler法是一種預測－校正方法：——Euler法預測——梯形法校正它保持了梯形法的誤差階數，但不是無條件穩定的。龍格-庫塔方法龍格-庫塔類算法采用區間內若干點的斜率的加權平均來近似整個區間的平均斜率，一般形式為如經典的4級4階（局部截斷誤差為）Runge-Kutta公式為（3）線性多步法線性步法具有一般形式（2）為隱式公式；為顯式公式。構造多步法公式有基于數值積分和Taylor展開兩種途徑。多步法（2）的局部截斷誤差為利用原微分方程后，成為因此利用Taylor公式，分別對和作Taylor展開，可確定線性多步（2）中的待定參數和，使它達到最高階精度（或指定精度）。預測－校正格式：不論單步法還是多步法，隱式公式比顯式公式的穩定性好，但隱式公式的計算比較困難。預測－校正格式是用顯式公式進行預測，再用隱式公式進行校正。常微分方程算法的MATLAB實現本課程設計采用兩個常微分方程實例，用MATLAB實現其數值解法。即用歐拉公式和四階龍格-庫塔法分別求解下列初值問題；運行結果：取步長h=0.1，計算（1）、（2）用歐拉公式計算的（1）、（2）從y0到y10分別為：y1=1.0000000000000001.0750000000000001.1441071428571431.2084631696428571.2688863281250001.3259862128906251.3802311034179691.4319897697961431.4815586464429321.5291801743643121.575055579595242y2=2.0000000000000002.1980198019801982.6207159177456213.3420138767581774.4944324549506526.2922054369309139.06817842381219913.32839647526759619.83005329247130229.69030078596532044.535451178947980用四階龍格-庫塔法公式計算的（1）、（2）從y0到y10分別為：y3=1.0000000000000001.0718306860207101.1382106030709361.2001659371598141.2584390526264251.3135871986541281.3660411042607671.4161417295045441.4641643601444171.5103349913532301.554841814242189y4=2.0000000000000002.0294857261237892.0776988952688942.1433760827529662.2249712844066192.3208061549556472.4291959592761182.5485393147928252.6773718703260642.8143907897117602.958458858618671五、結論算法分析：取步長為h=0.1，分別用歐拉及四階龍格-庫塔方法進行計算。歐拉初值公式為：y.n+1=y.n+h*f(t.n-1,y.n),y.0=a;四階龍格-庫塔公式為:k1=f(tn,yn),k2=f(tn+h/2,yn+h*k1/2),k3=f(tn+h/2,yn+h*k2/2),k4=f(tn+h,yn+h*k3),y.n+1=y.n+(h/6)*(k1+2*k2+2*k3+k4).從上面的運行結果可知，四階龍格庫塔方法具有高得多的精確度。但需要指出的是，龍格庫塔方法是基于泰勒展開法，因而求方程的解具有足夠的光滑性。如果解的光滑性差，使用四階龍格庫塔方法求得數值解的精確度，可能不如改進的歐拉方法精確度高。因此，在實際的計算中，要根據具體問題的特性，選擇合適的算法。六、參考文獻[1]朱曉臨主編《數值分析》中國科學技術大學出版社，2010[2]關治，陸金甫編著《數值分析基礎》高等教育出版社，2008附錄：源程序disp('取步長h=0.1，計算（1）、（2）')y1=zeros(11,1);y2=zeros(11,1);y1(1)=1;y2(1)=2;k=1;fori=1:10y1(k+1)=y1(k)+0.1*(0.9*y1(k)/(1+2*k*0.1));y2(k+1)=y2(k)+(k*0.1)*y2(k)/(1+(k*0.1)*(k*0.1));k=k+1;enddisp('用歐拉公式計算的（1）、（2）從y0到y10分別為：')y1y2y3=zeros(11,1);y4=zeros(11,1);y3(1)=1;y4(1)=2;k=1;h=0.1;fori=1:10k1=0.9*y3(k)/(1+2*h*k);k2=0.9*(y3(k)+h*k1/2)/(1+2*(h*k+h/2));k3=0.9*(y3(k)+h*k2/2)/(1+2*(h*k+h/2));k4=0.9*(y3(k)+h*k3)/(1+2*(h*k+h));y3(k+1)=y3(k)+h*(k1+2*k2+2*k3+k4)/6;m1=(k*h)*y4(k)/(1+(k*h*k*h));m2=(h*k+h/2)*(y4(k)+h*m1/2)/(1+(h*k+h/2)*(h*k+h/2));m3=(h*k+h/2)*(y4