抛物线型方程的向前差分格式

---

前言

本文主要介绍求解抛物型方程的最简单显格式（explicit format），也就是向前差分法，包括其理论基础及相应的MATLAB代码实现。

---

一、问题介绍

本文主要研究最简单的抛物型方程的有界域上的初边值问题，以常系数热传导方程：$$Lu=\frac{\partial u}{\partial t}-a\frac{{\partial }^{2}u}{\partial x^2}=f,a为常数，(x,t)\in \Omega ，\Omega 为（0<x<l)\times (0<t\le T)$$

对应的初边值条件为：
初始条件：$u{|}_{t=0}=\phi (x),0\le x\le l$
边界条件：$u{|}_{x=0}=\alpha (t),u{|}_{x=l}=\beta (t),0\le t\le T$

二、向前差分格式

1.网格剖分

对区域$\Omega$进行矩形网格剖分：用平行直线族$x=jh(j=0,1,2,\dots ,M),y=kt(k=0,1,2,\dots ,N)$对区域$\Omega$作矩形网格剖分，网点$(x_j,t_k)简记为(j,k)$,其中$t和h$为正常数，分别成为时间步长和空间步长，且$h=\frac{1}{M},t=\frac{1}{N}$。

2.差商代替导数

将上述方程利用Taylor展开，用差商去代替导数，可以得到以下方程：
$$u_h^{k+1}=[(1-2r)I+rC]u_h^k+t{f}_h^k$$
其中$u_h^k=[u_1^k,u_2^k,...,u_{M-1}^k{]}^T,f_h^k=[f_1^k,f_2^k,...,f_{M-1}^k{]}^T$
$C=\left[\begin{array}{ccccc}0& 1& & & \\ 1& 0& 1& & \\ & \ddots & \ddots & \ddots & \\ & & 1& 0& 1\\ & & & 1& 0\end{array}\right]$$为M-1阶矩阵,r=\frac{at}{h^2}称为网比$

3.截断误差及稳定性

向前差分格式的截断误差为$R_j^k=O(h^2+t)$，且可以通过矩阵方法或者傅里叶方法得到当$网比r\le \frac{1}{2}$时，向前差分格式是稳定的。

---

三.实例及MATLAB代码实现

例题：$$Lu=\frac{\partial u}{\partial t}-\frac{1}{16}\cdot \frac{{\partial }^{2}u}{\partial x^2}=0$$ Ω = { ( x , t ) ∣ 0 < x < 1 ; 0 < t < 1 } \Omega=\left \{ (x,t)|0 <x<1;0 <t<1\right \} Ω={(x,t)∣0<x<1;0<t<1}初边值条件如下：
$$\{\begin{array}{cc}& u(x,0)=2sin(2\pi x),0<x<1\\ & u(0,t)=u(1,t)=0,0<t<1\end{array}$$

对于此问题的数值解的求解方法，本文用MATLAB实现如下：

function [Uh,e,r]=simplestformat(f,M,N,a,b,c,d,k,c1,c2)
%向前差分格式
%f为对应的函数f（x,t)
%M,N分别为x轴和t轴方向的剖分份数
%对应的区域为[a,b]×[c,d]
%k为方程中对应的常数a
%c1,c2为对应的边值条件
%Uh为数值解，e为误差，r为网比
%步长
h=(b-a)/M; 
t=(d-c)/N;
r=k*t/h^2;%网比
if r>0.5
    error("网比大于0.5，不稳定")
end
Uh=zeros(N+1,M+1);%数值解
uh=zeros(N+1,M+1);%精确解
%计算精确解uh
for i=1:N+1
    for j=1:M+1
        uh(i,j)=2*exp((-pi^2)*((i-1)*t)/4)*sin(2*pi*(j-1)*h);%
    end 
end
%赋边界条件
Uh(2:N+1,1)=c1;
Uh(2:N+1,M+1)=c2;
%赋初始条件
Uh(1,1:M+1)=fg(a:h:b);
%计算内点上的数值解Uh
for i=2:N+1
    for j=2:M
        Uh(i,j)=(1-2*r)*Uh(i-1,j)+r*(Uh(i-1,j-1)+Uh(i-1,j+1))+t*f;
    end
end
e=norm(Uh-uh);%误差

初值计算函数：

function y=fg(x)%初值条件的处理
 y=2.*sin(2*x.*pi);%问题一
end

如果需要作出数值解的图像，还可以附上下面的代码：

%数值解绘图
x=a:h:b;
y=c:t:d;
[X,Y]=meshgrid(x,y);
mesh(X,Y,Uh);
colormap("hot")

此处给出M=N=5时的运算结果如下：

e=0.0949720366610601

相应的数值解图像为：

如果需要实现向后差分格式和CN格式数值算法，可以在本算法的基础上，进行修改即可；同时如果需要用本代码求解其他问题，只需要将其中的精确解函数、右端函数$f$以及用于初值计算的函数$fg$作出相应的修改即可。

由于本人水平有限，如有错误敬请批评指正！

You can also refer to my know-it-all blog链接: link