The one-dimensional diffusion equation is:

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

We will discretize the second-order derivative with a Central Difference scheme: a combination of Forward Difference and Backward Difference of the first derivative.

$$u_{i+1}=u_i+\triangle x\frac{\partial u}{\partial x}\big|_i+\frac{\triangle x^2}{2}\frac{\partial^2 u}{\partial x^2}\big|_i+\frac{\triangle x^3}{3!}\frac{\partial^3 u}{\partial x^3}\big|_i+O{\triangle x^4}$$

$$ u_{i-1}=u_i-\triangle x\frac{\partial u}{\partial x}\big|_i+\frac{\triangle x^2}{2}\frac{\partial^2 u}{\partial x^2}\big|_i-\frac{\triangle x^3}{3!}\frac{\partial^3 u}{\partial x^3}\big|_i+O{\triangle x^4}$$

If we add these two expansions, you can see that the odd-numbered derivative terms will cancel each other out.

$$u_{i+1}+u_{i-1}=2u_i+\triangle x^2\frac{\partial^2 u}{\partial x^2}\big|_i+O(\triangle x^4)$$
the result is

$$ \frac{\partial^2 u}{\partial x^2}=\frac{u_{i+1}-2u_i+u_{i-1}}{\triangle x^2}+O(\triangle x^2) $$


We can now write the discretized version of the diffusion equation in 1D

$$ \frac{u_i^{n+1}-u_i^n}{\triangle t}=\mu \frac{u_{i+1}^n-2u_i^n+u_{i-1}^n}{\triangle x^2} $$

we re-arrange the equation solving for our unknown

$$ u_i^{n+1}=u_i^n+\mu\frac{\triangle t}{\triangle x^2}\left( u_{i+1}^{n}-2u_i^n+u_{i-1}^n \right) $$

we need an initial condition at $t=0 , u=1$ in the interval $0.5\leq x\leq 1$ and $u=2$everywhere else. $\mu=0.3$


In particular, the error bounds can be estimated by the numericalsolutions obtained on grids with different mesh sizes.For example, for the explicit scheme of the heat equation, let $\tau = \mu h^2,\mu \leq \frac{1}{2}$, then we have $U_{(h)j}^m=u_{(h)j}^m+O(h^2)+O(h^4)$

$$U_j^{(1)m}\triangleq \frac{4U_{(\frac{h}{2}2j)}^{4m}-U_{(h)j}^m}{3}=u_{(h)j}^m+O(h^4)$$


import numpy as np
import matplotlib.pyplot as plt

nx = 41 # 网格节点数
dx = 2 / (nx - 1)  # 网格尺寸
nt = 80 # 时间步长
nu = 0.3 # 扩散系数
sigma = 0.2 # 定义中间变量,用于确定时间步长dt
dt = dx**2 * sigma / nu


u = np.ones(nx)
u[int(0.5/dx):int(1/dx)] = 2
un = np.ones(nx)
plt.ion()
plt.show()
for n in range(nt):
    plt.cla() 
    un = u.copy()
    uu = np.ones(nx*2)
    for i in range(1,nx-1):
        u[i] = un[i] + dt * nu / dx**2 * (un[i+1] - 2 * un[i] + u[i-1])
    for i in range(i,nx-1):
        uu[2*i-1]=u[i]
    for i in range(i,nx-1):
        uu[2*i]=(uu[2*i-1]+uu[2*i+1])/2
    for i in range(i,nx-1):
        u[i]=(4*uu[2*i]-uu[2*i-1])/3

    plt.plot(np.linspace(0,2,nx),u,'r',linewidth=3)
    plt.ylim(0.9,2.1)
    plt.title("time: %.4f s"% (n*dt))
    plt.pause(0.2)  
plt.ioff()     
plt.show()