量化分析师的 Python 日记【第 10 天 Q Quant 兵器谱 -之偏微分方程 1】

从今天开始我们将进入一个系列 —— 偏微分方程。作为这一系列的开篇,我们以热传导方差为引子,引出:

*如何提一个偏微分方程的初边值问题;

*利用差分格式将偏微分方程离散化;

*显示差分格式;

*显示差分格式的条件稳定性

最后一点将作为伏笔,引出我们下一天的学习:无条件稳定格式。

量化分析师的 Python 日记【第 10 天 Q Quant 兵器谱 -之偏微分方程 1】

热传导方程

我们这里使用 1D 热传导方程作为例子:

量化分析师的 Python 日记【第 10 天 Q Quant 兵器谱 -之偏微分方程 1】

其中:

*𝜅 称为热传导系数

*[2] 称为方程的初值条件(Initial Condition)

*[3][4] 称为方程的边值条件(Boundaries Condition)。这里我们使用Dirichlet 条件

我们可以看一下初值条件的形状 :

量化分析师的 Python 日记【第 10 天 Q Quant 兵器谱 -之偏微分方程 1】

量化分析师的 Python 日记【第 10 天 Q Quant 兵器谱 -之偏微分方程 1】

显式差分格式

这里的基本思想是用差分格式替换对应的微分形式,并且期盼两种格式的 " 误差 " 在网格足够密的情况下会趋于 0。我们分别在时间方向以及空间方向做差分格式:

量化分析师的 Python 日记【第 10 天 Q Quant 兵器谱 -之偏微分方程 1】

合并在一起,我们就得到了原始微分方程的差分格式:

量化分析师的 Python 日记【第 10 天 Q Quant 兵器谱 -之偏微分方程 1】

这里我们使用差分网格上的近似值 $U_{j,k}$ 代替 $u_{j,k}$,得到新的方程:

量化分析师的 Python 日记【第 10 天 Q Quant 兵器谱 -之偏微分方程 1】

到这里我们得到一个迭代方程组:

量化分析师的 Python 日记【第 10 天 Q Quant 兵器谱 -之偏微分方程 1】

其中 𝜌下面我们使用 Python 代码实现上面的过程。

首先定义基本变量:

*N 空间方向的网格数

*M 时间方向的网格数

*T 最大时间期限

*X 最大空间范围

*U 用来存储差分网格点上值得矩阵

量化分析师的 Python 日记【第 10 天 Q Quant 兵器谱 -之偏微分方程 1】

量化分析师的 Python 日记【第 10 天 Q Quant 兵器谱 -之偏微分方程 1】

这里我们做正向迭代:迭代时 k = 0, 1 ... M-1, 代表我们从 0 时刻运行至 `T```

量化分析师的 Python 日记【第 10 天 Q Quant 兵器谱 -之偏微分方程 1】

我们可以画出不同时间点的结果:

量化分析师的 Python 日记【第 10 天 Q Quant 兵器谱 -之偏微分方程 1】

也可以通过三维立体图看一下整体的热传导过程:

量化分析师的 Python 日记【第 10 天 Q Quant 兵器谱 -之偏微分方程 1】

量化分析师的 Python 日记【第 10 天 Q Quant 兵器谱 -之偏微分方程 1】

量化分析师的 Python 日记【第 10 天 Q Quant 兵器谱 -之偏微分方程 1】

量化分析师的 Python 日记【第 10 天 Q Quant 兵器谱 -之偏微分方程 1】

组装起来

就像在前一天二叉树建模中介绍的一样,我们这里会以面向对象的方式重新封装分散的代码,方便复用。首先是方程的描述:

量化分析师的 Python 日记【第 10 天 Q Quant 兵器谱 -之偏微分方程 1】

下面的是显式差分格式的描述:

量化分析师的 Python 日记【第 10 天 Q Quant 兵器谱 -之偏微分方程 1】

有了以上的部分,现在整个过程可以简单的通过初始化和一行关于 roll_back 的调用完成:

量化分析师的 Python 日记【第 10 天 Q Quant 兵器谱 -之偏微分方程 1】

我们可以获取与之前相同的图像:

量化分析师的 Python 日记【第 10 天 Q Quant 兵器谱 -之偏微分方程 1】

量化分析师的 Python 日记【第 10 天 Q Quant 兵器谱 -之偏微分方程 1】

什么时候显式格式会失败?

显式格式不能任意取时间和空间的网格点数,即 M 与 N 不能随意取值。我们称显式格式为条件稳定。特别地,需要满足所谓 CFL 条件(Courant–Friedrichs–Lewy):

量化分析师的 Python 日记【第 10 天 Q Quant 兵器谱 -之偏微分方程 1】

例如:

*M = 2500

*N = 25

则:

量化分析师的 Python 日记【第 10 天 Q Quant 兵器谱 -之偏微分方程 1】

*M = 1200

*N = 25

则:

量化分析师的 Python 日记【第 10 天 Q Quant 兵器谱 -之偏微分方程 1】

下面的代码计算在第二种情形下的网格点计算过程:

量化分析师的 Python 日记【第 10 天 Q Quant 兵器谱 -之偏微分方程 1】

我们可以通过下图看到,在CFL 条件无法满足的情况下,数值误差累计的结果(特别注意后面的锯齿):

量化分析师的 Python 日记【第 10 天 Q Quant 兵器谱 -之偏微分方程 1】

量化分析师的 Python 日记【第 10 天 Q Quant 兵器谱 -之偏微分方程 1】

今天的日记到此为止,这个问题我们会在下一篇中进行讨论,引出无条件稳定格式:隐式差分格式(Implicit)。

量化分析师的 Python 日记【第 10 天 Q Quant 兵器谱 -之偏微分方程 1】

点击阅读原文,登陆或者注册,克隆源代码到自己的研究环境,就可以自己上手啦!

来源链接:mp.weixin.qq.com