目錄
- 拋物方程的有限差分方法
- 馮諾依曼分析法
- 追趕法
- Crank-Nicolson方法
- 邊界條件
- 馮諾依曼分析法
拋物型方程的有限差分方法
最為常見(jiàn)的拋物型偏微分方程就是熱傳導(dǎo)方程
它本質(zhì)上其實(shí)是一個(gè)初邊值問(wèn)題翁逞，直觀來(lái)看是因?yàn)樘悠牵刃枰踔禇l件讹俊，又需要邊值條件。
對(duì)于偏導(dǎo)數(shù)的差分格式，一般分為隱式差分與顯式差分兩種，下面分別進(jìn)行介紹。
顯式差分格式
對(duì)于u的一階導(dǎo)數(shù)牵囤，可以進(jìn)行如下的近似：
其中，為設(shè)定的空間步長(zhǎng)滞伟。如果設(shè)定h為單位長(zhǎng)度揭鳞，則
其中，是后向差分的符號(hào)诗良，而不是梯度運(yùn)算汹桦。
由此進(jìn)而可以推導(dǎo)出的表達(dá)式為：
然后，再將代回鉴裹，便可得到最終的差分格式：
該式子中均是對(duì)空間的差分舞骆，所以并沒(méi)有變化钥弯。
在對(duì)于時(shí)間的差分中，我們便有兩種不同的選擇督禽，這也對(duì)應(yīng)了顯式和隱式差分格式脆霎。
對(duì)于顯式差分，其時(shí)間的差分方式如下：
其中狈惫，為設(shè)定的時(shí)間步長(zhǎng)睛蛛。可以引入變量
我們定義就是在
處的數(shù)值解胧谈，那么自然的忆肾，我們可以將它按照時(shí)間的順序排列，得到下面的一個(gè)形式:
馮諾依曼分析法
它可以應(yīng)用在各種常見(jiàn)的偏微分方程中菱肖。所謂的馮諾依曼分析法客冈，其本質(zhì)就是用向量，將所有的空間離散點(diǎn)拼起來(lái)稳强。那么這樣的話场仲，這個(gè)矩陣所對(duì)應(yīng)的表達(dá)式就僅僅與時(shí)間有關(guān)了。這樣自然分析就會(huì)方便很多退疫。
上述的差分形式渠缕，可以寫(xiě)成矩陣形式，如下：
寫(xiě)的緊湊一點(diǎn)褒繁，就是公式,那么為了讓這個(gè)公式收斂亦鳞，我們顯然需要讓
的最大特征值小于1，最后一般需要
實(shí)際的求解中棒坏，可以直接進(jìn)行遞推蚜迅，然后對(duì)于首尾直接代入邊界條件即可。
#一維熱傳導(dǎo)
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.animation as animation
a = 1
#第三類(lèi)邊界條件
c = 1
#dx太小或dt太大可能出錯(cuò)
dx = 0.01
dt = 0.00005
x = np.arange(0, 1, dx)
t = np.arange(0, 1, dt)
u = np.zeros((len(x), len(t)))
#初始條件
u[:, 0] = 0
#邊界條件
m1 = lambda t : 1 + 0.0 * np.sin(t)
m2 = lambda t : 2 - 0.0 * np.sin(10 * t)
#構(gòu)造對(duì)角矩陣1俊抵，-2，1
#其中diag為對(duì)角矩陣坐梯，1表示上移一格
A = -2 * np.eye((len(x))) + np.diag([1] * (len(x) - 1), 1) + np.diag([1] * (len(x) - 1), -1)
for n in range(len(t) - 1):
u[:, n + 1] = u[:, n] + ((a**2 /(dx**2)) * np.dot(A, u[:, n]) + f.T) * dt
#第一類(lèi)邊界條件
#用邊界條件確定第一個(gè)和最后一個(gè)
u[0, n + 1] = m1(n + 1)
u[-1, n + 1] = m2(n + 1)
figure = plt.figure()
ax = Axes3D(figure)
T, X = np.meshgrid(t, x)
ax.plot_surface(X, T, u, cmap='rainbow')
plt.xlabel('x')
plt.ylabel('t')
# plt.plot(x, u[:, -1])
# plt.grid()
plt.show()
若得到的解出現(xiàn)震蕩徽诲，則可以適當(dāng)增加dx或減小dt，最后得到的求解結(jié)果如下圖：
隱式差分格式
與顯式差分格式相比所不同的是對(duì)時(shí)間的差分格式吵血，其差分格式如下
那么也容易通過(guò)化簡(jiǎn)得到我們的矩陣方程:
寫(xiě)成緊湊形式便是谎替，，我們可以通過(guò)計(jì)算得到A的特征值為
蹋辅，所以與顯式差分相比該格式是無(wú)條件穩(wěn)定的钱贯。但是并不能直接按時(shí)間進(jìn)行遞推求解，需要不斷求解矩陣方程侦另，而
是以三對(duì)角矩陣秩命，較為稀疏尉共，直接對(duì)其進(jìn)行求逆會(huì)造成效率低下，于是需要利用追趕法進(jìn)行求解弃锐，可以大幅提升求解效率袄友。
追趕法
對(duì)于一矩陣方程：
其中為一非奇異三對(duì)角矩陣，可以對(duì)
進(jìn)行LU分解霹菊，得到一上三角矩陣和一下三角矩陣剧蚣。之后可引入中間變量
,則有
之后已知便可求得
,已知
后便可求得
。具體代碼實(shí)現(xiàn)如下：
import numpy as np
from scipy import linalg
#可解決大規(guī)模的三對(duì)角的矩陣方程
#如在偏微分方程的求解中
'''
b c 0 0
a b c 0
0 a b c
0 0 a b
其中b的絕對(duì)值
大于等于a與c的絕對(duì)值之和
'''
A = -2 * np.eye((5)) + np.diag([1] * (5 - 1), 1) + np.diag([1] * (5 - 1), -1)
b = np.array([1,1,1,1,1])
# print(b[1])
def Thomas(La, Mb, Uc, b):
"""
La -- [lower item for tri-diagonal matrix] -- [a, a, a]
Mb -- [mian item for tri-diagonal matrix] -- [b, b, b, b]
Uc -- [upper item for tri-diagonal matrix] -- [c, c, c]
b -- [AX = b, where A is the tri-diagonal matrix] [1,1,1,1,1]
"""
n = len(Mb)
Uc[0] = Uc[0] / Mb[0]
for i in range(1, n-1):
Uc[i] = Uc[i] / (Mb[i] - La[i - 1] * Uc[i - 1])
b[0] = b[0] / Mb[0]
for i in range(1, n):
b[i] = (b[i] - La[i-1]*b[i-1]) / (Mb[i] - La[i-1] * Uc[i-1])
ls = list(range(n-1))[::-1]
for i in ls:
b[i] = b[i] - Uc[i]*b[i+1]
return b
之后便可利用追趕法和初始條件旋廷，遍歷時(shí)間節(jié)點(diǎn)鸠按，逐次求解每個(gè)時(shí)間點(diǎn)的溫度。
Crank-Nicolson方法
因?yàn)橹暗膬蓚€(gè)方法它們的誤差是,這其實(shí)會(huì)有問(wèn)題就是主要的誤差其實(shí)都在時(shí)間步長(zhǎng)上饶碘。我們需要把時(shí)間步長(zhǎng) 取得很小目尖，可能才能夠達(dá)到我們要的精度。因此
方法可以很好地彌補(bǔ)這個(gè)缺陷熊镣。
這里我們的使用向后差分卑雁，和上面那一個(gè)一樣。而我們的
使用的是混合差分绪囱。具體來(lái)說(shuō)测蹲，就是用:
最后按時(shí)間進(jìn)行整合后得到：
大體上的形式就是
其中：
可以計(jì)算得到的特征值為：
其特征值也是小于1的，故也是無(wú)條件穩(wěn)定鬼吵。同時(shí)其截?cái)嗾`差為,這比上面兩種方法的誤差要小很多了
