用Python數值求解偏微分方程

1 引言

微分方程是描述一個系統的狀態隨時間和空間演化的最基本的數學工具之一，其在物理、經濟、工程、社會等各方面都有及其重要的應用。然而，只有很少的微分方程可以解析求解，尤其對於偏微分方程，能解析求解的種類更是寥寥可數。更多的微分方程可以採用數值法進行求解，只要精度足夠高，就可以滿足科學和工程上的需求。

數值求解微分方程的基本思路是先把時間和空間離散化，然後將微分化為差分，建立遞推關係，然後利用計算機強大的重複計算能力，快速得到任意格點處的值。Python的Numpy、Scipy工具包可以很好地實現此功能，Matplotlib工具包則可以將求解結果畫為非常直觀的圖形。

接下來，我們先以常微分方程的數值求解為例，引入差分的思想，再將其推廣到偏微分方程中。

2 常微分方程的差分求解

一般地，一階常微分方程可以寫為

$\begin{equation} \left\{ \begin{aligned} &\frac{dy}{dt}=f(y,t)\\ &y(t_0)=y_0\\ \end{aligned} \right. \end{equation}$

首先，將連續的變數

和

離散化，連續的函式

和

化為離散的序列

和

，則上述微分方程可以化為差分方程【1】

$\begin{equation} \left\{ \begin{aligned} &\frac{y_{n+1}-y_{n}}{\Delta t}=\frac{f(y_n,t_n)+f(y_{n+1},t_{n+1})}{2}\\ &y(0)=y_0\\ \end{aligned} \right. \end{equation}$

從而我們得到遞推關係

$\begin{equation} \left\{ \begin{aligned} &y_{n+1}=y_n+\frac{\Delta t}{2}(f(y_n,t_n)+f(y_{n+1},t_{n+1}))\\ &y(0)=y_0 \end{aligned} \right. \end{equation}$

有了遞推關係和初始條件以後，就可以利用Python的強大計算功能，得到任意的

的值了，下面我們透過一個具體的例子來說明。

2.1 RC迴路放電問題

對於一個

迴路，我們有

$\begin{equation} I R+\frac{Q}{C}=0 \end{equation}$

其中，

分別為電流，電阻，電量和電容，利用

，並定義

$\tau\equiv RC$

，我們得到一個含初始條件的一階常微分方程

$\begin{equation} \left\{ \begin{aligned} &\frac{dQ}{dt}=-\frac{Q}{\tau}\\ &Q(t_0)=Q_0 \end{aligned} \right. \end{equation}$

這個方程當然可以解析求解，得到

$Q=Q_0 e^{-t/\tau}$

。另一方面按照差分法，可以得到遞推關係

$\begin{equation} \left\{ \begin{aligned} &Q_{n+1}=Q_n-\frac{1}{2}(\frac{Q_n}{\tau}+\frac{Q_{n+1}}{\tau})\Delta t\\ &Q(0)=Q_0 \end{aligned} \right. \end{equation}$

下面我們用Python進行數值求解，並和解析結果進行比較。

import

numpy

import

matplotlib。pyplot

plt

2。0

#設定常數

0。5

#設定步長

1000

#設定分割段數

0。0

#設定初始時間

1。0

#設定初始電量

#先定義三個空列表

［］

#用來盛放差分得到的q值

qt0

［］

#用來盛放解析得到的q值

time

［］

#用來盛放時間值

for

range

（

）：

#qn+1的近似值

0。5

（

）

#差分遞推關係

。

exp

（

）

#解析關係

。

append

（

）

#差分得到的q值列表

qt0

。

append

（

）

#解析得到的q值列表

time

。

append

（

）

#時間列表

plt

。

plot

（

time

，

‘o’

，

label

‘Euler-Modify’

）

#差分得到的電量隨時間的變化

plt

。

plot

（

time

，

qt0

，

‘r-’

，

label

‘Analytical’

）

#解析得到的電量隨時間的變化

plt

。

xlabel

（

‘time’

）

plt

。

ylabel

（

‘charge’

）

plt

。

xlim

（

，

）

plt

。

ylim

（

0。2

，

1。0

）

plt

。

legend

（

loc

‘upper right’

）

plt

。

show

（）

圖1：用改進的尤拉法差分得到的結果和解析結果的比較

在圖1中給出了差分法得到的結果與解析法得到結果的比較，發現兩者符合得很好，這說明對於這個問題，改進的尤拉法已經可以給出足夠精確的結果。需要注意的是，這個微分方程本身比較簡單，可以解析求解，而對於複雜得多的微分方程，沒法解析求解，但是上述數值求解差分方法仍然是適用的。

3 偏微分方程的差分求解

有了差分代替微分的思想，接下來我們將其推廣到偏微分方程的求解中。以一般二階拋物型偏微分方程為例，一般的可以寫為

$\begin{equation} \left\{ \begin{aligned} &\frac{\partial u(x,t)}{\partial t}=\lambda \frac{\partial ^2 u(x,t) }{\partial x^2} \qquad 0\leq t\leq T, 0\leq x \leq l \\ &u(x,0)=f(x)\\ &u(0,t)=g_1(t)\\ &u(l,t)=g_2(t) \end{aligned} \right. \end{equation}$

仍然是將時間和空間離散化，將微分化為差分，即

$\begin{equation} \left\{ \begin{aligned} &\frac{\partial^2 u(x,t)}{\partial x^2}|_{i,k}=\frac{u_{i-1,k}-2u_{i,k}+u_{i+1,k}}{h^2} \\ &\frac{\partial{u(x,t)}}{\partial t}|_{i,k}=\frac{u_{i,k+1}-u_{i,k}}{\tau}\\ \end{aligned} \right. \end{equation}$

其中

和

$\tau$

分別為空間步長和時間步長，

和

分別標記空間指標和時間指標，則我們得到差分方程

$\begin{equation} \left\{ \begin{aligned} &\frac{u_{i,k+1}-u_{i,k}}{\tau}=\lambda \frac{u_{i-1,k}-2u_{i,k}+u_{i+1,k}}{h^2} \\ &u_{i,0}=f(ih) \qquad i=1, 2, ..., N-1, N=l/h \\ &u_{0,k}=g_1(k\tau)\qquad k=0, 1, ..., M, M=T/\tau \\ &u_{N,k}=g_2(k\tau) \end{aligned} \right. \end{equation}$

由此得到遞推關係

$\begin{equation} \left\{ \begin{aligned} &u_{i,k+1}=\frac{\lambda\tau}{h^2}u_{i+1,k}+(1-2\frac{\lambda \tau}{h^2})u_{i,k}+\frac{\lambda\tau}{h^2}u_{i-1,k} \\ &u_{i,0}=f(ih) \qquad i=1, 2, ..., N-1, N=l/h \\ &u_{0,k}=g_1(k\tau)\qquad k=0, 1, ..., M, M=T/\tau \\ &u_{N,k}=g_2(k\tau) \end{aligned} \right. \end{equation}$

下面我們考察一個具體的例子，一維熱傳導方程的求解。

3.1 一維熱傳導方程的求解

一維熱傳導方程是一個典型的拋物型二階偏微分方程。設

表示在時間

，空間

處的溫度，則根據傅立葉定律（單位時間內流經單位面積的熱量和該處溫度的負梯度成正比），可以匯出熱傳導方程【2】

$\begin{equation} \frac{\partial u}{\partial t}=\lambda \frac{\partial^2 u}{\partial x^2} \end{equation}$

其中

$\lambda \equiv \kappa/c\rho$

稱為熱擴散率，

$\kappa, c, \rho$

分別為熱導率，比熱和質量密度，都是由系統本身確定的常量。

為了具體，設

$\lambda=1, l=3, T=1$

，設邊界條件為

$\begin{equation} \left\{ \begin{aligned} &\frac{\partial u(x,t)}{\partial t}= \frac{\partial ^2 u(x,t) }{\partial x^2} \qquad 0\leq t\leq 1000, 0\leq x \leq 1 \\ &u(x,0)=4x(3-x)\\ &u(0,t)=0\\ &u(3,t)=0 \end{aligned} \right. \end{equation}$

設步長為：

$h=0.1, \tau=0.0001$

，從而

$\tau\lambda/h^2=1/100，N=l/h=30, M=T/\tau=10000$

，所以遞推關係為

$\begin{equation} \left\{ \begin{aligned} &u_{i,k+1}=\frac{1}{100}u_{i+1,k}+\frac{49}{50}u_{i,k}+\frac{1}{100}u_{i-1,k} \qquad i=1, 2, ...,29 \qquad k=0, 1, ...,9999 \\ &u_{i,0}=4ih(3-ih) \qquad i=0, 1, ...,29 \\ &u_{0,k}=0 \qquad k=0, 1, ...,100000 \\ &u_{N,k}=0 \end{aligned} \right. \end{equation}$

圖2：遞推關係示意圖

圖2直觀地給出了差分法求解偏微分方程的過程。先把時空座標都離散化，根據遞推關係，由下一行的三個藍點的值可以給出上一行的一個紅點的值，由於邊界條件和初始條件（即最下方和兩邊的綠線）已知，所以按這個遞推關係可以得到網格中的所有值。下面我們用Python程式碼來實現。

import numpy as np

import matplotlib。pyplot as plt

h = 0。1#空間步長

N =30#空間步數

dt = 0。0001#時間步長

M = 10000#時間的步數

A = dt/（h**2） #lambda*tau/h^2

U = zeros（［N+1，M+1］）#建立二維空陣列

Space = arange（0，（N+1）*h，h）#建立空間等差數列，從0到3，公差是h

#邊界條件

for k in arange（0，M+1）：

U［0，k］ = 0。0

U［N，k］ = 0。0

#初始條件

for i in arange（0，N）：

U［i，0］=4*i*h*（3-i*h）

#遞推關係

for k in arange（0，M）：

for i in arange（1，N）：

U［i，k+1］=A*U［i+1，k］+（1-2*A）*U［i，k］+A*U［i-1，k］

上述程式碼中，我們首先把位於0-3中的空間等分為30份，位於0-1的時間等分為10000份，然後透過設定初始條件、邊界條件和遞推關係並藉助for迴圈就得到了1個30*10000的二維陣列，裡面放著每個離散的時空點的溫度值。

為了直觀地展現溫度隨時空的變化關係，接下來開始畫圖，首先畫出不同時刻溫度隨空間座標的變化

#不同時刻的溫度隨空間座標的變化

plt。plot（Space，U［：，0］， ‘g-’， label=‘t=0’，linewidth=1。0）

plt。plot（Space，U［：，3000］， ‘b-’， label=‘t=3/10’，linewidth=1。0）

plt。plot（Space，U［：，6000］， ‘k-’， label=‘t=6/10’，linewidth=1。0）

plt。plot（Space，U［：，9000］， ‘r-’， label=‘t=9/10’，linewidth=1。0）

plt。plot（Space，U［：，10000］， ‘y-’， label=‘t=1’，linewidth=1。0）

plt。ylabel（‘u（x，t）’， fontsize=20）

plt。xlabel（‘x’， fontsize=20）

plt。xlim（0，3）

plt。ylim（-2，10）

plt。legend（loc=‘upper right’）

plt。show（）

圖3：不同時刻的溫度隨空間座標的變化

從圖3可以看到，溫度關於

呈現軸對稱分佈，這是由初始條件造成的。另外，對每一點的空間座標，隨著時間的推移，溫度越來越低。

接下來，我們來畫出溫度等高線來描述溫度隨任意時空點的變化

#溫度等高線隨時空座標的變化，溫度越高，顏色越偏紅

extent = ［0，1，0，3］#時間和空間的取值範圍

levels = arange（0，10，0。1）#溫度等高線的變化範圍0-10，變化間隔為0。1

plt。contourf（U，levels，origin=‘lower’，extent=extent，cmap=plt。cm。jet）

plt。ylabel（‘x’， fontsize=20）

plt。xlabel（‘t’， fontsize=20）

plt。show（）

圖4：溫度隨時空座標的變化，溫度越高，顏色越紅

在圖4中，利用顏色的深淺來標記溫度，溫度越高，顏色越紅。從中同樣可以看到，溫度隨空間的分佈關於

軸對稱，而且隨著時間的推移，溫度越來越低。

4 總結

在本文中，我們利用Python數值求解了常微分方程和偏微分方程，基本思想是先將連續的座標離散化，然後將微分化為差分，由差分方程得到遞推關係，然後利用計算機強大的重複計算能力得到任意格點處的函式值。雖然上面只算了兩個例子，但是這種方法完全可以推廣到任意偏微分方程的求解中。

在量子色動力學（QCD）中，由於強相互作用具有漸進自由的特性，所以在低能情況下沒辦法像QED那樣使用微擾論計算，這時就要採用格點QCD的方法計算。其基本思想也是將時空離散化，然後從第一性原理的路徑積分出發去計算。由於時空被離散化了，相當於人為地引入了一個最小時空距離

，在傅立葉變換到動量空間後相當於引入了一個最大的動量截斷

$\Lambda$

，所以計算結果不會出現紫外發散，從而可以算到很高的精度，在一些情況下，格點的計算結果甚至比實驗更精確。所以，將連續引數離散化，把微分化為差分的思想，是極其重要的。

【1】不同的演算法對於方程右邊具體取什麼形式並不一樣，從而精度也不一樣。例如，尤拉法右邊取得是

；改進的尤拉法右邊取的是

$(f(y_n,t_n)+f(y_{n+1},t_{n+1}))/2$

；二階Runge-Kutta法右邊取的是

$f(y_{n+1/2},t_{n+1/2})$

；四階Runge-Kutta法右邊取的是

，其中

，

$k_2=f(y_n+k_1\Delta t/2,t_{n+1/2})$

，

$k_3=f(y_n+k_2\Delta t/2,t_{n+1/2})$

，

$k_4=f(y_n+k_3\Delta t,t_{n+1})$

。這些演算法的差別在於計算精度不同，並不改變差分的本質思想。為了具體，我們這裡採用改進的尤拉法。

【2】傅立葉定律告訴我們單位時間透過單位面積的熱量和該處的溫度負梯度成正比，即

$J=-\kappa\nabla u$

，其中

是熱流，即單位時間透過單位面積的熱量，

$\kappa$

為熱導率，

為溫度。能量守恆定律告訴我們單位時間流出某閉合曲面的熱量等於其內部減少的能量，即

$\oint_S J \cdot dS=-\frac{\partial}{\partial t}\int_V Cu dV$

，其中

$C=\rho c$

表示單位體積的熱容，而

$\rho$

和

分別為質量密度和單位質量的熱容（比熱），聯合散度定理

$\oint_S J \cdot dS=\int_V \nabla\cdot J dV$

，我們得到

$\nabla\cdot J=-C\frac{\partial u}{\partial t}$

，再將傅立葉定律帶入，就得到了熱傳導方程

$\frac{\partial u}{\partial t}=\lambda \nabla^ 2u$

，其中

$\lambda \equiv \kappa/C=\kappa/c\rho$

稱為熱擴散率。在一維的情況下，熱傳導方程就退化到了正文中的形式，即

$\frac{\partial u}{\partial t}=\lambda\frac{\partial ^2u}{\partial x^2}$

。

用Python數值求解偏微分方程

HER2是個什麼鬼？為什麼說“她”是乳腺癌的“幫兇”？

哈利波特中佛地魔和鄧布利多誰厲害

隨便看看

小米手機發出嘯叫？

盔犀鳥目前在世界上還有多少隻？

s引數表示什麼意思？

羊油釣鯽魚絕招？

用Python數值求解偏微分方程

HER2是個什麼鬼？為什麼說“她”是乳腺癌的“幫兇”？

哈利波特中佛地魔和鄧布利多誰厲害

猜你喜歡

世界上著名的數列有哪些？

什麼叫和數列？

帶不等式約束凸最佳化問題的求解——內點法(七)

隨便看看

小米手機發出嘯叫？

盔犀鳥目前在世界上還有多少隻？

s引數表示什麼意思？

羊油釣鯽魚絕招？