數(shù)值分析方法面向“四新”人才培養(yǎng)普通高等教育系列教材第六章常微分方程數(shù)值解法目錄/Contents6.1

認(rèn)識微分方程

6.2

微分方程初值問題的Euler方法

6.3

微分方程初值問題數(shù)值解的誤差與穩(wěn)定性分析

6.4

微分方程初值問題的Runge-Kutta法

6.5

非線性微分方程組初值問題的龍格-庫塔法

目錄/Contents6.6

線性多步方法

6.7

微分方程組的剛性問題

6.8二階微分方程的邊值問題

6.9

微分方程計算機(jī)實驗

含未知函數(shù)及其導(dǎo)數(shù)的方程叫做常微分方程如果未知函數(shù)是多元函數(shù)且方程中含有未知函數(shù)的偏導(dǎo)數(shù)則稱其為偏微分方程常微分方程中未知函數(shù)導(dǎo)數(shù)的階數(shù)稱為微分方程的階數(shù)含有任意常數(shù)的函數(shù)滿足微分方程且任意常數(shù)的個數(shù)與微分方程的階數(shù)相等稱為該常微分方程的通解既滿足微分方程且適合初值條件的解稱為該微分方程初值問題的特解6.1

認(rèn)識微分方程6.1.1微分方程模型舉例

因此求解滿足

解的問題稱為微分方程的初值問題。例6.1.2振動模型

介質(zhì)中質(zhì)量為m的質(zhì)點(diǎn)，假定處在彈性約束之下作一維振動(即僅需一個位置參數(shù)就可完全描述質(zhì)點(diǎn)狀態(tài)的運(yùn)動)，通常以彈簧作為這類一維彈性振動的代表模型（圖6-1）。已知質(zhì)點(diǎn)在介質(zhì)中運(yùn)動所受阻力與質(zhì)點(diǎn)速度成正比，根據(jù)牛頓第二定律，

SIS模型只把人群分為S類和I類，患者病愈后并未產(chǎn)生免疫力，仍未S類成員。SIS模型中各類人員的相互作用關(guān)系如圖6-2所示。NSNI

SIR模型考慮S類、I類和R類人群，且I類人病愈后具有免疫能力，不會被傳染而進(jìn)入I類。S，I，R三類人口的關(guān)系如圖6-3所示。NSNINR

方程組（6.1.5）常稱為SIR模型。例6.1.4洛倫茲方程

考慮底部受熱的均勻厚度的水平流體層。設(shè)底部溫度高于頂部溫度，且流體受到重力作用。當(dāng)?shù)撞苛黧w受熱膨脹時，下面流體密度小于上面流體密度，于是在重力場的作用下，上面流體有下沉的趨勢，而下面流體有上浮的趨勢。當(dāng)溫度梯度較小時，由于流體粘性的阻礙，流體仍處于靜止?fàn)顟B(tài)，在流體內(nèi)部只有熱傳導(dǎo)，溫度隨高度呈線性變化。然而當(dāng)溫度梯度較大時，流體的靜止平衡狀態(tài)失穩(wěn)，出現(xiàn)對流現(xiàn)象。這一流體的運(yùn)動可以用一個復(fù)雜的偏微分方程表示。

6.1.2微分方程數(shù)值解一般地，對于一階常微分方程的初值問題

（6.1.8）對于一個由多個變量構(gòu)成微分方程系統(tǒng)，或稱微分方程組，本質(zhì)上解法是一樣的。一般地，一個由n個方程組成的方程組，具有如下形式

謝謝數(shù)值分析方法主編

李冬果李林高磊首都醫(yī)科大學(xué)生物醫(yī)學(xué)工程學(xué)院智能醫(yī)學(xué)工程學(xué)學(xué)系面向“四新”人才培養(yǎng)普通高等教育系列教材第六章常微分方程數(shù)值解法目錄/Contents6.1

認(rèn)識微分方程

6.2

微分方程初值問題的Euler方法

6.3

微分方程初值問題數(shù)值解的誤差與穩(wěn)定性分析

6.4

微分方程初值問題的Runge-Kutta法

6.5

非線性微分方程組初值問題的龍格-庫塔法

目錄/Contents6.6

線性多步方法

6.7

微分方程組的剛性問題

6.8二階微分方程的邊值問題

6.9

微分方程計算機(jī)實驗

6.2微分方程初值問題的Euler方法6.2.1顯式Euler公式

xEuler法y精確解01.000001.000000.11.100001.089960.21.180941.15940.31.242271.207070.41.282771.230960.51.300611.228220.61.293251.194970.71.257361.125890.81.188541.013400.91.080820.8457811.00.9258020.6016496.2.2隱式Euler公式與改進(jìn)Euler公式

例6.2.2

用改進(jìn)Euler公式（6.2.10）求解初值問題（6.2.3）。xEuler法y精確解改進(jìn)的Euler法y01.000001.0000010.11.100001.089961.090470.21.180941.15941.160470.31.242271.207071.20880.41.282771.230961.233480.51.300611.228221.231740.61.293251.194971.199760.71.257361.125891.132360.81.188541.013401.022160.91.080820.8457810.8578951.00.9258020.6016490.619366謝謝數(shù)值分析方法主編

李冬果李林高磊首都醫(yī)科大學(xué)生物醫(yī)學(xué)工程學(xué)院智能醫(yī)學(xué)工程學(xué)學(xué)系面向“四新”人才培養(yǎng)普通高等教育系列教材第六章常微分方程數(shù)值解法目錄/Contents6.1

認(rèn)識微分方程

6.2

微分方程初值問題的Euler方法

6.3

微分方程初值問題數(shù)值解的誤差與穩(wěn)定性分析

6.4

微分方程初值問題的Runge-Kutta法

6.5

非線性微分方程組初值問題的龍格-庫塔法

目錄/Contents6.6

線性多步方法

6.7

微分方程組的剛性問題

6.8二階微分方程的邊值問題

6.9

微分方程計算機(jī)實驗

6.3初值問題數(shù)值解的誤差與穩(wěn)定性分析6.3.1誤差分析

（6.3.1）

6.3.2收斂性與穩(wěn)定性分析1.數(shù)值方法的收斂性

2.數(shù)值方法的穩(wěn)定性

謝謝數(shù)值分析方法主編

李冬果李林高磊首都醫(yī)科大學(xué)生物醫(yī)學(xué)工程學(xué)院智能醫(yī)學(xué)工程學(xué)學(xué)系面向“四新”人才培養(yǎng)普通高等教育系列教材第六章常微分方程數(shù)值解法目錄/Contents6.1

認(rèn)識微分方程

6.2

微分方程初值問題的Euler方法

6.3

微分方程初值問題數(shù)值解的誤差與穩(wěn)定性分析

6.4

微分方程初值問題的Runge-Kutta法

6.5

非線性微分方程組初值問題的龍格-庫塔法

目錄/Contents6.6

線性多步方法

6.7

微分方程組的剛性問題

6.8二階微分方程的邊值問題

6.9

微分方程計算機(jī)實驗

6.4.微分方程初值問題的Runge-Kutta法6.4.1Runge-Kutta法的基本思想與二階Runge-Kutta法

例6.4.1

用二階Runge-Kutta公式（6.4.2）求解初值問題（6.2.3）。

二階Runge-Kutta公式（6.4.2）所求解與改進(jìn)Euler方法所求解具有基本一致的誤差，圖6-4（右）為兩種算法與精確解比較的結(jié)果，其中菱形與圓點(diǎn)數(shù)據(jù)分別表示相同節(jié)點(diǎn)處精確解分別與二階Runge-Kutta公式（6.4.2）所求解和改進(jìn)Euler方法所求解之差。可以看出二種算法有微小差別。6.4.2三、四階Runge-Kutta法

計算的結(jié)果見圖6-5，其中實心方塊為精確解，實心三角為三階Runge-Kutta公式，空心三角為二階Runge-Kutta公式所求結(jié)果，空心方塊為改進(jìn)Euler方法所求解。

6.4.3隱式Runge-Kutta法

謝謝數(shù)值分析方法主編

李冬果李林高磊首都醫(yī)科大學(xué)生物醫(yī)學(xué)工程學(xué)院智能醫(yī)學(xué)工程學(xué)學(xué)系面向“四新”人才培養(yǎng)普通高等教育系列教材第六章常微分方程數(shù)值解法目錄/Contents6.1

認(rèn)識微分方程

6.2

微分方程初值問題的Euler方法

6.3

微分方程初值問題數(shù)值解的誤差與穩(wěn)定性分析

6.4

微分方程初值問題的Runge-Kutta法

6.5

非線性微分方程組初值問題的龍格-庫塔法

目錄/Contents6.6

線性多步方法

6.7

微分方程組的剛性問題

6.8二階微分方程的邊值問題

6.9

微分方程計算機(jī)實驗

6.5非線性微分方程組初值問題的

龍格-庫塔法

謝謝數(shù)值分析方法主編

李冬果李林高磊首都醫(yī)科大學(xué)生物醫(yī)學(xué)工程學(xué)院智能醫(yī)學(xué)工程學(xué)學(xué)系面向“四新”人才培養(yǎng)普通高等教育系列教材第六章常微分方程數(shù)值解法目錄/Contents6.1

認(rèn)識微分方程

6.2

微分方程初值問題的Euler方法

6.3

微分方程初值問題數(shù)值解的誤差與穩(wěn)定性分析

6.4

微分方程初值問題的Runge-Kutta法

6.5

非線性微分方程組初值問題的龍格-庫塔法

目錄/Contents6.6

線性多步方法

6.7

微分方程組的剛性問題

6.8二階微分方程的邊值問題

6.9

微分方程計算機(jī)實驗

6.6線性多步方法

6.6.1線性多步方法的構(gòu)造

（6.6.1）

6.6.2線性多步方法的應(yīng)用及預(yù)測-校正方法1.線性多步方法的應(yīng)用

2.預(yù)測-校正方法

2.線性多步方法的收斂性與穩(wěn)定性

謝謝數(shù)值分析方法主編

李冬果李林高磊首都醫(yī)科大學(xué)生物醫(yī)學(xué)工程學(xué)院智能醫(yī)學(xué)工程學(xué)學(xué)系面向“四新”人才培養(yǎng)普通高等教育系列教材第六章常微分方程數(shù)值解法目錄/Contents6.1

認(rèn)識微分方程

6.2

微分方程初值問題的Euler方法

6.3

微分方程初值問題數(shù)值解的誤差與穩(wěn)定性分析

6.4

微分方程初值問題的Runge-Kutta法

6.5

非線性微分方程組初值問題的龍格-庫塔法

目錄/Contents6.6

線性多步方法

6.7

微分方程組的剛性問題

6.8二階微分方程的邊值問題

6.9

微分方程計算機(jī)實驗

6.7微分方程組的剛性問題

求解微分方程組時，經(jīng)常出現(xiàn)解的分量數(shù)量級差別很大的情形，這給數(shù)值求解帶來了較大困難，這種問題稱為剛性問題。

步數(shù)nxyz10202-7.0354920.1041-6.828563-261.236134014.-8989.144-7.13065*10^147.47272*10^26-3.39266*10^255-4.98155*10^1044.42368*10^172-2.00838*10^171

謝謝數(shù)值分析方法主編

李冬果李林高磊首都醫(yī)科大學(xué)生物醫(yī)學(xué)工程學(xué)院智能醫(yī)學(xué)工程學(xué)學(xué)系面向“四新”人才培養(yǎng)普通高等教育系列教材第六章常微分方程數(shù)值解法目錄/Contents6.1

認(rèn)識微分方程

6.2

微分方程初值問題的Euler方法

6.3

微分方程初值問題數(shù)值解的誤差與穩(wěn)定性分析

6.4

微分方程初值問題的Runge-Kutta法

6.5

非線性微分方程組初值問題的龍格-庫塔法

目錄/Contents6.6

線性多步方法

6.7

微分方程組的剛性問題

6.8二階微分方程的邊值問題

6.9

微分方程計算機(jī)實驗

6.8二階微分方程的邊值問題二階微分方程也可以看作是微分方程組，如（6.1.2）可以化作（6.1.10），可以利用微分方程組的數(shù)值方法求解實際應(yīng)用中較為常見的微分方程的初值問題，但對二階微分方程而言，所謂的邊值問題在科學(xué)技術(shù)問題中也常常出現(xiàn)，其求解問題與上述初值問題不同，本節(jié)簡單討論兩種本方法。6.8.1二階微分方程邊值問題的打靶法

6.8.2二階線性微分方程邊值問題的差分法

固體表面擴(kuò)散層液體生物膜

謝謝數(shù)值分析方法主編

李冬果李林高磊首都醫(yī)科大學(xué)生物醫(yī)學(xué)工程學(xué)院智能醫(yī)學(xué)工程學(xué)學(xué)系面向“四新”人才培養(yǎng)普通高等教育系列教材第六章常微分方程數(shù)值解法目錄/Contents6.1

認(rèn)識微分方程

6.2

微分方程初值問題的Euler方法

6.3

微分方程初值問題數(shù)值解的誤差與穩(wěn)定性分析

6.4

微分方程初值問題的Runge-Kutta法

6.5

非線性微分方程組初值問題的龍格-庫塔法

目錄/Contents6.6

線性多步方法

6.7

微分方程組的剛性問題

6.8二階微分方程的邊值問題

6.9

微分方程計算機(jī)實驗

6.9微分方程計算機(jī)實驗6.9.1顯式Euler公式和改進(jìn)Euler公式的實現(xiàn)我們可以通過以下代碼實現(xiàn)Euler法：defeuler(f,x0,y0,b,h):

y=[y0]

x=[x0]

while(x0+h)<=b:

y0=h*f(x0,y0)+y0

x0=x0+h

y.append(y0)

x.append(x0)

returny,x其中輸入?yún)?shù):函數(shù)f即為所求的初值問題中的導(dǎo)函數(shù)，x0為初值點(diǎn)，y0為初值點(diǎn)處對應(yīng)的函數(shù)值，b為積分區(qū)間的終點(diǎn)，而h則是步長。返回值為兩個列表，其中分別存儲計算出的微分方程解的函數(shù)值，和與該函數(shù)值對應(yīng)的自變量值。類似于Euler公式，將遞推公式進(jìn)行修改，就得到了改進(jìn)Euler法的求解函數(shù)defeuler_improved(f,x0,y0,b,h):

y=[y0]

x=[x0]

while(x0+h)<=b:

x1=x0+h

yp=y0+h*f(x0,y0)

yc=y0+h*f(x1,yp)

y0=(yp+yc)/2

x0=x1

y.append(y0)

x.append(x0)

returny,x其參數(shù)與之前Euler法的定義相同。例6.9.1用Euler公式（6.2.2）和改進(jìn)Euler公式（6.2.10）求解初值問題（6.2.3）。解：

為了實現(xiàn)指數(shù)運(yùn)算，需要導(dǎo)入其他工具包，這里使用numpy工具包。importnumpyasnp使用lambda定義導(dǎo)函數(shù)fun_6_2_1和解函數(shù)real_ans_6_2_1fun_6_2_1=lambdax,y:y-3*x/y**(1/3)

real_ans_6_2_1=lambdax:(9+12*x-5*np.exp(4*x/3))**(3/4)/2**(3/2)應(yīng)用之前定義的函數(shù)實現(xiàn)Euler公式和改進(jìn)Euler公式，將函數(shù)值分別返回給y1和y2[y1,x1]=euler(fun_6_2_1,0,1,1,0.1)

[y2,x2]=euler_improved(fun_6_2_1,0,1,1,0.1)計算真實的函數(shù)值y_real：

y_real=[real_ans_6_2_1(x)forxinx1]

打印結(jié)果：print("x\tEuler法y\t真實值\t改進(jìn)Euler法y")

forkinrange(11):

print(f"{x1[k]:.1f}\t{y1[k]:.6f}\t{y_real[k]:.6f}\t{y2[k]:.6f}")結(jié)果見表6-7。6.9.2四階Runge-Kutta法的實現(xiàn)類似于之前Euler和改進(jìn)Euler法的實現(xiàn)，可以根據(jù)公式（6.4.7）實現(xiàn)四階的Runge-Kutta法，可以如下定義函數(shù)defrk4(f,x0,y0,b,h):

x=[x0]

y=[y0]

while(x0+h)<=b:

x1=x0+h

k1=f(x0,y0)

k2=f(x0+h/2,y0+h/2*k1)

k3=f(x0+h/2,y0+h/2*k2)

k4=f(x1,y0+h*k3)

y0=y0+h*(k1+2*k2+2*k3+k4)/6

x0=x1

y.append(y0)

x.append(x0)

returny,x其參數(shù)與之前Euler法的定義相同。

可以使用matplotlib工具包來繪制函數(shù)圖像，為了方便數(shù)據(jù)處理，同時引入numpy工具包importnumpyasnp

importmatplotlib.pyplotasplt

x1_array=np.array(x1)

Y_array=np.array(Y)

markers=["s","o","s","^"]

markerfillstyles=["none","full","full","full"]

plt.figure(figsize=(10,6))

forkinrange(4):

plt.plot(x1_array[:200:5],Y_array[k,:200:5],

marker=markers[k],fillstyle=markerfillstyles[k],

linestyle='',color='blue')

plt.axhline(5,0,10,linestyle='-')

plt.axhline(2,0,10,linestyle='--')

plt.show()結(jié)果見圖6-6。6.9.3方程組的四階Runge-Kutta法實現(xiàn)根據(jù)公式（6.5.2），可以得到方程組的數(shù)值計算方法，可以編寫Python函數(shù)如下。求解由2個方程組成的方程組的函數(shù)rk4_2:defrk4_2(f,g,t0,x0,y0,b,h):

t=[t0]

x=[x0]

y=[y0]

while(t0+h)<=b:

t1=t0+h

k1=f(t0,x0,y0)

l1=g(t0,x0,y0)

k2=f(t0+h/2,x0+h/2*k1,y0+h/2*l1)

l2=g(t0+h/2,x0+h/2*k1,y0+h/2*l1)

k3=f(t0+h/2,x0+h/2*k2,y0+h/2*l2)

l3=g(t0+h/2,x0+h/2*k2,y0+h/2*l2)

k4=f(t1,x0+h*k3,y0+h*l3)

l4=g(t1,x0+h*k3,y0+h*l3)

x0=x0+h*(k1+2*k2+2*k3+k4)/6

y0=y0+h*(l1+2*l2+2*l3+l4)/6

t0=t1

x.append(x0)

y.append(y0)

t.append(t0)

returnx,y,t

同理可以計算無阻尼假設(shè)下的運(yùn)動軌跡：g2=lambdat