Python求解偏微分方程

Python求解偏微分方程

更新於 發佈於 閱讀時間約 4 分鐘

1 引言

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


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


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

2 常微分方程的差分求解

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

raw-image

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

raw-image


從而我們得到遞推關係

raw-image

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


2.1 RC迴路放電問題

對於一個 RC 迴路,我們有

raw-image

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

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

raw-image


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

raw-image

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

raw-image


下面我們用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()
raw-image

在COLAB跑了一下

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

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

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

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

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


avatar-img
于正龍(Ricky)的沙龍
36會員
53內容數
人工智能工作經驗跟研究
留言
avatar-img
留言分享你的想法!
根據泰勒定理,f(x)可以寫成右邊一長串的導數的組合 為了更好理解這個東西我們可以用python實作 首先定義f(x)和定義factorial怎麼算 然後寫泰勒定理 f(x) = f(a) + f'(a)(x-a) ....後面還有一串 注意公式往後面看其實是有規律的 例如從
一個罐頭其中高度為h,底部半徑為r,且 h/r=2 已知 容量V = 1(公升) 找到一組: h(高度) , r(底部半徑) 使得表面積最小 目的: 因表面積最小因此可以讓製造成本降低幫助企業省錢 算式如上 答案我也使用python驗證一下求出數值解 接下來換個角度跟工具
根據泰勒定理,f(x)可以寫成右邊一長串的導數的組合 為了更好理解這個東西我們可以用python實作 首先定義f(x)和定義factorial怎麼算 然後寫泰勒定理 f(x) = f(a) + f'(a)(x-a) ....後面還有一串 注意公式往後面看其實是有規律的 例如從
一個罐頭其中高度為h,底部半徑為r,且 h/r=2 已知 容量V = 1(公升) 找到一組: h(高度) , r(底部半徑) 使得表面積最小 目的: 因表面積最小因此可以讓製造成本降低幫助企業省錢 算式如上 答案我也使用python驗證一下求出數值解 接下來換個角度跟工具