一维非齐次热传导方程的紧致差分格式(附Matlab代码)

一维非齐次热传导方程的紧致差分格式(附Matlab代码)

考虑一维非齐次热传导方程$Dirichlet$初边值问题(第一类边界值问题):

$$\{\begin{array}{lcl}\frac{\partial u}{\partial t}=\frac{{\partial }^{2}u}{\partial x^2}& ,& 0\le x\le 1,0\le t\le 1\\ u(x,0)=e^x& ,& 0\le x\le 1\\ u(0,t)=e^t,u(1,t)=e^{1+t},& ,& 0\le t\le 1\end{array}$$

该问题的精确解为$u(x,t)=e^{x+t}$.

定义误差为$$E_{\infty }(h,\tau )=\underset{\genfrac{}{}{0ex}{}{1\le i\le m-1}{1\le k\le n}}{\max }|u(x_i,t_k)-u_k^k)|$$

针对以上定解问题建立一个具有$O({\tau }^2+h^4)$精度的无条件稳定的紧致差分格式,求上述问题的数值解并对数值解、精度和误差阶进行相应的数值分析.

建立差分格式

将$x$进行$m$等分，将$t$进行$n$等分, 记$h=\frac{1}{m},\tau =\frac{1}{n}$.

$$x_i=ih,0\le i\le m$$

$$t_k=k\tau ,0\le k\le n$$

设 w = { w i ∣ 0 ⩽ i ⩽ m } ∈ V h . w=\left\{w_{i} \mid 0 \leqslant i \leqslant m\right\} \in V_{h}. w={wi​∣0⩽i⩽m}∈Vh​. 定义算子
$$(\mathcal{A}w{)}_i=\{\begin{array}{ll}\frac{1}{12}\left(w_{i-1}+10w_i+w_{i+1}\right),& 1⩽i⩽m-1\\ w_i,& i=0,m\end{array}$$
通过泰勒展式等一系列高端操作,得到如下差分格式
$$\{\begin{array}{l}\mathcal{A}{\delta }_t{u}_i^{k+\frac{1}{2}}-a{\delta }_x^2{u}_i^{k+\frac{1}{2}}=\mathcal{A}f\left(x_i,t_{k+\frac{1}{2}}\right),1⩽i⩽m-1,0⩽k⩽n-1\\ u_i^0=\phi \left(x_i\right),0⩽i⩽m\\ u_0^k=\alpha \left(t_k\right),u_m^k=\beta \left(t_k\right),1⩽k⩽n\end{array}$$
记
$$u^k=\left(u_0^k,u_1^k,\cdots ,u_{m-1}^k,u_m^k\right)$$
由初值条件知第 0 层上的值 $u^0$ 已知. 设已确定出了第 $k$ 层的值 $u^k$. 则关于第 $k+1$ 层 值的差分格式为
$$\begin{array}{l}\mathcal{A}{\delta }_t{u}_i^{k+\frac{1}{2}}-a{\delta }_x^2{u}_i^{k+\frac{1}{2}}=\mathcal{A}f\left(x_i,t_{k+\frac{1}{2}}\right),1⩽i⩽m-1\\ u_0^{k+1}=\alpha \left(t_{k+1}\right),u_m^{k+1}=\beta \left(t_{k+1}\right)\end{array}$$
写成如下矩阵形式
$$\begin{array}{l}\left(\begin{array}{ccccc}\frac{5}{6}+r& \frac{1}{12}-\frac{1}{2}r& & & \\ \frac{1}{12}-\frac{1}{2}r& \frac{5}{6}+r& \frac{1}{12}-\frac{1}{2}r& & \\ & \ddots & \ddots & \ddots & \\ & & \frac{1}{12}-\frac{1}{2}r& \frac{5}{6}+r& \frac{1}{12}-\frac{1}{2}r\\ & & & \frac{1}{12}-\frac{1}{2}r& \frac{5}{6}+r\end{array}\right)\left(\begin{array}{c}{u}_1^{k+1}\\ u_2^{k+1}\\ ⋮\\ u_{m-2}^{k+1}\\ u_{m-1}^{k+1}\end{array}\right)\\ \\ =\left(\begin{array}{ccccc}\frac{5}{6}-r& \frac{1}{12}+\frac{1}{2}r& & & \\ \frac{1}{12}+\frac{1}{2}r& \frac{5}{6}-r& \frac{1}{12}+\frac{1}{2}r& & \\ & \ddots & \ddots & \ddots & \\ & & \frac{1}{12}+\frac{1}{2}r& \frac{5}{6}-r& \frac{1}{12}+\frac{1}{2}r\\ & & \frac{1}{12}+\frac{1}{2}r& \frac{5}{6}-r& \end{array}\right)\left(\begin{array}{c}{u}_1^k\\ u_2^k\\ ⋮\\ u_{m-2}^k\\ u_{m-1}^k\end{array}\right)\\ \\ +\left(\begin{array}{c}\left(\frac{1}{12}+\frac{1}{2}r\right)u_0^k-\left(\frac{1}{12}-\frac{1}{2}r\right)u_0^{k+1}+\tau \mathcal{A}f\left(x_1,t_{k+\frac{1}{2}}\right)\\ \tau \mathcal{A}f\left(x_2,t_{k+\frac{1}{2}}\right)\\ ⋮\\ \left(\frac{1}{12}+\frac{1}{2}r\right)u_m^k-\left(\frac{1}{12}-\frac{1}{2}r\right)u_m^{k+1}+\tau \mathcal{A}f\left(x_{m-1},t_{k+\frac{1}{2}}\right)\end{array}\right)\end{array}$$
其中，$1\le i\le m-1,1\le k\le n,\gamma =\frac{\tau }{h^2},f(x_i,t_k)=0.$

求解程序基于Matlab R2019b

主程序源代码

% 开发者: 图灵的猫
% 邮箱: turingscat@126.com
function [t,x,u] = fsolve(tau,h)
%   用紧致差分格式解求数值解
%   @t 时间向量
%   @x 空间向量
%   @tau 时间步长
%   @h 空间步长
%   @u 数值解
    N=1/tau;%t被分割的区间数
    M=1/h;%x被分割的区间数
    t=0:tau:1;
    x=0:h:1;
    u=ones(N+1,M+1);
    u(1,:)=exp(x);%初值
    %边值
    u(:,1)=exp(t);
    u(:,M+1)=exp(1+t);

    r=tau/h^2;%步长系数
    a1=ones(M-2,1)*(1/12-r/2);%下对角线
    b1=ones(M-1,1)*(5/6+r);%主对角线
    c1=ones(M-2,1)*(1/12-r/2);%上对角线
    A1=diag(b1,0)+diag(a1,-1)+diag(c1,1);%系数矩阵
    
    a2=ones(M-2,1)*(1/12+r/2);%下对角线
    b2=ones(M-1,1)*(5/6-r);%主对角线
    c2=ones(M-2,1)*(1/12+r/2);%上对角线
    A2=diag(b2,0)+diag(a2,-1)+diag(c2,1);%系数矩阵
%     for k=2:N+1
%         %右端项
%         f=A2*u(k-1,2:M)';
%         f(1)=f(1)+r/2*(u(k,1)+u(k-1,1));
%         f(M-1)=f(M-1)+r/2*(u(k,M+1)+u(k-1,M+1));
%         u(k,2:M)=(A1\f)';%由k-1层求第k层的值
%     end
    for k=2:N+1
        %右端项
        f=zeros(M-1,1);
        f(1)=(1/12+r/2)*u(k-1,1)-(1/12-r/2)*u(k,1);
        f(M-1)=(1/12+r/2)*u(k-1,M+1)-(1/12-r/2)*u(k,M+1);
        
        u(k,2:M)=(A1\(A2*u(k-1,2:M)'))'+(A1\f)';%由k-1层求第k层的值
    end
end

程序运行结果

当$\tau =\frac{1}{10},h=\frac{1}{100}$时,t=1处的数值解和精确解见下图,从图像上看很接近.

当$\tau =\frac{1}{10},h=\frac{1}{10}$时,取$x=0.5$, 不同$t$处的值见下表,
当层数越深时，误差越大，这是因为,每一次由$k$层求解$k+1$层时都有误差，随着$t$的增大，误差不断累积，越来越大,到t=1处误差变得最大.

取不同$\tau$ 和 $h$时，t=1处的误差曲线见下图,步长越小，误差也越小.

分别用向后Euler格式和紧致差分格式，取不同步长时最大误差和最大误差的比值见下表.
可知，用向后Euler格式， $\tau$缩小为原来的$\frac{1}{4}$，$h$缩小为原来的$\frac{1}{2}$，误差会缩小为原来的$\frac{1}{4}$，符合$O(\tau +h^2)$的截断误差; 而紧致差分格式，$\tau$缩小为原来的$\frac{1}{4}$，$h$缩小为原来的$\frac{1}{2}$，误差会缩小为原来的$\frac{1}{16}$，符合$O({\tau }^2+h^4)$的截断误差,精度高于向后Euler格式.

若需要绘图及误差分析等代码部分,请在评论区留下邮箱.

若想深入了解紧致差分格式的构造理论,请加好友探讨,互相学习,共同进步!

数值算例来源: 《微分方程数值解》-M.Ran