曲柄尼科尔森方案介绍

简介

曲柄尼科尔森方案（Crank-Nicolson method）是求解偏微分方程的一种方法，它是一种隐式的时间迭代方法，由美国物理学家 John Crank 和英国数学家 Phyllis Nicolson 所提出。

该算法本质上将欧拉方法的显式求解变为一个隐式求解，通过更加准确地模拟计算，从而减少误差。

算法

曲柄尼科尔森方案的迭代公式如下：

$$\frac{u^{n+1}j-u^n_j}{\Delta t}=\frac{1}{2}\left(\frac{\partial^2 u}{\partial x^2}\bigg|{n,j+1}+\frac{\partial^2 u}{\partial x^2}\bigg|_{n+1,j}\right)$$

其中，$u^n_j$ 表示第 $n$ 个时间步长和第 $j$ 个空间节点处的解，$\Delta t$ 为时间步长，$\bigg|_{n,j}$ 表示在时间 $n$ 和位置 $j$ 处的导数值。

可以将上述迭代公式表示为矩阵形式：

$$\begin{bmatrix} 1+\alpha & -\frac{\alpha}{2} & 0 &\dots & 0 \ -\frac{\alpha}{2} & 1+\alpha & -\frac{\alpha}{2} & \ddots & \vdots \ 0 & -\frac{\alpha}{2} & 1+\alpha & \ddots & 0 \ \vdots & \ddots & \ddots & \ddots & -\frac{\alpha}{2} \ 0 & \dots & 0 & -\frac{\alpha}{2} & 1+\alpha \end{bmatrix}\begin{bmatrix} u_1 \ u_2 \ \vdots \ \vdots \ u_n \end{bmatrix}=\begin{bmatrix} b_1 \ b_2 \ \vdots \ \vdots \ b_n \end{bmatrix}$$

其中，$\alpha=\frac{\Delta t}{(\Delta x)^2}$，$u_i$ 表示在位置 $i$ 处的解，$b_i$ 表示由上一时间步长计算出的值。

编程实现

以下是使用 Python 编写的曲柄尼科尔森方案求解一维抛物线方程的程序。其中，$f(x)$ 表示初始条件，$u0(x)$ 表示边界条件。

import numpy as np

def crank_nicolson(u0, f, a, b, tn, dx, dt):
    """
    使用曲柄尼科尔森方案求解一维抛物线方程

    :param u0: 指定边界条件
    :param f: 指定初始条件
    :param a: 区间左端点 a
    :param b: 区间右端点 b
    :param tn: 时间区间 [0, tn]
    :param dx: 空间步长
    :param dt: 时间步长
    :return: 计算得到的解
    """
    alpha = dt / (dx ** 2)
    N = int((b - a) / dx)
    M = int(tn / dt)
    u = np.zeros((M + 1, N + 1))
    x = np.linspace(a, b, N + 1)
    t = np.linspace(0, tn, M + 1)
    u[0, :] = f(x)
    u[:, 0] = u0(t)
    u[:, N] = 0  # 右端点固定为0
    A = np.zeros((N - 1, N - 1))
    A[0, 0], A[N - 2, N - 2] = 1 + alpha, 1 + alpha
    A[0, 1] = -alpha / 2
    A[N - 2, N - 3] = -alpha / 2
    for i in range(1, N - 2):
        A[i, i], A[i, i - 1], A[i, i + 1] = 1 + alpha, -alpha / 2, -alpha / 2
    for n in range(1, M + 1):
        b = u[n - 1, 1:N]
        b[0] += alpha / 2 * u[n, 0]
        b[N - 2] += alpha / 2 * u[n, N]
        u[n, 1:N] = np.dot(np.linalg.inv(A), b)
    return x, t, u

调用 crank_nicolson() 函数，并将解绘制成三维图像：

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

a, b, tn, dx, dt = 0, 1, 0.5, 0.1, 0.005

def f(x):
    return np.sin(np.pi * x)

def u0(t):
    return 0

x, t, u = crank_nicolson(u0, f, a, b, tn, dx, dt)

X, T = np.meshgrid(x, t)

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, T, u, cmap='cool')
ax.set_xlabel('x')
ax.set_ylabel('t')
ax.set_zlabel('u')
plt.show()

运行上述代码，得到以下图像：

结论

曲柄尼科尔森方案可以用于求解偏微分方程的数值解，可以较为准确地模拟计算，减少误差。在实际应用中，需要根据具体情况调整参数，以获得更加准确的解。