问题内容

用显格式，向隐格式Crank-Nicolson格式来求解，取h=0.1，r=为0.1进行计算，后并分析误差。（$\tau$为t步长，h为x步长，为m位置，k时刻的温度）

算法求解

显格式（explicit scheme）求解

参考课本可得显格式（explicit scheme）的迭代过程： $\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}1-2 r & r & & & \\ r & 1-2 r & r & & \\ & \ddots & \ddots & \ddots & \\ & & & r & 1-2 r & r \\ & & & & r & 1-2 r\end{array}\right)\left(\begin{array}{c}u_{1}^{k}+ru_{0}^k \\ u_{2}^{k} \\ \vdots \\ u_{m-2}^{k} \\ u_{m-1}^{k}+ru_{m}^k\end{array}\right)$ 实现程序：

%参考课本
%slove program by explicit scheme
%u(j,n+1)=u(j,n)+v*(u(j+1,n)-2*u(j,n)+u(j-1,n))
xl=0;
xr=1;
j=10;
dx=(xr-xl)/j;%步长
tf=0.1;
Nt=50;
dt=tf/Nt;
mu=dt/(dx)^2;
%make sure dt satisy stability condition
if mu>0.5
    error('mu shuold<0.5!')
end
%initial condition
x=xl:dx:xr;%grind point
f=sin(pi*x)+sin(2*pi*x);
%store the solution at all grid points for all time steps
u=zeros(j+1,Nt);
u_ture=zeros(j+1,Nt);
%find the approximate solution at each time step
for n=1:Nt
    t=n*dt;
    %boundary condition at left side
    gl=exp(-1*pi*pi*t).*sin(pi*xl)+exp(-4*pi*pi*t).*sin(2*pi*xl);
    %boundary condition at right side
    gr=exp(-1*pi*pi*t).*sin(pi*xr)+exp(-4*pi*pi*t).*sin(2*pi*xr);
    if n==1
        for i=2:j
            u(i,n)=f(i)+mu*(f(j+1)-2*f(j)+f(j-1));
        end
        u(1,n)=gl;
        u(j+1)=gr;
    else
        for i=2:j
           u(i,n)=u(i,n-1)+mu*(u(i+1,n-1)-2*u(i,n-1)+u(i-1,n-1));
        end
        u(1,n)=gl;
        u(j+1)=gr;
    end
    %calculate the analytic solution
    for i=1:j+1
    xi=xl+(i-1)*dx;
    u_ture(i,n)=exp(-1*pi*pi*t).*sin(pi*xi)+exp(-4*pi*pi*t).*sin(2*pi*xi);
    end

end
%plot the result
tt=dt:dt:Nt*dt;
figure(1)
colormap(gray);
surf(x,tt,u');
xlabel('x');
ylabel('t');
zlabel('u');
title('explicit scheme')

%polt the analytic result
figure(2)
colormap(jet);
surf(x,tt,u_ture');
xlabel('x');
ylabel('t');
zlabel('u');
title('analytic solution')

隐格式（implicit scheme）求解

参考课本可得隐格式（implicit scheme）的迭代过程： $\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)\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}u_{1}^{k-1}+r u_{0}^{k} \\ u_{2}^{k-1} \\ \vdots \\ u_{m-2}^{k-1} \\ u_{m-1}^{k-1}+r u_{m}^{k}\end{array}\right)$

参照上算法

Crank_Nicolsm格式求解

参考课本可得Crank_Nicolsm格式的迭代过程： $\left(\begin{array}{ccccc}1+r & -\frac{r}{2} & & & \\ -\frac{r}{2} & 1+r & -\frac{r}{2} & & \\ & \ddots & \ddots & \ddots & \\ & & -\frac{r}{2} & 1+r & -\frac{r}{2} \\ & & & -\frac{r}{2} & 1+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}1-r & \frac{r}{2} & & & \\ \frac{r}{2} & 1-r & \frac{r}{2} & & \\ & \ddots & \ddots & \ddots & \\ & & & \frac{r}{2} & 1-r & \frac{r}{2} \\ & & & \frac{r}{2} & 1-r\end{array}\right)\left(\begin{array}{c}u_{1}^{k}+\frac{r}{2}\left(u_{0}^{k}+u_{0}^{k+1}\right) \\ u_{2}^{k} \\ \vdots \\ u_{m-2}^{k} \\ u_{m-1}^{k}+\frac{r}{2}\left(u_{m}^{k}+u_{m}^{k+1}\right)\end{array}\right)$ 程序

参照上程序

实验结果

三格式的迭代求解结果

误差分析

本次实验以CN格式分析为例，按上述计算思想计算得误差数据汇总的得到表1：（$\boldsymbol{E}_{\infty}=\max \left|\boldsymbol{u}\left(\boldsymbol{x}_{i}, \boldsymbol{t}_{k}\right)-\boldsymbol{u}_{i}^{k}\right|$ ）

实验总结

本次使用三种不同的差分格式：Crank_Nicolsm格式、显格式（explicit scheme）和隐格式（implicit scheme）对一个详细的一维抛物型方程进行求解，通过对上述格式进行误差分析，可以得到CN格式的求解误差最小，同时改变步长计算误差，得到CN格式下的误差收敛速度为： 。但是实验在三种格式的收敛速度的进一步比较上，限于时间精力，没有进一步展开，有待进一步研究提高。