# 量化分析師的Python日記【第10天 Q Quant兵器譜 -之偏微分方程1】 > 來源：https://uqer.io/community/share/5530d9f1f9f06c8f3390465a > 從今天開始我們將進入一個系列 —— 偏微分方程。作為這一系列的開篇，我們以熱傳導方差為引子，引出： > > 1. 如何提一個偏微分方程的初邊值問題； > 1. 利用差分格式將偏微分方程離散化； > 1. 顯示差分格式； > 1. 顯示差分格式的條件穩定性。 > > 最后一點將作為伏筆，引出我們下一天的學習：無條件穩定格式。 ## 1. 熱傳導方程  其中： + `κ` 稱為熱傳導系數 + `[2]` 稱為方程的初值條件（Initial Condition) + `[3][4]` 稱為方程的邊值條件 （Boundaries Condition）。這里我們使用Dirichlet條件 我們可以看一下初值條件的形狀: ```py from matplotlib import pylab import seaborn as sns import numpy as np font.set_size(20) def initialCondition(x): return 4.0*(1.0 - x) * x xArray = np.linspace(0,1.0,50) yArray = map(initialCondition, xArray) pylab.figure(figsize = (12,6)) pylab.plot(xArray, yArray) pylab.xlabel('$x$', fontsize = 15) pylab.ylabel('$f(x)$', fontsize = 15) pylab.title(u'一維熱傳導方程初值條件', fontproperties = font) <matplotlib.text.Text at 0x12523810> ```  ## 2. 顯式差分格式 這里的基本思想是用差分格式替換對應的微分形式，并且期盼兩種格式的"誤差"在網格足夠密的情況下會趨于0。我們分別在時間方向以及空間方向做差分格式：  合并在一起，我們就得到了原始微分方程的差分格式：  這里我們使用差分網格上的近似值`Uj,k`代替`uj,k`，得到新的方程：  到這里我們得到一個迭代方程組：  其中。下面我們使用Python代碼實現上面的過程。 首先定義基本變量： + `N` 空間方向的網格數 + `M` 時間方向的網格數 + `T` 最大時間期限 + `X` 最大空間范圍 + `U` 用來存儲差分網格點上值得矩陣 ```py N = 25 # x方向網格數 M = 2500 # t方向網格數 T = 1.0 X = 1.0 xArray = np.linspace(0,X,N+1) yArray = map(initialCondition, xArray) starValues = yArray U = np.zeros((N+1,M+1)) U[:,0] = starValues ``` ```py dx = X / N dt = T / M kappa = 1.0 rho = kappa * dt / dx / dx ``` 這里我們做正向迭代：迭代時 `k=0,1...M?1`, 代表我們從0時刻運行至`T` ```py for k in range(0, M): for j in range(1, N): U[j][k+1] = rho * U[j-1][k] + (1. - 2*rho) * U[j][k] + rho * U[j+1][k] U[0][k+1] = 0. U[N][k+1] = 0. ``` 我們可以畫出不同時間點 `U(,˙τk)` 的結果： ```py pylab.figure(figsize = (12,6)) pylab.plot(xArray, U[:,0]) pylab.plot(xArray, U[:, int(0.10/ dt)]) pylab.plot(xArray, U[:, int(0.20/ dt)]) pylab.plot(xArray, U[:, int(0.50/ dt)]) pylab.xlabel('$x$', fontsize = 15) pylab.ylabel(r'$U(\dot, \tau)$', fontsize = 15) pylab.title(u'一維熱傳導方程', fontproperties = font) pylab.legend([r'$\tau = 0.$', r'$\tau = 0.10$', r'$\tau = 0.20$', r'$\tau = 0.50$'], fontsize = 15) <matplotlib.legend.Legend at 0x12577cd0> ```  也可以通過三維立體圖看一下整體的熱傳導過程： ```py tArray = np.linspace(0, 0.2, int(0.2 / dt) + 1) xGrids, tGrids = np.meshgrid(xArray, tArray) ``` ```py from mpl_toolkits.mplot3d import Axes3D from matplotlib import cm fig= pylab.figure(figsize = (16,10)) ax = fig.add_subplot(1, 1, 1, projection = '3d') surface = ax.plot_surface(xGrids, tGrids, U[:,:int(0.2 / dt) + 1].T, cmap=cm.coolwarm) ax.set_xlabel("$x$", fontdict={"size":18}) ax.set_ylabel(r"$\tau$", fontdict={"size":18}) ax.set_zlabel(r"$U$", fontdict={"size":18}) ax.set_title(u"熱傳導方程 $u_\\tau = u_{xx}$" , fontproperties = font) fig.colorbar(surface,shrink=0.75) <matplotlib.colorbar.Colorbar instance at 0xf6eb878> ```  ## 3. 組裝起來 就像在前一天二叉樹建模中介紹的一樣，我們這里會以面向對象的方式重新封裝分散的代碼，方便復用。首先是方程的描述： ```py class HeatEquation: def __init__(self, kappa, X, T, initialConstion = lambda x:4.0*x*(1.0-x), boundaryConditionL = lambda x: 0, boundaryCondtionR = lambda x:0): self.kappa = kappa self.ic = initialConstion self.bcl = boundaryConditionL self.bcr = boundaryCondtionR self.X = X self.T = T ``` 下面的是顯式差分格式的描述： ```py class ExplicitEulerScheme: def __init__(self, M, N, equation): self.eq = equation self.dt = self.eq.T / M self.dx = self.eq.X / N self.U = np.zeros((N+1, M+1)) self.xArray = np.linspace(0,self.eq.X,N+1) self.U[:,0] = map(self.eq.ic, self.xArray) self.rho = self.eq.kappa * self.dt / self.dx / self.dx self.M = M self.N = N def roll_back(self): for k in range(0, self.M): for j in range(1, self.N): self.U[j][k+1] = self.rho * self.U[j-1][k] + (1. - 2*self.rho) * self.U[j][k] + self.rho * self.U[j+1][k] self.U[0][k+1] = self.eq.bcl(self.xArray[0]) self.U[N][k+1] = self.eq.bcr(self.xArray[-1]) def mesh_grids(self): tArray = np.linspace(0, self.eq.T, M+1) tGrids, xGrids = np.meshgrid(tArray, self.xArray) return tGrids, xGrids ``` 有了以上的部分，現在整個過程可以簡單的通過初始化和一行關于`roll_back`的調用完成： ```py ht = HeatEquation(1.,1.,1.) scheme = ExplicitEulerScheme(2500,25, ht) scheme.roll_back() ``` 我們可以獲取與之前相同的圖像： ```py tGrids, xGrids = scheme.mesh_grids() fig= pylab.figure(figsize = (16,10)) ax = fig.add_subplot(1, 1, 1, projection = '3d') cutoff = int(0.2 / scheme.dt) + 1 surface = ax.plot_surface(xGrids[:,:cutoff], tGrids[:,:cutoff], scheme.U[:,:cutoff], cmap=cm.coolwarm) ax.set_xlabel("$x$", fontdict={"size":18}) ax.set_ylabel(r"$\tau$", fontdict={"size":18}) ax.set_zlabel(r"$U$", fontdict={"size":18}) ax.set_title(u"熱傳導方程 $u_\\tau = u_{xx}$" , fontproperties = font) fig.colorbar(surface,shrink=0.75) <matplotlib.colorbar.Colorbar instance at 0x12d69e60> ```  ## 4. 什么時候顯式格式會失敗？ 顯式格式不能任意取時間和空間的網格點數，即`M`與`N`不能隨意取值。我們稱顯式格式為條件穩定。特別地，需要滿足所謂CFL條件（Courant–Friedrichs–Lewy）：  例如： + `M` = 2500 + `N` = 25 則：  + `M` = 1200 + `N` = 25 則：  下面的代碼計算在第二種情形下的網格點計算過程： ```py ht = HeatEquation(1.,1.,1.) scheme = ExplicitEulerScheme(1200,25, ht) scheme.roll_back() ``` 我們可以通過下圖看到，在CFL條件無法滿足的情況下，數值誤差累計的結果（特別注意后面的鋸齒）： ```py tGrids, xGrids = scheme.mesh_grids() fig= pylab.figure(figsize = (16,10)) ax = fig.add_subplot(1, 1, 1, projection = '3d') cutoff = int(0.2 / scheme.dt) + 1 surface = ax.plot_surface(xGrids[:,:cutoff], tGrids[:,:cutoff], scheme.U[:,:cutoff], cmap=cm.coolwarm) ax.set_xlabel("$x$", fontdict={"size":18}) ax.set_ylabel(r"$\tau$", fontdict={"size":18}) ax.set_zlabel(r"$U$", fontdict={"size":18}) ax.set_title(u"熱傳導方程 $u_\\tau = u_{xx}$, $\\rho = 0.521$" , fontproperties = font) fig.colorbar(surface,shrink=0.75) <matplotlib.colorbar.Colorbar instance at 0x10f51b48> ```  今天的日記到此為止，這個問題我們會在下一篇中進行討論，引出無條件穩定格式：隱式差分格式（Implicit）。