数值分析中，Crank-Nicolson方法是有限差分方法中的一种，用于数值求解热方程以及形式类似的偏微分方程。它在时间方向上是隐式的二阶方法，数值稳定。该方法诞生于20世纪，由John Crank与Phyllis Nicolson发展。

        对于扩散方程（包括许多其他方程），可以证明Crank-Nicolson方法无条件稳定。但是，如果时间步长与空间步长平方的比值过大（一般地，大于1/2），近似解中将存在虚假的振荡或衰减。基于这个原因，当要求大时间步或高空间分辨率的时候，往往会采用数值精确较差的后向欧拉方法进行计算，这样即可以保证稳定，又避免了解的伪振荡。

        这个方法本质上就是用不同的差分来近似$u_{xx}$和$u_t$。这里$u_t$使用向后差分，$u_{xx}$使用混合差分。

$$\frac{U_j^{m+1}-U_j^m}{\tau}=\frac{1}{2}\left( \frac{U_{j+1}^m-2U_j^m+U_{j-1}^m}{h^2}+ \frac{U_{j+1}^{m+1}-2U_j^{m+1}+U_{j-1}^{m+1}}{h^2}\right)$$

按照时间的顺序去排列，可以得到下面的形式

$$-\mu U_{j-1}^{m+1}+(2+2\mu)U_j^{m+1}-\mu U_{j+1}^{m+1}=\mu U_{j-1}^m+(2-2\mu)U_j^m+\mu U_{j+1}^m$$

写成矩阵的形式为$Aw^{m+1}=Bw^m$,这里

$$A=\left(\begin{matrix}2+2\mu & -\mu & 0 & \cdots & 0\\-\mu & 2+2\mu & -\mu & \ddots & 0 \\0 & -\mu & 2+2\mu & \ddots & 0 \\\vdots & \ddots & \ddots & \ddots & \vdots \\0 & \cdots & 0 & -\mu & 2+2\mu\end{matrix}\right)$$


$$B=\left(\begin{matrix}2-2\mu & \mu & 0 & \cdots & 0\\\mu & 2-2\mu & \mu & \ddots & 0 \\0 & \mu & 2-2\mu & \ddots & 0 \\ \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & \cdots & 0 & \mu & 2-2\mu\end{matrix}\right) $$

一维热传导方程为：

$$ \frac{\partial u}{\partial t}=c\frac{\partial^2 u}{\partial x^2} $$

这里$u(0,x) = sin(\pi x); \quad u(t, 0)=u(t,1)=0; \quad c=1$

import numpy as np
import matplotlib.pyplot as plt
from scipy.sparse import diags
from scipy.linalg import solve

Nx = 35
h = 1 / (Nx - 1)

N_temporales = 400
# eps = h**2
eps = 0.001

s = eps / h**2

T = np.zeros((N_temporales, Nx))

T[0] = np.sin(np.pi * np.arange(Nx)*h)

M_derecha = diags((s, (2-2*s), s), (-1, 0, 1), shape=(Nx-2, Nx-2))
M_izquierda = diags((-s, (2 + 2*s), -s), (-1, 0, 1), shape=(Nx-2, Nx-2))

for i in range(1, N_temporales):
    b = np.matmul(M_derecha.toarray(), T[i-1, 1:-1])
    T[i, 1:-1] = solve(M_izquierda.toarray(), b)

x = np.arange(Nx) * h
time = np.arange(N_temporales) * eps

plt.figure(1)
plt.clf()

for i in range(0, N_temporales, 25):
    label = "T={:.3f}".format(eps * i)
    plt.plot(x, T[i], label=label)
    plt.pause(0.2)  

plt.xlabel('x')
plt.ylabel('Temperatura')
plt.legend()
plt.show()

plt.figure(2)
plt.clf()

X, Time = np.meshgrid(x, time)
plt.pcolormesh(X, Time, T)

plt.xlabel("x")
plt.ylabel("tiempo [unidades]")

plt.show()

plt.figure(3)
plt.clf()

plt.plot(time, T[:, int(0.5*(Nx-1))])
plt.xlabel("tiempo")
plt.ylabel('Temperatura')

plt.show()