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

一维热传导方程差分方法的 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.$$

其中 是正常数，、、 和 为已知函数。 为初始条件， 和 为边界条件。

向前差分格式

以空间步长 、时间步长 分别将 轴上区间 、 轴上区间 分成 、 等分，可得

$$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,$$

其中 ，. 以 表示网比。

以上格式可改写差分格式

$$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.$$

先取 ，利用 和边值 算出第一层值 ，再取 ，利用 和边值便可算出 。如此下去，便可求出所有 。

设

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

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

其中矩阵 、向量 的定义如下，注意向量 的首尾元素已包含了边界条件。

$$\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.$$

方程的真解：.

计算数值解

% 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