用显格式，向隐格式Crank-Nicolson格式来求解，取h=0.1，r=为0.1进行计算，后并分析误差。（ τ \tau τ为t步长，h为x步长，为m位置，k时刻的温度）
参考课本可得显格式（explicit scheme）的迭代过程： ( u 1 k + 1 u 2 k + 1 ⋮ u m − 2 k + 1 u m − 1 k + 1 ) = ( 1 − 2 r r r 1 − 2 r r ⋱ ⋱ ⋱ r 1 − 2 r r r 1 − 2 r ) ( u 1 k + r u 0 k u 2 k ⋮ u m − 2 k u m − 1 k + r u m k ) \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) ⎝⎜⎜⎜⎜⎜⎛u1k+1u2k+1⋮um−2k+1um−1k+1⎠⎟⎟⎟⎟⎟⎞=⎝⎜⎜⎜⎜⎛1−2rrr1−2r⋱r⋱⋱r1−2rrr1−2r⎠⎟⎟⎟⎟⎞⎝⎜⎜⎜⎜⎜⎛u1k+ru0ku2k⋮um−2kum−1k+rumk⎠⎟⎟⎟⎟⎟⎞ 实现程序：
%参考课本 %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）的迭代过程： ( 1 + 2 r − r − r 1 + 2 r − r ⋱ ⋱ ⋱ − r 1 + 2 r − r − r 1 + 2 r ) ( u 1 k u 2 k ⋮ u m − 2 k u m − 1 k ) \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) ⎝⎜⎜⎜⎜⎛1+2r−r−r1+2r⋱−r⋱−r⋱1+2r−r−r1+2r⎠⎟⎟⎟⎟⎞⎝⎜⎜⎜⎜⎜⎛u1ku2k⋮um−2kum−1k⎠⎟⎟⎟⎟⎟⎞ = ( u 1 k − 1 + r u 0 k u 2 k − 1 ⋮ u m − 2 k − 1 u m − 1 k − 1 + r u m k ) =\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) =⎝⎜⎜⎜⎜⎜⎛u1k−1+ru0ku2k−1⋮um−2k−1um−1k−1+rumk⎠⎟⎟⎟⎟⎟⎞
参照上算法参考课本可得Crank_Nicolsm格式的迭代过程： ( 1 + r − r 2 − r 2 1 + r − r 2 ⋱ ⋱ ⋱ − r 2 1 + r − r 2 − r 2 1 + r ) ( u 1 k + 1 u 2 k + 1 ⋮ u m − 2 k + 1 u m − 1 k + 1 ) \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) ⎝⎜⎜⎜⎜⎛1+r−2r−2r1+r⋱−2r⋱−2r⋱1+r−2r−2r1+r⎠⎟⎟⎟⎟⎞⎝⎜⎜⎜⎜⎜⎛u1k+1u2k+1⋮um−2k+1um−1k+1⎠⎟⎟⎟⎟⎟⎞ = ( 1 − r r 2 r 2 1 − r r 2 ⋱ ⋱ ⋱ r 2 1 − r r 2 r 2 1 − r ) ( u 1 k + r 2 ( u 0 k + u 0 k + 1 ) u 2 k ⋮ u m − 2 k u m − 1 k + r 2 ( u m k + u m k + 1 ) ) =\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) =⎝⎜⎜⎜⎜⎛1−r2r2r1−r⋱2r⋱⋱2r2r1−r1−r2r⎠⎟⎟⎟⎟⎞⎝⎜⎜⎜⎜⎜⎛u1k+2r(u0k+u0k+1)u2k⋮um−2kum−1k+2r(umk+umk+1)⎠⎟⎟⎟⎟⎟⎞ 程序
参照上程序本次实验以CN格式分析为例，按上述计算思想计算得误差数据汇总的得到表1：（ E ∞ = max ∣ u ( x i , t k ) − u i k ∣ \boldsymbol{E}_{\infty}=\max \left|\boldsymbol{u}\left(\boldsymbol{x}_{i}, \boldsymbol{t}_{k}\right)-\boldsymbol{u}_{i}^{k}\right| E∞=max∣∣u(xi,tk)−uik∣∣ ）
本次使用三种不同的差分格式：Crank_Nicolsm格式、显格式（explicit scheme）和隐格式（implicit scheme）对一个详细的一维抛物型方程进行求解，通过对上述格式进行误差分析，可以得到CN格式的求解误差最小，同时改变步长计算误差，得到CN格式下的误差收敛速度为： 。但是实验在三种格式的收敛速度的进一步比较上，限于时间精力，没有进一步展开，有待进一步研究提高。
