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

本文最后更新于：1 个月前

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

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

$$ { \[\begin{array}{lcl} \dfrac{\partial{u}}{\partial{t}}=\dfrac{\partial^{2}{u}}{\partial{x}^{2}} &,&0 \leq x \leq 1,0 \leq t \leq 1 \\ u(x,0)=e^x &, & 0 \leq x \leq 1 \\ u(0,t)=e^t,u(1,t)=e^{1+t},&, &0 \leq t \leq 1 \end{array}\]

. $$

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

定义误差为\[ E_{\infty}(h,\tau)=\max \limits_{1 \leq i \leq m-1 \atop 1 \leq k \leq n } |u(x_{i},t_{k})-u_{k}^k)| \]

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

建立差分格式

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

\[x_i=ih,0 \leq i \leq m\]

\[t_k=k \tau,0 \leq k \leq n\]

设\(w=\left\{w_{i} \mid 0 \leqslant i \leqslant m\right\} \in V_{h}.\) 定义算子 \[ (\mathcal{A} w)_{i}=\left\{\begin{array}{ll} \frac{1}{12}\left(w_{i-1}+10 w_{i}+w_{i+1}\right), & 1 \leqslant i \leqslant m-1 \\ w_{i}, & i=0, m \end{array}\right. \] 通过泰勒展式等一系列高端操作,得到如下差分格式 \[ \left\{ \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), \quad 1 \leqslant i \leqslant m-1, \quad 0 \leqslant k \leqslant n-1 \\ u_{i}^{0}=\varphi\left(x_{i}\right), \quad 0 \leqslant i \leqslant m \\ u_{0}^{k}=\alpha\left(t_{k}\right), \quad u_{m}^{k}=\beta\left(t_{k}\right), \quad 1 \leqslant k \leqslant n \end{array} \right. \] 记 \[ 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), \quad 1 \leqslant i \leqslant m-1 \\ u_{0}^{k+1}=\alpha\left(t_{k+1}\right), \quad 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} \\ \vdots \\ 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} \\ \vdots \\ 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) \\ \vdots \\ \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 \leq i \leq m-1,\ 1 \leq k \leq n,\ \gamma=\dfrac{\tau}{h^2},\ f(x_i,t_k)=0.\)

求解程序基于Matlab R2019b

主程序源代码 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
% 开发者: 捉不到小恐龙
% 邮箱: 517093978@qq.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处的数值解和精确解见下图,从图像上看很接近。

t=1处的数值解(\tau=\dfrac{1}{10},h=\dfrac{1}{100})和精确解

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

x=0.5时，不同t处的数值解、精确解和误差

取不同\(\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