差分法

不同的差分法方法在计算机上有着不同的解法，更加精确的差分方法也对应着更高的精度（更容易收敛），作为计算机仿真偏微分方程的重要内容，值得记录一下。并且记录下各个方法的\(A\)矩阵，便于后续的使用。（本差分的思想都是依据泰勒公式，且主要是记录几种差分方法，其原理就不做介绍）

差分法

计算机无法处理连续的问题，因此需要借助差分将连续的问题离散化。对于常见的偏微分或者微分方程计算中，差分法的体现就在于用\(dx\)将整个定义域细分，用每个\(dx\)开头、结尾或开头和结尾平均的斜率来代替偏导或导数。当\(dx\)取值合适的时候，能起到很好的计算机仿真作用。

向前差分

我们举例就采用热传导中的傅里叶定律： \[ \frac {\partial T}{\partial t}=\alpha \frac {\partial^2 T}{\partial^2 x} \] 首先，对应情景中的\(x、t\)，我们取\(x=mdx, t=ndt\),依据上述的替代思想，有： \[ \frac{T^{n+1}_m-T^n_m}{dt}=\alpha\frac{T^n_{m+1}-2T^n_m+T^n_{m-1}}{dx^2} \] 向前的关键在于左边的偏导数是用\(\frac{T^{n+1}_m-T^n_m}{dt}\)替代的，取的是\(T^{n+1}_m-T^n_m\)，而不是\(T^n_m-T^{n-1}_m\)。这有什么影响呢？对于向前和向后来说，他们的精度是一样的。但是由于显式和隐式的区别，前者能迭代算出，而后者不行，需要使用追赶法（或者直接matlab把矩阵除过去）。因此他们的计算编程方法不同。

向前差分写成矩阵形式： \[ \begin{bmatrix} T^{n+1}_1\\ T^{n+1}_2 \\ \vdots\\ T^{n+1}_{m-1}\\T^{n+1}_{m} \end{bmatrix} = \begin{bmatrix} T^{n}_1\\ T^{n}_2 \\ \vdots\\ T^{n}_{m-1}\\T^{n}_{m} \end{bmatrix} + \frac{\alpha dt}{dx^2} \begin{bmatrix} -2 & 1 & \cdots & \cdots \\ 1 & -2 & 1 & \cdots \\ & \ddots & \ddots & \ddots \\ \cdots & 1 &-2 &1 \\\ \cdots & \cdots &1 & -2 \end{bmatrix} \begin{bmatrix} T^{n}_1\\ T^{n}_2 \\ \vdots\\ T^{n}_{m-1}\\T^{n}_{m} \end{bmatrix} \] 进行计算机仿真时，只需要设定\(\alpha、dx、dt\)与\(A\)矩阵—— \[ A= \begin{bmatrix} -2 & 1 & \cdots & \cdots \\ 1 & -2 & 1 & \cdots \\ & \ddots & \ddots & \ddots \\ \cdots & 1 &-2 &1 \\\ \cdots & \cdots &1 & -2 \end{bmatrix} \] 即可通过每一个\(T^{n}_{m}\)迭代计算出\(T^{n+1}_{m}\)。

向后差分

有了向前差分的先见，我们直接给出傅里叶定律向后差分的形式： \[ \frac{T^{n}_m-T^{n-1}_m}{dt}=\alpha\frac{T^n_{m+1}-2T^n_m+T^n_{m-1}}{dx^2} \] 这一步向后的关键就在于此时的\(T^{n}_m和T^{n-1}_m\)没有办法写成\(T^{n}_m=XT^{n-1}_m\)，形成迭代关系（因而称为隐式）。这就需要用到不同的计算方法，但重点是在于\(A\)矩阵的变化。

我们设\(\lambda=\frac{\alpha dt}{dx^2}\)此时我们需要将式子改写成： \[ -\lambda T^n_{m-1}+(1+2\lambda)T^n_m-\lambda T^n_{m+1}=T^{n-1}_m \] 即： \[ \begin{bmatrix} 1+2\lambda & -\lambda & \cdots & \cdots \\ -\lambda & 1+2\lambda & -\lambda & \cdots \\ & \ddots & \ddots & \ddots \\ \cdots & -\lambda & 1+2\lambda & -\lambda \\ \cdots & \cdots & -\lambda & 1+2\lambda \end{bmatrix} \begin{bmatrix} T^{n}_1\\ T^{n}_2 \\ \vdots\\ T^{n}_{m-1}\\T^{n}_{m} \end{bmatrix} = \begin{bmatrix} T^{n-1}_1\\ T^{n-1}_2 \\ \vdots\\ T^{n-1}_{m-1}\\T^{n-1}_{m} \end{bmatrix} \] 其中设置： \[ A=\begin{bmatrix} 1+2\lambda & -\lambda & \cdots & \cdots \\ -\lambda & 1+2\lambda & -\lambda & \cdots \\ & \ddots & \ddots & \ddots \\ \cdots & -\lambda & 1+2\lambda & -\lambda \\ \cdots & \cdots & -\lambda & 1+2\lambda \end{bmatrix} \]

看起来像个反显式是吗？把A矩阵除过去或者使用追赶法，就可以迭代出新的温度啦。

二级差分

不知道它的官方叫法，姑且这么喊着。

这种差分法，主要是对等式右边进行离散化的优化。

直接给出式子： \[ \frac{T^{n+1}_m-T^n_m}{dt}=\frac{\alpha}{2}(\frac{T^n_{m+1}-2T^n_m+T^n_{m-1}}{dx^2}+\frac{T^{n+1}_{m+1}-2T^{n+1}_m+T^{n+1}_{m-1}}{dx^2}) \] 可以看到，等式右边分别在\(n和n+1\)时刻取了斜率代替导数，并将二者平均，和前面仅用一处进行替代的方法显得更加精准。因此这也是一种很好的差分方法。

想解这种差分法，关键还是形成\(n+1\)时刻与\(n\)时刻的迭代关系。在等式的右边，二者已经出现，结合左边来看，形成了高度的对称性，因此将他们分别移项和合并同类项即可得（我们设置\(\beta=\frac{\alpha dt}{2dx^2}\)）： \[ \begin{bmatrix} -(1+2\beta) & -\beta & \cdots & \cdots \\ \beta & -(1+2\beta) & -\beta & \cdots \\ & \ddots & \ddots & \ddots \\ \cdots & \beta & -(1+2\beta) & -\beta \\ \cdots & \cdots & \beta & -(1+2\beta) \end{bmatrix} \begin{bmatrix} T^{n+1}_1\\ T^{n+1}_2 \\ \vdots\\ T^{n+1}_{m-1}\\T^{n+1}_{m} \end{bmatrix} = \begin{bmatrix} -1+2\beta & -\beta & \cdots & \cdots \\ -\beta & -1+2\beta & -\beta & \cdots \\ & \ddots & \ddots & \ddots \\ \cdots & -\beta & -1+2\beta & -\beta \\ \cdots & \cdots & -\beta & -(1+2\beta) \end{bmatrix} \begin{bmatrix} T^{n}_1\\ T^{n}_2 \\ \vdots\\ T^{n}_{m-1}\\T^{n}_{m} \end{bmatrix} \]

无非是矩阵多了些，作为一个迭代求解，等式右边是完全已知的。再把等式左边的系数矩阵除过去，即可迭代出新一时刻的温度。

注

具体解相应的热力学问题时，需要把\(A\)矩阵（即各类系数矩阵的一个称呼吧）的第一行和最后一行进行相应边界条件的替换，才能实现真正的使用。

后记

学会离散化的求解，矩阵化的运算，也算踏入计算机运算的门槛了吧。