# 量化分析師的Python日記【第11天 Q Quant兵器譜之偏微分方程2】 > 來源：https://uqer.io/community/share/5534ad3ff9f06c8f33904689 > 這是量化分析師的偏微分方程系列的第二篇，在這一篇中我們將解決上一篇顯式格式留下的穩定性問題。本篇將引入隱式差分算法，讀者可以學到： > > 1. 隱式差分格式描述 > 1. 三對角矩陣求解 > 1. 如何使用`scipy`加速算法實現 > > 在完成兩天的基礎學習之后，在下一天中，我們將把已經學到的知識運用到金融定價領域最重要的方程之一：Black - Shcoles - Merton 偏微分方差 ```py from matplotlib import pylab import seaborn as sns import numpy as np np.set_printoptions(precision = 4) font.set_size(20) def initialCondition(x): return 4.0*(1.0 - x) * x ``` ## 1. 隱式差分格式 像上一天一樣，我們從差分格式的數學表述開始。隱式格式與顯式格式的區別，在于我們時間方向選擇的基準點。顯式格式使用`k`，而隱式格式選擇`k+1`:  剩下的推到過程我完全一樣，我們看到無論隱式格式還是顯式格式，它們的截斷誤差是一樣的：  用離散值`Uj,k`替換`uj,k`，我們得到差分方程：  最后，到這里我們得到一個迭代方程組：  其中。 ```py N = 500 # x方向網格數 M = 500 # 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 ``` ### 1.1 矩陣求解(`TridiagonalSystem`) 雖然看上去形式只是變了一點，但是求解的問題有很大的變化。在每個時間點上，我們需要求解如下的一個線性方程組：  這里` A`為：  幸運的是，這個是個三對角矩陣，可以很簡單的利用Gauss消去法求解。我們這里不會詳細討論算法的描述，細節都可以在下面的python類`TridiagonalSystem`中了解到： ```py class TridiagonalSystem: def __init__(self, udiag, cdiag, ldiag): ''' 三對角矩陣： udiag -- 上對角線 cdiag -- 對角線 ldiag -- 下對角線 ''' assert len(udiag) == len(cdiag) assert len(cdiag) == len(ldiag) self.udiag = udiag self.cdiag = cdiag self.ldiag = ldiag self.length = len(self.cdiag) def solve(self, rhs): ''' 求解以下方程組 A \ dot x = rhs ''' assert len(rhs) == len(self.cdiag) udiag = self.udiag.copy() cdiag = self.cdiag.copy() ldiag = self.ldiag.copy() b = rhs.copy() # 消去下對角元 for i in range(1, self.length): cdiag[i] -= udiag[i-1] * ldiag[i] / cdiag[i-1] b[i] -= b[i-1] * ldiag[i] / cdiag[i-1] # 從最后一個方程開始求解 x = np.zeros(self.length) x[self.length-1] = b[self.length - 1] / cdiag[self.length - 1] for i in range(self.length - 2, -1, -1): x[i] = (b[i] - udiag[i]*x[i+1]) / cdiag[i] return x def multiply(self, x): ''' 矩陣乘法： rhs = A \dot x ''' assert len(x) == len(self.cdiag) rhs = np.zeros(self.length) rhs[0] = x[0] * self.cdiag[0] + x[1] * self.udiag[0] for i in range(1, self.length - 1): rhs[i] = x[i-1] * self.ldiag[i] + x[i] * self.cdiag[i] + x[i+1] * self.udiag[i] rhs[self.length - 1] = x[self.length - 2] * self.ldiag[self.length - 1] + x[self.length - 1] * self.cdiag[self.length - 1] return rhs ``` ### 1.2 隱式格式求解 ```py for k in range(0, M): udiag = - np.ones(N-1) * rho ldiag = - np.ones(N-1) * rho cdiag = np.ones(N-1) * (1.0 + 2. * rho) mat = TridiagonalSystem(udiag, cdiag, ldiag) rhs = U[1:N,k] x = mat.solve(rhs) U[1:N, k+1] = x U[0][k+1] = 0. U[N][k+1] = 0. ``` ```py from lib.utilities import plotLines plotLines([U[:,0], U[:, int(0.10/ dt)], U[:, int(0.20/ dt)], U[:, int(0.50/ dt)]], xArray, title = u'一維熱傳導方程', xlabel = '$x$', ylabel = r'$U(\dot, \tau)$', legend = [r'$\tau = 0.$', r'$\tau = 0.10$', r'$\tau = 0.20$', r'$\tau = 0.50$']) ```  ```py from lib.utilities import plotSurface tArray = np.linspace(0, 0.2, int(0.2 / dt) + 1) tGrids, xGrids = np.meshgrid(tArray, xArray) plotSurface(xGrids, tGrids, U[:,:int(0.2 / dt) + 1], title = u"熱傳導方程 $u_\\tau = u_{xx}$，隱式格式（$\\rho = 50$）", xlabel = "$x$", ylabel = r"$\tau$", zlabel = r"$U$") ```  ## 2. 繼續組裝 像我們在顯示格式那一節介紹的同樣做法，我們把之前的代碼整合起來，歸集與一個完整的類`ImplicitEulerScheme`中： ```py from lib.utilities import HeatEquation ``` 上面的代碼（使用`library`功能，關于該功能的具體介紹請見[幫助 — Library是干什么的](https://app.wmcloud.com/mercury/help/faq/#Library是干什么的))導入我們在上一期中已經定義過的類`HeatEquation`，避免代碼重復。 ```py class ImplicitEulerScheme: 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): udiag = - np.ones(self.N-1) * self.rho ldiag = - np.ones(self.N-1) * self.rho cdiag = np.ones(self.N-1) * (1.0 + 2. * self.rho) mat = TridiagonalSystem(udiag, cdiag, ldiag) rhs = self.U[1:self.N,k] x = mat.solve(rhs) self.U[1:self.N, k+1] = x self.U[0][k+1] = self.eq.bcl(self.xArray[0]) self.U[self.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 ``` 然后我們可以使用下面的三行簡單調用完成功能： ```py ht = HeatEquation(1.,X, T) scheme = ImplicitEulerScheme(M,N, ht) scheme.roll_back() scheme.U array([[ 0.0000e+00, 0.0000e+00, 0.0000e+00, ..., 0.0000e+00, 0.0000e+00, 0.0000e+00], [ 7.9840e-03, 7.2843e-03, 6.9266e-03, ..., 3.8398e-07, 3.7655e-07, 3.6926e-07], [ 1.5936e-02, 1.4567e-02, 1.3852e-02, ..., 7.6795e-07, 7.5308e-07, 7.3851e-07], ..., [ 1.5936e-02, 1.4567e-02, 1.3852e-02, ..., 7.6795e-07, 7.5308e-07, 7.3851e-07], [ 7.9840e-03, 7.2843e-03, 6.9266e-03, ..., 3.8398e-07, 3.7655e-07, 3.6926e-07], [ 0.0000e+00, 0.0000e+00, 0.0000e+00, ..., 0.0000e+00, 0.0000e+00, 0.0000e+00]]) ``` ## 3. 使用 `scipy`加速 軟件工程行業里有句老話，叫做：“不要重復發明輪子！”。實際上，之前的代碼里面，我們就造了自己的輪子：`TridiagonalSystem`。三對角矩陣作為最最常見的稀疏矩陣，關于它的線性方程組求解算法實際上早已為業界熟知，也已經有很多庫內置了工業級別強度實現。這里我們取`scipy`作為例子，來展示使用外源庫實現的好處： + 更加穩健的算法： 知名庫算法由于使用者廣泛，有更大的概率發現一些極端情形下的bug。庫作者可以根據用戶反饋，及時調整算法； + 更高的性能： 由于庫的使用更為廣泛，庫作者有更大的動力去使用各種技術去提高算法的性能：例如使用更高效的語言實現，例如C。scipy中的情形就是一例。 + 持續的維護： 庫的受眾范圍廣，社區的力量會推動庫作者持續維護。 下面的代碼展示，如何使用`scipy`中的`solve_banded`算法求解三對角矩陣： ```py import scipy as sp from scipy.linalg import solve_banded A = np.zeros((3, 5)) A[0, :] = np.ones(5) * 1. # 上對角線 A[1, :] = np.ones(5) * 3. # 對角線 A[2, :] = np.ones(5) * (-1.) # 下對角線 b = [1.,2.,3.,4.,5.] x = solve_banded ((1,1), A,b) print 'x = A^-1b = ',x x = A^-1b = [ 0.1833 0.45 0.8333 0.95 1.9833] ``` 我們使用上面的算法替代我們之前的`TridiagonalSystem`， ```py import scipy as sp from scipy.linalg import solve_banded for k in range(0, M): udiag = - np.ones(N-1) * rho ldiag = - np.ones(N-1) * rho cdiag = np.ones(N-1) * (1.0 + 2. * rho) mat = np.zeros((3,N-1)) mat[0,:] = udiag mat[1,:] = cdiag mat[2,:] = ldiag rhs = U[1:N,k] x = solve_banded ((1,1), mat,rhs) U[1:N, k+1] = x U[0][k+1] = 0. U[N][k+1] = 0. ``` ```py plotLines([U[:,0], U[:, int(0.10/ dt)], U[:, int(0.20/ dt)], U[:, int(0.50/ dt)]], xArray, title = u'一維熱傳導方程，使用scipy', xlabel = '$x$', ylabel = r'$U(\dot, \tau)$', legend = [r'$\tau = 0.$', r'$\tau = 0.10$', r'$\tau = 0.20$', r'$\tau = 0.50$']) ```  同樣的我們定義一個新類`ImplicitEulerSchemeWithScipy`使用`scipy`的算法： ```py class ImplicitEulerSchemeWithScipy: 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): udiag = - np.ones(self.N-1) * self.rho ldiag = - np.ones(self.N-1) * self.rho cdiag = np.ones(self.N-1) * (1.0 + 2. * self.rho) mat = np.zeros((3,self.N-1)) mat[0,:] = udiag mat[1,:] = cdiag mat[2,:] = ldiag rhs = self.U[1:self.N,k] x = solve_banded((1,1), mat, rhs) self.U[1:self.N, k+1] = x self.U[0][k+1] = self.eq.bcl(self.xArray[0]) self.U[self.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 ``` 下面的代碼，比較了兩種做法的性能。可以看到僅僅簡單的替代三對角矩陣算法，我們就獲得了接近8倍的性能提升 ```py import time startTime = time.time() loop_round = 10 # 不使用scipy for k in range(loop_round): ht = HeatEquation(1.,X, T) scheme = ImplicitEulerScheme(M,N, ht) scheme.roll_back() endTime = time.time() print '{0:<40}{1:.4f}'.format('執行時間(s) -- 不使用scipy.linalg: ', endTime - startTime) # 使用scipy startTime = time.time() for k in range(loop_round): ht = HeatEquation(1.,X, T) scheme = ImplicitEulerSchemeWithScipy(M,N, ht) scheme.roll_back() endTime = time.time() print '{0:<40}{1:.4f}'.format('執行時間(s) -- 使用scipy.linalg: ', endTime - startTime) 執行時間(s) -- 不使用scipy.linalg: 12.1589 執行時間(s) -- 使用scipy.linalg: 1.6224 ``` ## 4. 尾聲 到這里為止，我們已經結束了偏微分方差差分格式的基礎學習。這是一個很大的學科，這兩天也只能做到“管中窺豹”。但是有了以上的基礎知識，讀者已經有了足夠的積累，可以處理一些金融工程中會實際遇到的方程。在下一天中，我們將把這兩天學習到的知識運用到金融工程史上最重要的方程：Black - Scholes - Merton 偏微分方程。