第二章常微分方程的數值求解方法1、引言2、數值求解的基本思想3、歐拉方法4、R-K方法5、線性多步法6、預估校正方法7、步長選擇8、常微分方程組和高階微分方程9、邊值問題10、有限差分方法內容提綱1、引言微分方程：定義：包含自變量、未知函數及未知函數的導數或微分的方程稱為微分方程。“常”和“偏”：在微分方程中,自變量的個數只有一個,

稱為常微分方程。自變量的個數為兩個或兩個以上的微分方程叫偏微分方程。微分方程中出現的未知函數最高階導數的階數稱為微分方程的階數。如果未知函數y及其各階導數。都是一次的，則稱它是線性的，否則稱為非線性的。解析求解?在高等數學中，對于常微分方程的求解，給出了一些典型方程求解析解的基本方法，如可分離變量法、常系數齊次線性方程的解法、常系數非齊次線性方程的解法等。

但能求解的常微分方程仍然是有限的（1%?)。

大多數的常微分方程是不可能給出解析解。譬如：這個一階微分方程就不能用初等函數及其積分來表達它的解。一個具體的例子

解析求解無能為力例：多體運動問題多體運動問題是物理研究中的一個普遍問題。設有n個物體（粒子）在靜電場或引力場中運動，因而有6n個未知數xi(坐標)，vi（速度），（i=1,2,…,n）。根據牛頓定律或庫倫定律，應該滿足6n個一階常微分方程組：其中mi是第i個物體（或粒子）的質量，ei在靜電場中是第i個粒子的電荷，在引力場中事實上，當n≥3時，無法找到解析解。數值求解根據常微分方程的定解條件，常微分方程可分為初值問題和邊值問題。如果定解條件是描述函數在初始點的狀態，則稱為初值問題。一個典型的常微分方程初值問題的形式如下：

具有以下形式的常微分方程及定解條件，如

則為常微分方程的邊值問題。常微分方程初值問題的數值求解

由此可見，常微分方程數值求解的基本出發點就是將求解區間離散化。

3、歐拉（Euler）方法這是一種最為基礎、簡單的方法，它的精確度不高，因此在實際計算中并不常使用。由于其能反映出很多復雜方法的特征，被廣泛用于微分方程解法的講解。根據計算方法不同，分以下三種：歐拉（Euler）方法：根據前一個點的函數值，可以計算后一個點的函數值。

三種簡單的數值解法1、化導數為差商；2、數值積分法；3、Taylor展開

數值積分法

Talyor級數展開

歐拉公式及其幾何意義

歐拉方法的穩定性分析由歐拉公式可以看出，在常微分方程初值問題的數值求解過程中，每一步計算都包含了前一步的計算結果，因此前面的計算誤差將會影響后面的數值計算結果。需要對數值求解方法的誤差進行系統分析，考察計算過程的誤差積累會不會掩蓋真解，這就是常微分方程數值解的穩定性問題。需要對所采用的數值方法進行穩定性分析，確定其穩定性區域，保證該方法得到的數值解是合理的。

總體誤差和局部誤差的關系

由此可得：

收斂性和穩定性是兩個不同的概念。收斂性是指數值方法的截斷誤差對計算結果的影響穩定性指的是某一步的計算誤差對計算結果的影響此外，穩定性與步長密切相關。

數值方法的穩定性

0-1-2ReImg歐拉方法的絕對穩定區域絕對穩定區域越大，這個方法的適應性就越強。若絕對穩定區域包含復平面的整個左半平面，就稱這個數值方法是A穩定的。210ReImg后退歐拉方法的絕對穩定區域：向后歐拉方法是A穩定的，但是仍受到迭代的限制，即0<hL<1。

總體截斷誤差局部截斷誤差穩定性歐拉方法向后歐拉方法A穩定改進歐拉方法A穩定

以上三種方法的誤差和穩定性總結如下：RUNGE-KUTTA方法

該方程組存在無窮多個解。所有滿足該方程組的解都為二階Runge-Kutta公式。

這就是改進歐拉公式。因此，改進歐拉公式也是二階Runge-Kutta公式。

具體推導參考《計算物理學》李茂枝等編著（2024）

四階Runge-Kutta公式的經典形式

一步法在計算時只用到了前面一步的近似值，但要提高精度，需要增加中間函數值的計算，這就加大了計算量。R-K方法的優缺點是否可以尋找一種精度較高、計算量適中的方法呢？注意到：計算yn+1時，yn，yn-1，yn-2等等均已知，那么可否利用這些值提高精度而不顯著加大計算量呢？優點：一步法，在給定初值后可以逐步計算下去；精度高；在計算過程中便于改變步長缺點：計算量大線性多步法

四階Adams方法

Adams內插公式

44內插值迭代求解45Adams預估--矯正公式

步長的自動選擇46小步長高精度4748一個例子49步數和步長步數步長10.0735892100.0397480200.00475879300.00169146400.0123375500.356060600.0280562700.00326283790.000106687800.000120953900.005208841000.1201341060.069488250階數和步長

二分法選取步長52常微分方程組和高階微分方程53

然后采用合適的數值方法求解

邊值問題54在實際問題中也常常碰到所謂的兩點邊值問題

要得到方程（1）的解的存在唯一性非常困難。打靶法（shootingmethod）是數值求解常微分方程兩點邊值問題的常用方法。

打靶法（shootingmethod）55

打靶法的基本思想是將微分方程的邊值問題轉化為一系列初值問題。先考慮一個初值問題，即

打靶法就是將微分方程的邊值問題轉化為初值問題。通過適當選取和調整初值條件，求解一系列初值問題，使之逼近給定的邊界條件。如果將描述的曲線視做彈道，那么求解過程就是不斷調整試射條件使之打到預定的靶子。打靶法的主要思想炮擊中的例子：試射兩發，調整后齊射αβ有限差分方法有限差分法是一種求解微分方程數值解的近似方法，其主要思想是將微分方程中的微分項直接進行差分近似，從而將微分方程轉化為代數方程組求解。隨著計算機技術和數值算法的快速發展，有限差分法已成為數值求解微分方程的重要方法之一，在眾多領域得到了廣泛的應用和發展。在有限差分方法中，我們放棄了微分方程中獨立變量可以取連續值的特征，而關注獨立變量離散取值后對應的函數值。原則上，這種方法仍然可以達到任意滿意的計算精度，因為方程的連續數值解可以通過減小獨立變量離散取值的間隔，或者通過離散點上的函數值插值來近似得到。63有限差分法包含兩部分：1、用差分代替微分方程中的微分，將連續變化的變量離散化，從而得到差分方程組的數學形式。2、求解差分方程組。在第一步中，通過所謂的網絡分割法，將函數定義域分成大量相鄰而不重合的子區域，通常采用的是規則的分割方式。網絡線劃分的交點稱為節點。

具體方法64

將函數定義域離散化后，需要求得特定問題在所有這些節點上的近似值，因此數值求解的關鍵就是要找到適當的數值計算方法。設一個函數在x點上的一階和二階微商，可以近似地用它的近鄰兩點的函數值的差分來表示，即差分形式65一、二階微商66

利用上面幾個公式，可以構造出微分方程的差分形式。

小結1、線性多步法(避免大量迭代）

2、預估校正方法（合理選擇步長）

3、步長事后估計4、兩點邊值問題

（打靶法，求解大量實際問題）5、常微分方程組和高階微分方程

高階方程總可以化作一階方程組

一階方程組的求解與一階方程的求解類似6、有限差分方法—更實用的方法6869求解一個實際問題70R-K公式的經典形式

71Matlab形式K1=h*feval(@fx,x(i),y(i));K2=h*feval(@fx,x(i)+h/2,y(i)+K1/2);K3=h*feval(@fx,x(i)+h/2,y(i)+K2/2);K4=h*feval(@fx,x(i)+h,y(i)+K3);y(i+1)=y(i)+(K1+2*K2+2*K3+K4)/6;72初值和相應參數區間[a,b]

[0,1]

a=0;b=1;步長h=0.2(可調)總共的計算次數N

N=(b-a)/h初始值：

x0=0;y0=1;73N次求解fori=1:N;x(i)=a+(i-1)*h;K1=h*feval(@fx,x(i),y(i));K2=h*feval(@fx,x(i)+h/2,y(i)+K1/2);K3=h*feval(@fx,x(i)+h/2,y(i)+K2/2);K4=h*feval(@fx,x(i)+h,y(i)+K3);y(i+1)=y(i)+(K1+2*K2+2*K3+K4)/6;end74基本大功告成！???Errorusing==>plotVectorsmustbethesamelengths.Errorin==>RKat23plot(x,y);x(N+1)=b;精確解：y=sqrt(1+2.*x);75結果令人滿意76增加求解區間[0,1][0,8]77改變h到0.0278Adams方法的修正方案y(i+1)=c(i+1)-19.*(c(i+1)-p(i+1))./270;p(i+1)=y(i)+h.*(55.*fn-59.*fn_1+37.*fn_2-9.*fn_3)./24;m(i+1)=p(i+1)+251.*(c(i)-p(i))./270;q(i+1)=feval(@fx,x(i+1),m(i+1));c(i+1)=y(i)+h.*(9.*q(i+1)+19.*fn-5.*fn_1+fn_2)./24;79Fn,fn-1,fn-2,fn-3?

根據定義：fn=feval(@fx,x(i),y(i));fn_1=feval(@fx,x(i-1),y(i-1));fn_2=feval(@fx,x(i-2),y(i-2));fn_3=feval(@fx,x(i-3),y(i-3));如何求解前三個點？R-K方法！80R-K方法求解前三個點fori=1:3;x(i)=a+(i-1).*h;K1=h.*feval(@fx,x(i),y(i));K2=h.*feval(@fx,x(i)+h/2,y(i)+K1/2);K3=h.*feval(@fx,x(i)+h/2,y(i)+K2/2);K4=h.*feval(@fx,x(i)+h,y(i)+K3);y(i+1)=y(i)+(K1+2*K2+2*K3+K4)./6;end81就這么簡單么？???Undefinedfunctionormethod'c'forinputargumentsoftype'double'.Errorin==>Adamsat36m(i+1)=p(i+1)+251.*(c(i)-p(i))./270;i=4,p(4)沒有定義

p(4)=0同樣，c(4)也沒有定義c(4)=08283在區間[0,8]上結果如何？84Physicistscandoanything!h=0.01例：BicycleRacingWhatfactorsdeterminetheultimatespeedofabicycle?Howtoestimatethespeedforarealisticcase?Asimplercase:frictionisignoredEquationofmotion:Constantpowermodel:《ComputationalPhysics》N.J.GiordanoandH.Nakanishi,(2007)選取：m=70Kg,v0=4m/s,deltat=0.1s.ConsidertheairresistanceC:dragcoefficientA:frontalareaoftheobjectρ:densityofairIntheabove