写在前面：炎炎夏日，用代码给你的crush送上一朵雪花吧。一起感受一下偏微分方程的浪漫。

有限差分法相关参考资料

以下资料为我学习的来源：
链接1: 中科院的有限差分法的详细解释课程
链接2: 论文：Modeling and numerical simulations of dendritic crystal growth
链接3: 知乎 对应上文解释 代码来源！

先上手看代码，然后理解数理概念

要注释

from numpy import *

# 定义参数

# 本模块对应的知识点1【划分网格对于求解二维偏微分方程的作用】
# 网格的行数，列数以及时间迭代的总步数
M,  N,  K   = 300,  300,  4000 
# 空间时长步长
Dx, Dy, Dt  = 0.01, 0.01, 1e-5


# 材料参数
# 包含物理和材料属性参数，包括相场动力学时间常数、各向异性能量强度、相场方向依赖的强度、晶格取向因子、参考取向、界面宽度参数、温度敏感度、平衡温度和热耦合系数。
# 为求解数学物理方程中的控制物理生成物形态的参数
tau = 3e-4
eps_bar = 0.01
sigma = 0.02
J = 4.
theta_0 = 0.2
alpha = 0.9
gamma = 10.
T_eq = 1.
kappa = 1.8


# 初始条件设置
p = zeros((M, N))
T = zeros((M, N))

# 半径为sqrt(5.0)的圆内，被视为初始固化区域的一部分，这部分我们认为相变已经完成，即雪花已经开始产生
for i in range(M):
    for j in range(N):
        if (i - M/2)**2 + (j - N/2)**2 < 5.0:  
            p[i, j] = 1.0 # 已经开始产生相变，设为1

# 定义二维拉普拉斯算子(扩散、波动等现象的偏微分方程时非常重要的算子)
# 知识点2【临近点对于求解偏微分方程的作用】
# 知识点3【有限差分方法中的中心差分公式】
def Lap(p):
    # 提取内部和邻近点的值
    p_i_j  = delete(delete(p, [0, -1], axis=0), [0, -1], axis=1)
    p_im_j = delete(delete(p, [0, -1], axis=0), [-1,-2], axis=1)
    p_ip_j = delete(delete(p, [0, -1], axis=0), [0,  1], axis=1)
    p_i_jm = delete(delete(p, [0, -1], axis=1), [0,  1], axis=0)
    p_i_jp = delete(delete(p, [0, -1], axis=1), [-1,-2], axis=0)
    # 计算拉普拉斯算子,使用限差分方法中的中心差分公式
    Lap_p  = (p_im_j + p_ip_j + p_i_jm + p_i_jp - 4*p_i_j)/Dx**2
    # 处理边界条件
    Lap_pj = vstack((Lap_p[0,:], Lap_p, Lap_p[-1,:]))
    return hstack((Lap_pj[:,0].reshape(N,1), Lap_pj, Lap_pj[:,-1].reshape(N,1)))

# 计算相场p的演化，这里可以看知乎的那篇推导原文，十分详细
def Phase_field(p, T):
    theta = arctan2(gradient(p, Dy, axis=1), gradient(p, Dx, axis=0))
    eps = eps_bar * (1. + sigma * cos(J * (theta - theta_0)))
    g = -eps * eps_bar * sigma * J * sin(J * (theta - theta_0)) * gradient(p, Dy, axis=1)
    h = -eps * eps_bar * sigma * J * sin(J * (theta - theta_0)) * gradient(p, Dx, axis=0)
    m = alpha/pi * arctan(gamma * (T_eq - T))
    # 结合界面能、界面动力学以及热效应等因素，形成相场的演化方程。
    term_1 = - p*(p - 1.)*(p - 0.5 + m)
    term_2 = - gradient(g, Dx, axis=0)
    term_3 = gradient(h, Dy, axis=1)
    term_4 = eps**2 * Lap(p)
    p_ev = Dt / tau * (term_1 + term_2 + term_3 + term_4)
    return p + p_ev

# 计算温度场T的演化
def Temp(T, p_new, p_old):
    T_ev = Dt*Lap(T) + kappa*(p_new - p_old)
    return T + T_ev

# 相场更新
p_hist = []
T_hist = []
p_old = p; T_old = T
for t_step in range(K):
    p_new = Phase_field(p_old, T_old)
    T_new = Temp(T_old, p_new, p_old)
    p_old = p_new
    T_old = T_new
    if t_step % 50 == 0:
        p_hist.append(p_new)
        T_hist.append(T_new)
        print('step finished:', t_step,'/',str(K))

有限差分法的理解

Q: 什么是有限差分法？

在我的理解里，有限差分法就是对于微分方程/偏微分方程求解的一种数理方法，本质就是把对于微分的求解转化为求解大量代数方程组。转化的过程本质上是一种近似。

举一个最简单的例子，假设我们有一个连续的函数，我们想要计算它在某个特定点的导数。为了使用有限差分法，我们将函数在该点附近进行离散化处理。首先，我们选择一个很小的步长，称为差分间距（finite difference interval）或网格间距（grid spacing）。然后，在该点的左右两侧选择一些离散点。我们根据这些离散点的函数值来估计函数在该点的导数。具体来说，我们可以使用中心差分（central difference）或前向差分（forward difference）等方法来计算导数。

代码中涉及的知识点

1. 划分网格对于求解二维偏微分方程的作用

网格划分是数值方法中处理偏微分方程的基础，它将连续的物理空间离散化成有限的点集，使得连续问题可以在计算机上用离散方式处理。
在这段代码中，M和N定义了在两个空间维度上的网格点数，这决定了模型的空间分辨率。网格点越多，空间分辨率越高，模拟结果越精确，但计算量也越大。

2. 临近点对于求解偏微分方程的作用

将求解域进行离散化是偏微分方程的关键。如果是一维条件下进行计算求解，我们可以在脑子里想象一条线，知道后一个点，就可以根据初始值迭代。如果是二维条件下进行计算求解，那么每个点周围的临近点应该会有四个，在脑子里把他们连起来，那么会形成一个矩形，使用任何四个结点的值就可以算出最后一个。这就是为什么二维中心差分公式会使用如下表达式迭代：

此外，代码中的空间步长(Dx, Dy)和时间步长(Dt)决定了数值模拟的精度和稳定性。步长越小，理论上数值解越接近真实解，但计算时间也更长。时间步长还与数值方法的稳定性有关。特别是对于显式方法，如果步长过大，可能导致数值解的不稳定或发散。

3. 有限差分方法中的中心差分公式

中心差分因为相比前向差分/后向差分多一重精度，因此，经常被用来作为差分公式的公式。
Lap§函数计算二维拉普拉斯算子，这是通过有限差分的中心差分公式实现的。拉普拉斯算子在物理中描述了某个量（如温度、压力等）的空间扩散。在边界处采用简单的延伸处理，这相当于假设边界外的值等于边界上的值（Neumann边界条件）。

总结

通过学习使用代码产生雪花，我对于有限差分法有了进一步的理解。
另外标题是吸引大家一起学习有限差分法的，不建议给crush发雪花代码，带他/她去哈尔滨看雪效果更好。