前言
微積分(其實是教basic fourier analysis)期末專題。
我研究的主題是熱傳導方程式 ，並使用python模擬熱傳導的模型。
介紹
傅立葉冷卻定律描述熱在介質中傳播的規律，可以用兩種形式表達:
微分形式，關注局部的熱傳導率
積分形式，關注流出/入部分介質表面的熱
我使用的熱傳導方程就是傅立葉冷卻定律 的微分形式，來描述區域內溫度的演化。
正文
有限一維模型-數學推導
考慮一個理想的長棍，若$t=0$時，長棍上各點的溫度可用$f(x)$表達。
若可用一函數$u(x,t)$表達長棍一點的溫度於特定時間，求$u(x,t)$。
可知熱傳導方程式為: $\frac{\partial u}{\partial t}=\alpha \nabla^2 u$
因此處模型只有一維，故式中laplacian可改為: $\frac{\partial u}{\partial t}=\alpha \frac{\partial^2 u}{\partial x^2}$
在物理意義上為溫度隨時間改變的速度正比於相鄰溫差之差。
使用Neumann boundary condition: $\frac{\partial u}{\partial x}(0,t)=\frac{\partial u}{\partial x}(l,t)=0$
當然也可以使用Dirichlet boundary condition: $u(0,t)=u(l,t)=0$，但這裡只展示Neumann條件下的推導而已。
用分離變數法，假設$u(x,t)=X(x)T(t)$。
從$\frac{\partial u}{\partial t}=\alpha \nabla^2u$可知，$X(x)T’(t)=\alpha X’’(x)T(t)$。
令$\lambda=\alpha T(t)/T’(t)= X(x)/X’’(x)$。
先解出$X(x)$的部分，$\lambda X’’(x)-X(x)=0$。
利用特稱方程解此ODE，$\lambda m^2-1=0\iff m=\pm\sqrt\frac1\lambda$，所以$X(x)$的通解為$X(x)=c_1e^{\frac x{\sqrt\lambda}}+c_2e^{-\frac x{\sqrt\lambda}}$
檢查邊界條件。
$X’(x)=\frac{c_1}{\sqrt\lambda}e^{\frac x{\sqrt\lambda}}-\frac{c_2}{\sqrt\lambda}e^{-\frac x{\sqrt\lambda}}$.
$X’(0)=\frac{c_1}{\sqrt\lambda}-\frac{c_2}{\sqrt\lambda}=0\Rightarrow c_1=c_2$
$X’(l)=\frac{c_1}{\sqrt\lambda}e^{\frac l{\sqrt\lambda}}-\frac{c_2}{\sqrt\lambda}e^{-\frac l{\sqrt\lambda}}=0=\frac{c_1}{\sqrt\lambda}(e^{\frac{l}{\sqrt\lambda}}-e^{-\frac{l}{\sqrt\lambda}})\Rightarrow e^{\frac{l}{\sqrt\lambda}}=e^{-\frac{l}{\sqrt\lambda}}$
分別討論$\lambda$的正負:
若$\lambda>0$，無解。
若$\lambda=0$，$X(x)=0, u(x,t)=0$，此為平凡解。
若$\lambda<0$, $e^{i\frac{l}{\sqrt{|\lambda|}}}=e^{-i\frac{l}{\sqrt{|\lambda|}}}$。
利用歐拉公式，上式化簡為$2isin(\frac l{\sqrt{|\lambda|}})=0\iff|\lambda|=\frac{l^2}{n^2\pi^2}$。
由此可知，$|\lambda|$必滿足$|\lambda|=\frac{l^2}{n^2\pi^2}$。
回頭解$T(t)$，$\alpha T(t)/T’(t)=\lambda\Rightarrow\alpha T(t)=\lambda T’(t)\Rightarrow\alpha dt=\frac{\lambda}{T}dT\Rightarrow\int\frac \alpha\lambda dt=\int\frac1{T}dT=\frac\alpha\lambda t=ln(T)\Rightarrow ce^{t\frac \alpha\lambda}=T(t)$
最後，利用上述的解得到$u(x,t)=X(x)T(t)=K_ne^{-\alpha t\frac{n^2\pi^2}{l^2}}cos(x\frac{n\pi}{l})$
但注意此處有自由變數$n$，而根據此PDE的線性性質(若有$n$個解$f_i$，則$\sum_{i=1}^n f_i$構成解空間)，所以通解要寫成
$u(x,t)=\sum^{\infty}_{n=0} X_n(x)T_n(t)=\sum^\infty_{n=0}K_ncos(x\frac{n\pi}{l})e^{-\alpha t\frac{n^2\pi^2}{l^2}}$
當$t=0$時，$u(x,0)=f(x)=\sum^{\infty}_{n=0} X_n(x)T_n(t)=\sum^\infty_{n=0}K_ncos(x\frac{n\pi}{l})$
若改成用Dirichlet boundary condition，解則變為$u(x,t)=\sum^{\infty}_{n=0} X_n(x)T_n(t)=\sum^\infty_{n=0}K_nsin(x\frac{n\pi}{l})e^{-\alpha t\frac{n^2\pi^2}{l^2}}$
正是因為這兩種解都牽涉到將函數轉成用三角函數表達，而促使傅立葉發展出傅立葉級數。
有限一維模型-程式模擬
我是用python的matplotlib的套件來畫圖，想看這套件的用法請參考我以前的文章。
定積分用黎曼和計算。
2D version:1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
eps = 1e-8
def sec(x):
return 1 / np.cos(x);
l = np.pi
maxn = 1000
K = [0] * maxn
# original function
def f(x):
#return np.piecewise(x, [np.logical_and(np.less_equal(x, 3), np.greater_equal(x,2)), np.logical_not(np.logical_and(np.less_equal(x, 3), np.greater_equal(x,2)))], [1, 0])
return np.cos(x)
#return np.cos(np.sin(x * np.pi / l))
#return 1/3 * (x**3) - l/2 * (x**2)
#return (1/5) * (x**5) - ((7 + l) / 4) * (x**4) + ((10 + 7 * l) / 3) * (x**3) - 5*l *(x**2)
#return np.sin(np.pi * x**2 / (2 * l**2))
# calculate coefficient cos series
def calKn_cos():
dl = l / maxn
for n in range(maxn):
for i in range(maxn):
K[n] += f(dl * i) * np.cos(np.pi * n * dl * i / l)
K[n] *= (2 * dl / l)
K[0] /= 2
# calculate coefficient sin series
def calKn_sin():
dl = l / maxn
for n in range(maxn):
for i in range(maxn):
K[n] += f(dl * i) * np.sin(np.pi * n * dl * i / l)
K[n] *= (2 * dl / l)
# calculate coefficient
def calKn(mode = "cos"):
if mode == "cos":
calKn_cos()
elif mode == "sin":
calKn_sin()
else:
print("error argument")
# solution function
def u(x, t, alpha = 1, maxn = 1000):
s = 0
for n in range(maxn):
s += K[n] * np.cos(x * n * np.pi / l) * np.exp(-alpha * t * (n * np.pi / l) ** 2)
return s
fig = plt.figure()
X = np.linspace(0.5, l - 0.5, 1024)
calKn()
plt.plot(X, u(X, 0, alpha = 0.5), label = 't=0')
plt.plot(X, u(X, 1, alpha = 0.5), label = 't=1')
plt.plot(X, u(X, 2, alpha = 0.5), label = 't=2')
plt.plot(X, u(X, 4, alpha = 0.5), label = 't=4')
plt.legend(loc = 0)
plt.show()
3D version:1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
eps = 1e-8
def sec(x):
return 1 / np.cos(x);
l = np.pi
maxn = 1000
K = [0] * maxn
# original function
def f(x):
#return np.piecewise(x, [np.logical_and(np.less_equal(x, 3), np.greater_equal(x,2)), np.logical_not(np.logical_and(np.less_equal(x, 3), np.greater_equal(x,2)))], [1, 0])
return np.cos(x)
#return np.cos(np.sin(x * np.pi / l))
#return 1/3 * (x**3) - l/2 * (x**2)
#return (1/5) * (x**5) - ((7 + l) / 4) * (x**4) + ((10 + 7 * l) / 3) * (x**3) - 5*l *(x**2)
#return np.sin(np.pi * x**2 / (2 * l**2))
def calKn_cos():
dl = l / maxn
for n in range(maxn):
for i in range(maxn):
K[n] += f(dl * i) * np.cos(np.pi * n * dl * i / l)
K[n] *= (2 * dl / l)
K[0] /= 2
def calKn_sin():
dl = l / maxn
for n in range(maxn):
for i in range(maxn):
K[n] += f(dl * i) * np.sin(np.pi * n * dl * i / l)
K[n] *= (2 * dl / l)
def calKn(mode = "cos"):
if mode == "cos":
calKn_cos()
elif mode == "sin":
calKn_sin()
else:
print("error argument")
def u(x, t, alpha = 1, maxn = 1000):
s = 0
for n in range(maxn):
s += K[n] * np.cos(x * n * np.pi / l) * np.exp(-alpha * t * (n * np.pi / l) ** 2)
return s
fig = plt.figure()
ax = plt.axes(projection = '3d')
X = np.linspace(0.5, l - 0.5, 1024)
T = np.linspace(0, 10, 512)
calKn()
X, T = np.meshgrid(X, T)
U = u(X, T, alpha = 0.5)
ax.set_xlabel('x')
ax.set_ylabel('t')
ax.set_zlabel('temperature')
ax.plot_surface(X, T, U, cmap = 'rainbow')
plt.show()
圖片:
$f(x) = cos(x), l = 2\pi, \alpha = 0.5$
$f(x) = cos(x), l = \pi, \alpha = 0.5$
$f(x) = cos(sin(x * \pi/l)), l = 5, \alpha = 0.5$
$f(x) = \frac15x^5 - \frac{7+l}4x^4+\frac{10+7l}3x^3 - 5lx^2, l=5, \alpha = 0.5$
$f(x) = \begin{cases}1&,|x-2.5|<0.5\\\\0&,\text{else}\end{cases}, l=5, \alpha = 0.5$
原函數為方波。可以很清楚的看見吉布斯現象。
無限長一維模型-數學推導
具邊界條件的模型，距離邊界越近的地方，誤差值越大，若改成使用”無限長”的長棍，則在有限區間內幾乎無誤差值。
因現在模型為無限長，故沒有邊界條件，但熱方程與初始條件仍一樣。
同樣，假設$u(x,t)=X(x)T(t)$，同樣也會得到$\lambda=\frac{X(x)}{X’’(x)}=\frac{\alpha T(t)}{T’(t)}$，也同樣要解$\lambda X’’(x)-X(x)=0$。
但與有限長的長棍不同處即在此。
由前一例可知$\lambda < 0$，$X(x)=c_1e^{i\frac{x}{\sqrt{|\lambda|}}}+c_2e^{-i\frac{x}{\sqrt{|\lambda|}}}$
因無邊界條件，所以需要用其他方法化簡$X(x)$。
把$X(x)$寫成$\frac{c_1+c_2}2e^{i\frac{x}{\sqrt{|\lambda|}}}+\frac{c_1+c_2}2e^{-i\frac{x}{\sqrt{|\lambda|}}}+\frac{c_1-c_2}2e^{i\frac{x}{\sqrt{|\lambda|}}}-\frac{c_1-c_2}2e^{-i\frac{x}{\sqrt{|\lambda|}}}$
得到$X(x)=(c_1+c_2)cos(\frac{x}{\sqrt{|\lambda|}})+i(c_1-c_2)sin(\frac{x}{\sqrt{|\lambda|}})=Acos(\frac{x}{\sqrt{|\lambda|}})+Bsin(\frac{x}{\sqrt{|\lambda|}})$
$T(t)$則與前例一樣，$T(t)=ce^{t\frac \alpha\lambda}$
得到特解$u_{\lambda}(x,t)=(A_{\lambda}cos(\frac{x}{\sqrt{|\lambda|}})+B_{\lambda}sin(\frac{x}{\sqrt{|\lambda|}}))e^{-\frac{t\alpha}{|\lambda|}}$。
為了讓他好看點，令$p=\frac1{\sqrt{|\lambda|}}$。
則特解變為: $u_p(x,t)=(A_{p}cos(px)+B_{p}sin(px))e^{-p^2t\alpha}$
依據這個PDE的線性性質，即乘上一常數($\Delta p$)後仍為解答，$\sum^\infty_{n=0}\Delta p u_{n\Delta p}(x,t)$為通解。
$\lim_{\Delta p\to 0}\sum^\infty_{n=0}\Delta p u_{n\Delta p}(x,t)=\int^\infty_0u_p(x,t)dp=\int^\infty_0(A_{p}cos(px)+B_{p}sin(px))e^{-p^2t\alpha}dp=u(x,t)$
當$t=0$時，$u(x,0)=f(x)=\int^\infty_0(A_{p}cos(px)+B_{p}sin(px))dp$，其中，$A_p=\frac1\pi\int^\infty_{-\infty}f(x)cos(px)dx$ and $B_p=\frac1\pi\int^\infty_{-\infty}f(x)sin(px)dx$
可以自己驗證看看，這個解仍符合原先的PDE。
上式可進一步變換成複數型傅立葉轉換，過程與一般的傅立葉積分變成複數型傅立葉轉換一樣，故在此省略，直接寫出公式:
$u(x,t)=\frac1{2\pi}\int^{\infty}_{-\infty}F(p)e^{ipx}e^{-t\alpha p^2}dp$ ，其中， $F(p)=\int^{\infty}_{-\infty}f(v)e^{-ipv}dv$
一般會認為到此結束，但實際上因為$e^{-t\alpha p^2}$這項的存在，導致可以”交換積分順序”。
為何如此?
根據Fubini’s theorem，只要 $\int^{\infty}_{-\infty}(\int^{\infty}_{-\infty}|(e^{ipx}e^{-t\alpha p^2})f(v)e^{-ipv}|dv)dp$ 與 $\int^{\infty}_{-\infty}(\int^{\infty}_{-\infty}|(e^{ipx}e^{-t\alpha p^2})f(v)e^{-ipv}|dp)dv$ 兩者都收斂，則
$\int^{\infty}_{-\infty}(\int^{\infty}_{-\infty}(e^{ipx}e^{-t\alpha p^2})f(v)e^{-ipv}dv)dp$ 與 $\int^{\infty}_{-\infty}(\int^{\infty}_{-\infty}(e^{ipx}e^{-t\alpha p^2})f(v)e^{-ipv}dp)dv$同樣收斂，且結果相同。
檢查 $\int^{\infty}_{-\infty}(\int^{\infty}_{-\infty}|e^{ipx}e^{-t\alpha p^2}f(v)e^{-ipv}|dv)dp$ 是否收斂:
$\int^{\infty}_{-\infty}(\int^{\infty}_{-\infty}|e^{ipx}e^{-t\alpha p^2}f(v)e^{-ipv}|dv)dp=\\\\\int^{\infty}_{-\infty}e^{-t\alpha p^2}(\int^{\infty}_{-\infty}|e^{ip(x-v)}f(v)|dv)dp=\\\\\int^{\infty}_{-\infty}e^{-t\alpha p^2}(\int^{\infty}_{-\infty}|f(v)cos(p(x-v))+if(v)sin(p(x-v))|dv)dp=\\\\\int^{\infty}_{-\infty}e^{-t\alpha p^2}(\int^{\infty}_{-\infty}\sqrt{f^2(v)cos^2(p(x-v))+f^2(v)sin^2(p(x-v))}dv)dp=\\\\\int^{\infty}_{-\infty}e^{-t\alpha p^2}(\int^{\infty}_{-\infty}|f(v)|dv)dp$
檢查 $\int^{\infty}_{-\infty}(\int^{\infty}_{-\infty}|e^{ipx}e^{-t\alpha p^2}f(v)e^{-ipv}|dp)dv$ 是否收斂的步驟與上面相似，故省略。
兩者收斂的條件相同，都是需要$f(x)$在$\Bbb R$上絕對收斂，而這也正是傅立葉轉換存在的條件之一。
注意一般使用傅立葉轉換不能交換積分順序，原因是因為 $\frac{1}{2\pi}\int^{\infty}_{-\infty}(\int^{\infty}_{-\infty}|f(v)e^{-i\omega v} e^{i\omega x}| dv) d\omega$ 並不收斂(實際上兩種皆不收斂)。
利用交換積分順序
$u(x,t)=\frac1{2\pi}\int^{\infty}_{-\infty}f(v)(\int^{\infty}_{-\infty}e^{ip(x-v)}e^{-t\alpha p^2}dp)dv$。
注意中間的積分項
$\int^{\infty}_{-\infty}e^{ip(x-v)}e^{-t\alpha p^2}dp=\int^\infty_{-\infty}e^{-t\alpha p^2}cos(p(x-v))+ie^{-t\alpha p^2}sin(p(x-v))\,dp$
其中虛部為奇函數，故可消去，變為$\int^\infty_{-\infty}e^{-t\alpha p^2}cos(p(x-v))dp$
這個不定積分無法用初等函數寫出，但其定積分有很多做法: 曲線積分、積分內取微分、級數展開。
這裡用”積分內取微分”來解:
簡化上面的積分為 $\int^\infty_0e^{-s^2}cos(bs)ds$，
$I(b)=\int^\infty_0e^{-s^2}cos(bs)ds$
$I’(b)=\int^\infty_0\frac{\partial}{\partial b}(e^{-s^2}cos(bs))ds=\int^\infty_0-se^{-s^2}sin(bs)ds=-\frac12 bI(b)$
解此ODE:
$I’(b)=-\frac12 bI(b)\Rightarrow I(b)=\frac{\sqrt\pi}2e^{-\frac{b^2}4}$
令上述積分的$t\alpha p^2=s^2$ and $b=\frac{x-v}{\sqrt{t\alpha}}$，所以，$\int^\infty_{-\infty}e^{-t\alpha p^2}cos(p(x-v))dp=\sqrt\frac{\pi}{t\alpha}e^{-\frac{(x-v)^2}{4t\alpha}}$
答案變為
$u(x,t)=\frac1{2\pi}\int^{\infty}_{-\infty}f(v)(\sqrt\frac{\pi}{t\alpha}e^{-\frac{(x-v)^2}{4t\alpha}})dv=\frac1{2\sqrt{\pi t \alpha}}\int^{\infty}_{-\infty}f(v)e^{-\frac{(x-v)^2}{4t\alpha}}dv$
因為$t$在分母，會令$t=0$時程式實作不方便，可以令$z=\frac{v-x}{2\sqrt{t\alpha}}$，則最終答案:
無限長一維模型-程式模擬
這裡牽涉到上下界無限的瑕積分，我寫了兩種方法:
第一是直接截斷，將上下界當作一個很大的值，再直接當定積分做；
第二是利用變數變換，考慮$\int_{\Bbb R}f(x)dx$，令$x=tan k$，則剛剛的積分可寫為$\int^{\pi/2}_{-\pi/2}f(tan k)\cdot sec^2(k)dk$，如此，上下界無限的積分就可變為一般的定積分，
(若要更精確，還可以將$[\pi/2-\Delta k, \pi/2]$以及$[-\pi/2,-\pi/2+\Delta k]$這兩個被捨去的區間一律當作$f(tan(\pm\pi/2\mp\Delta k))\cdot sec^2(\pm\pi/2\mp\Delta k)$計算)
(不過在我實際測試中，肉眼看不太出來兩種方法的差別，我猜測是因為收斂速度很快導致的。)
2D version:1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
eps = 1e-8
def sec(x):
return 1 / np.cos(x);
# original function
def f(x):
#return np.piecewise(x, [np.less_equal(abs(x), 1), np.logical_not(np.less_equal(abs(x), 1))], [1, 0])
return np.exp(-x**2)
#return np.sinc(x)
# tan(k)-subtitute version
def sub_u(x, t, alpha = 1, maxn = 10000):
dk = np.pi / maxn
coe = dk / np.sqrt(np.pi)
s = 0
for n in range(1, maxn):
s += f(x + 2 * np.tan(-np.pi / 2 + dk * n) * np.sqrt(t * alpha)) * np.exp(-np.tan(-np.pi / 2 + dk * n) ** 2) * (sec(-np.pi / 2 + dk * n) ** 2)
s += f(x + 2 * np.tan(-np.pi / 2 + dk) * np.sqrt(t * alpha)) * np.exp(-np.tan(-np.pi / 2 + dk) ** 2) * (sec(-np.pi / 2 + dk) ** 2)
s += f(x + 2 * np.tan(np.pi / 2 - dk) * np.sqrt(t * alpha)) * np.exp(-np.tan(np.pi / 2 - dk) ** 2) * (sec(np.pi / 2 - dk) ** 2)
s *= coe
return s
# non-subtitute version
def nonsub_u(x, t, alpha = 1, L = -100, R = 100, maxn = 10000):
dk = (R - L) / maxn
coe = dk / np.sqrt(np.pi)
s = 0
for n in range(1, maxn - 1):
s += f(x + 2 * (L + dk * n) * np.sqrt(t * alpha)) * np.exp(-(L + dk * n)**2)
s *= coe
return s
fig = plt.figure()
X = np.linspace(-10, 10, 1024)
plt.plot(X, sub_u(X, 0, alpha = 0.5), label = 't=0')
plt.plot(X, sub_u(X, 1, alpha = 0.5), label = 't=1')
plt.plot(X, sub_u(X, 2, alpha = 0.5), label = 't=2')
plt.plot(X, sub_u(X, 4, alpha = 0.5), label = 't=4')
plt.legend(loc = 0)
plt.show()
3D version:1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
eps = 1e-8
def sec(x):
return 1 / np.cos(x);
# original function
def f(x):
#return np.piecewise(x, [np.less_equal(abs(x), 1), np.logical_not(np.less_equal(abs(x), 1))], [1, 0])
#return np.exp(-x**2)
return np.sinc(x)
# tan(k)-subtitute version
def sub_u(x, t, alpha = 1, maxn = 10000):
dk = np.pi / maxn
coe = dk / np.sqrt(np.pi)
s = 0
for n in range(1, maxn):
s += f(x + 2 * np.tan(-np.pi / 2 + dk * n) * np.sqrt(t * alpha)) * np.exp(-np.tan(-np.pi / 2 + dk * n) ** 2) * (sec(-np.pi / 2 + dk * n) ** 2)
s += f(x + 2 * np.tan(-np.pi / 2 + dk) * np.sqrt(t * alpha)) * np.exp(-np.tan(-np.pi / 2 + dk) ** 2) * (sec(-np.pi / 2 + dk) ** 2)
s += f(x + 2 * np.tan(np.pi / 2 - dk) * np.sqrt(t * alpha)) * np.exp(-np.tan(np.pi / 2 - dk) ** 2) * (sec(np.pi / 2 - dk) ** 2)
s *= coe
return s
# non-subtitute version
def nonsub_u(x, t, alpha = 1, L = -100, R = 100, maxn = 10000):
dk = (R - L) / maxn
coe = dk / np.sqrt(np.pi)
s = 0
for n in range(1, maxn - 1):
s += f(x + 2 * (L + dk * n) * np.sqrt(t * alpha)) * np.exp(-(L + dk * n)**2)
s *= coe
return s
fig = plt.figure()
ax = plt.axes(projection = "3d")
X = np.linspace(-5, 5, 1024)
T = np.linspace(0, 10, 512)
X, T = np.meshgrid(X, T)
U = sub_u(X, T, alpha = 0.5)
ax.set_xlabel('x')
ax.set_ylabel('t')
ax.set_zlabel('temperature')
ax.plot_surface(X, T, U, cmap = "rainbow")
plt.show()
圖片:
$f(x)=e^{-x^2}, \alpha = 0.5$
直接截斷計算積分:
變數替換:
$f(x)=\begin{cases}1&,|x|<1\\\\0&,\text{else}\end{cases}, \alpha = 0.5$
直接截斷計算積分:
變數替換:
$f(x)=sinc(x), \alpha = 0.5$
直接截斷計算積分:
變數替換:
無限長一維模型-離散化 1
既然都有CFT了，那試試看DFT能不能也產生同樣的效果。
將原函數取樣成一個個數據點，接著就可以利用DFT以及非正常(?)IDFT來處理:
直接寫出答案: $u(n,t)=\frac1{L}\sum^{L-1}_{k=0}(\sum^{L-1}_{m=0}f[m]e^{-i\frac{2\pi}{L}km})e^{i\frac{2\pi}{L}nk}e^{-t\alpha k^2}$
其中$\hat{f}[k] = \sum^{L-1}_{m=0}f[m]e^{-i\frac{2\pi}{L}km}$就是DFT。
無限長一維模型-離散化 1-程式模擬:
離散型就沒畫成3D圖了，只有2D疊層圖。
1 | import numpy as np |
圖片:
$f(x)=e^{-x^2},\alpha = 0.5$
$f(x)=\begin{cases}1&,|x|<1\\\\0&,\text{else}\end{cases}, \alpha = 0.5$
$f(x)=sinc(x), \alpha = 0.5$
雖然看起來蠻正確的，但對比連續的模型誤差蠻大的，看下圖:
$f(x)=e^{-x^2},\alpha = 0.5$
離散模擬:
連續模擬:
我猜測是因為隨隨便便就把$e^{-t\alpha k^2}$離散化導致的。
無限長一維模型-離散化 2
仔細想想剛剛那樣不嚴謹的轉換有點糟糕，認真分析一下:
最開始的解為:
$u(x,t)=\frac1{2\pi}\int^{\infty}_{-\infty}(\int^{\infty}_{-\infty}f(v)e^{-ipv}dv)e^{ipx}e^{-t\alpha p^2}dp$
左思右想，這不就是$u(x, t)=\mathscr F^{-1}[\mathscr F[f(x)]\cdot e^{-t\alpha p^2}]$。
改寫一下，
先計算出
所以答案又可以寫成
先寫個連續的摺積試試看效果:
可以看見都沒有剛剛隨便離散導致的誤差問題。
那既然，連續的摺積可以了，那離散的呢?
先切樣本出來，假設我們是在$[l,r]$上切出$L$個樣本。
，其中
將摺積離散化後: $u(n,t)=\Delta k \sum_{m=0}^{2L-1}f[n-m]G_t[m]$
裡面的$\Delta k$是歸一化係數，一定要乘，否則數據會失準。
無限長一維模型-離散化 2-程式模擬
實作上有個小缺陷，當$t=0$時，會產生divided by zero exception所以我用$t=0.01$去逼近。
那如果從數學的角度看，當$\lim_{t\to0}$時，$\mathscr F^{-1}[e^{-t\alpha p^2}]=\delta(x)$
$(f * \delta)(x)=f(x)$，所以數學上是沒問題的。
很巧合地(應該是必然的結果)，$\lim_{t\to0}\frac1{\sqrt{\pi t\alpha}}e^{-\frac{x^2}{4t\alpha}}$，可以寫成另一種形式(令$a=2\sqrt{t\alpha}$): $\lim_{a\to0}\frac1{a\sqrt{\pi}}e^{-\frac{x^2}{a^2}}$，而這正是$\delta(x)$在分佈意義上的定義。
1 | import numpy as np |
圖片:
$f(x)=e^{-x^2}, l = -3, r = 3, L = 100, \alpha = 0.5$
$f(x)=\begin{cases}1&,|x|<1\\\\0&,\text{else}\end{cases}, l = -3, r = 3, L = 100, \alpha = 0.5$
$f(x)=sinc(x), l = -3, r = 3, L = 100, \alpha = 0.5$
顯然結果比隨便離散化好多了。
無限長一維模型-離散化-FFT ver.
那既然有摺積了，何不來點FFT加速?
，其中 
$u(n,t)=(f*G_t)[n]=\mathscr F^{-1}[\mathscr F[f]\cdot\mathscr F[G_t]][n]$
將FFT套用進DFT，就可以加速整體運算。
無限長一維模型-離散化-FFT ver.-程式模擬
這裡我參考正向的蝶形網路，先將數據點重新排列，再不斷的計算小部分的傅立葉轉換。
1 | import numpy as np |
圖片:
$f(x)=e^{-x^2}, l = -3, r = 3, L = 256, \alpha = 0.5$
$f(x)=\begin{cases}1&,|x|<1\\\\0&,\text{else}\end{cases}, l = -3, r = 3, L = 256, \alpha = 0.5$
$f(x)=sinc(x), l = -3, r = 3, L = 256, \alpha = 0.5$
source code
https://github.com/OEmiliatanO/heat_transfer_simulation
後記
這次的專題做的特別久，大概花了我3、4天的時間，但也收穫頗豐。
複習了以前做過的ODE解法，也學到了一點PDE解法、Leibniz integral rule、還有數值方法的技巧。
雖然推導過程不算太嚴謹，有一些步驟我沒有辦法判斷是否合法，例如從有限變成無限時，要乘上一常數，然後取極限，且這個常數是解的參數，從直覺上來講，這感覺就非法的，但結果是對的。
不確定其他類似的PDE也可以這樣做，希望有人可以指點迷津。
最後的離散化應該是我苦惱最久的，前面的連續版本只花一個晚上就推導出來了，之後就開始寫code、debug，寫完之後才處理離散化的問題，接著想了一整天才想到怎麼好好離散化。
blog的上一篇跟這篇隔了整整一年呢。
2023.5.17 更: 網站的render好像有改版，有一些數學公式都壞了，所以換成圖片。