一维热传导方程的差分格式

一维热传导方程差分方法的 MATLAB 编程实现。

一维热传导方程

$$\left\{\begin{aligned} &\frac{\partial u}{\partial t}=a \frac{\partial^{2} u}{\partial x^{2}}+f(x, t), && 0<x<1,\ 0<t \leqslant T,\\ & u(x, 0)=\phi(x), && 0 \leqslant x \leqslant 1, \\ & u(0, t)=\alpha(t), u(1, t)=\beta(t), && 0 < t \leqslant T. \end{aligned}\right.$$

其中 $a$ 是正常数，$f(x,t)$、$\phi(x)$、$\alpha(x)$ 和 $\beta(x)$ 为已知函数。$u(x, 0)=\phi(x)$ 为初始条件，$u(0,t)=\alpha(t)$ 和 $u(1,t)=\beta(t)$ 为边界条件。

向前差分格式

以空间步长 $h=1/M$、时间步长 $\tau=T/N$ 分别将 $x$ 轴上区间 $[0,1]$、$t$ 轴上区间 $[0,T]$ 分成 $M$、$N$ 等分，可得

$$x_{i}=j h, \quad 0 \leqslant j \leqslant M,$$ $$t_{n}=n \tau, \quad 0 \leqslant n \leqslant N.$$

一维热传导方程的向前差分格式

$$\frac{u_{j}^{n+1}-u_{j}^{n}}{\tau}=a \frac{u_{j+1}^{n}-2 u_{j}^{n}+u_{j-1}^{n}}{h^{2}}+f(x_j,t_n),$$ $$u_{j}^{0}=\phi_{j}=\phi\left(x_{j}\right), u_{0}^{n}=u_{M}^{n}=0,$$

其中 $j=1,2, \cdots, M-1$，$n=0,1, \cdots, N-1$. 以 $r=a \tau/h^{2}$ 表示网比。

以上格式可改写差分格式

$$u_{j}^{n+1}=r u_{j+1}^{n}+(1-2 r) u_{j}^{n}+r u_{j-1}^{n}+\tau f_{j},$$ $$1 \leqslant j \leqslant M-1,\ 0 \leqslant n \leqslant N-1.$$

先取 $n=0$，利用 $u_{j}^{0}$ 和边值 $u_{0}^{n}=u_{M}^{n}=0$ 算出第一层值 $u_j^1$，再取 $n=2$ ，利用 $u_j^1$ 和边值便可算出 $u_j^2$。如此下去，便可求出所有 $u_j^n$。

设 $\boldsymbol{u}^{n}=(u_{1}^{n}, u_{2}^{n}, \cdots, u_{M-2}^{n}, u_{M-1}^{n})^{\mathrm{T}}, \quad 0 \leqslant n \leqslant N.$

差分格式写成矩阵的形式：

$$\boldsymbol{u}^{n+1}=\boldsymbol{A}\boldsymbol{u}^{n}+\boldsymbol{f}^n, \quad 0 \leqslant n \leqslant N-1.$$

其中矩阵 $\boldsymbol{A}$、向量 $\boldsymbol{f}^n$ 的定义如下，注意向量 $\boldsymbol{f}^n$ 的首尾元素已包含了边界条件。

$$\boldsymbol{A}=\left(\begin{array}{ccccc}1-2 r & r & & & \\ r & 1-2 r & r & & \\ & \ddots & \ddots & \ddots & \\ & & r & 1-2 r & r \\ & & & r & 1-2 r \end{array}\right),$$ $$\boldsymbol{f}^n=\left(\begin{array}{c}\tau f\left(x_{1}, t_{n}\right)+r u_{0}^{n} \\ \tau f\left(x_{2}, t_{n}\right) \\ \vdots \\ \tau f\left(x_{M-2}, t_{n}\right) \\ \tau f\left(x_{M-1}, t_{n}\right)+r u_{M}^{n}\end{array}\right).$$

数值例子

$$\left\{\begin{aligned} &\frac{\partial u}{\partial t}=\frac{\partial^{2} u}{\partial x^{2}}, \quad 0<x<1, \quad 0 <t \leqslant 1, \\ &u(x, 0)=\mathrm{e}^{x}, \quad 0 \leqslant x \leqslant 1, \\ &u(0, t)=\mathrm{e}^{t}, \quad u(1, t)=\mathrm{e}^{1+t}, \quad 0<t \leqslant 1. \end{aligned}\right.$$

方程的真解：$u(x,t)=\mathrm{e}^{x+t}$.

计算数值解

% fdm_heat.m
% forward difference scheme for heat equation
% u_t=u_{xx}, (x,t) in (0,1)x(0,1],
% u(x,0)=exp(x), x in [0,1],
% u(0,t)=exp(t), u(1,t)=exp(1+t), t in (0,1]
% exact solution: u(x,t)=exp(x+t)
clear all; close all;
a=1;
h=0.05; x=[0:h:1];
tau=0.00125; t=[0:tau:1];
r=a*tau/h^2;
M=length(x)-1; N=length(t)-1;
% constructing the coefficient matrix
e=r*ones(M-1,1);
A=spdiags([e 1-2*e e],[-1 0 1],M-1,M-1);
% setting initial and boundary conditions
u=zeros(M+1,N+1);
u(:,1)=exp(x);
u(1,:)=exp(t); 
u(end,:)=exp(1+t);
for n=1:N
    u(2:M,n+1)=A*u(2:M,n);
    u(2,n+1)=u(2,n+1)+r*u(1,n);
    u(M,n+1)=u(M,n+1)+r*u(end,n);
end
% plot the figure
mesh(t(1:20:end),x,u(:,1:20:end))
set(gca,'fontsize',12)
xlabel('t','fontsize', 14) 
ylabel('x','fontsize',14)
zlabel('u','fontsize',14)

% calculating maximum error
[T X]=meshgrid(t,x);
ue=exp(X+T);
Error=max(max(abs(ue-u)))

% print -dpng -r600 fdm_heat.png
% print -depsc2 fdm_heat.eps

输出结果

Error =
   2.1748e-04