Python求解偏微分方程

1 引言

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

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

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

2 常微分方程的差分求解

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

上述微分方程可以化為差分方程

從而我們得到遞推關係

有了遞推關係和初始條件以後，就可以利用Python的強大運算功能來求解了

2.1 RC迴路放電問題

對於一個 RC 迴路，我們有

其中，I R Q C 分別為電流，電阻，電量和電容，利用 I=dQ/dt並定義

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

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

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

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

import numpy as np
import matplotlib.pyplot as plt
rc = 2.0 #設置常數
dt = 0.5 #設置步長
n = 1000 #設置分割段數
t = 0.0 #設置初始時間
q = 1.0 #設置初始電量

qt=[] #用來存放差分得到的q值
qt0=[] #用來存放解析得到的q值
time = [] #用來存放時間值

for i in range(n):
    t = t + dt
    q1 = q - q*dt/rc #qn+1的近似值
    q = q - 0.5*(q1*dt/rc + q*dt/rc) #差分遞推關係
    q0 = np.exp(-t/rc) #解析關係
    qt.append(q) #差分得到的q值列表
    qt0.append(q0) #解析得到的q值列表
    time.append(t) #時間列表

plt.plot(time,qt,'o',label='Euler-Modify') #差分得到的電量隨時間的變化
plt.plot(time,qt0,'r-',label='Analytical') #解析得到的電量隨時間的變化
plt.xlabel('time')
plt.ylabel('charge')
plt.xlim(0,20)
plt.ylim(-0.2,1.0)
plt.legend(loc='upper right')
plt.show()

在COLAB跑了一下

差分法得到的結果與解析法得到結果的比較，發現兩者符合得很好

這說明對於這個問題，改進的歐拉法已經可以給出足夠精確的結果。

要注意的是，這個微分方程本身比較簡單，可以解析求解，

而對於複雜得多的微分方程，沒辦法解析求解，

但是上述數值求解差分法仍然是適用的。