有限差分法求解偏微分方程MATLAB

南京理工大学

课程考核论文

课程名称： 高等数值分析 论文题目： 有限差分法求解偏微分方程 * 名： * * 学 号： ************ 成 绩：

任课教师评语： 签名： 年 月 日 有限差分法求解偏微分方程

一、主要内容

1.有限差分法求解偏微分方程，偏微分方程如一般形式的一维抛物线型方程：

u2u2f(x,t)tx(其中为常数)

具体求解的偏微分方程如下：

u2utx20u(x,0)sin(x)u(0,t)u(1,t)00x1

t02.推导五种差分格式、截断误差并分析其稳定性；

3.编写MATLAB程序实现五种差分格式对偏微分方程的求解及误差分析； 4.结论及完成本次实验报告的感想。

二、推导几种差分格式的过程：

有限差分法（finite-difference methods）是一种数值方法通过有限个微分方程近似求导从而寻求微分方程的近似解。有限差分法的基本思想是把连续的定解区域用有限个离散点构成的网格来代替；把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似；把原方程和定解条件中的微商用差商来近似，积分用积分和来近似，于是原微分方程和定解条件就近似地代之以代数方程组，即有限差分方程组，解此方程组就可以得到原问题在离散点上的近似解。

推导差分方程的过程中需要用到的泰勒展开公式如下：

f(x0)f(x0)f(n)(x0)2f(x)f(x0)(xx0)(xx0)......(xx0)no((xx0)n1)1!2!n! (2-1)

求解区域的网格划分步长参数如下：

tk1tk (2-2) xxhk1k2.1 古典显格式

2.1.1 古典显格式的推导

由泰勒展开公式将u(x,t)对时间展开得 u(xi,t)u(xi,tk)(当ttk1时有

u(xi,tk1)u(xi,tk)(u)i,k(tk1tk)o((tk1tk)2)t (2-4)

uu(xi,tk)()i,ko(2)tu)i,k(ttk)o((ttk)2) (2-3) t

得到对时间的一阶偏导数

u(xi,tk1)u(xi,tk)u()i,k=o() (2-5) t由泰勒展开公式将u(x,t)对位置展开得

u12u u(x,tk)u(xi,tk)()i,k(xxi)(2)i,k(xxi)2o((xxi)3) (2-6)

x2!x当xxi1和xxi1时，代入式(2-6)得

u12uu(xi1,tk)u(xi,tk)()i,k(xi1xi)(2)i,k(xi1xi)2o((xi1xi)3)x2!x  (2-7) 2u(x,t)u(x,t)(u)(xx)1(u)(xx)2o((xx)3)i1kiki,ki1ii,ki1ii1ix2!x2因为xk1xkh，代入上式得

u12uu(xi1,tk)u(xi,tk)()i,kh(2)i,kh2o(h3)x2!x  (2-8) 2u(x,t)u(x,t)(u)h1(u)h2o(h3)i1kiki,ki,kx2!x2得到对位置的二阶偏导数

u(xi1,tk)2u(xi,tk)u(xi1,tk)2u2o(h) (2-9) (2)i,k2xh将式(2-5)、(2-9)代入一般形式的抛物线型偏微分方程得

u(xi,tk1)u(xi,tk)u(x,t)2u(xi,tk)u(xi1,tk)2i1kf(x,t)o(h) ik2h(2-10)

为了方便我们可以将式(2-10)写成

uik1uik uik1uikuik12uikuik1kfi (2-11) 2hh2uki12uikuik1fik (2-12)

最后得到古典显格式的差分格式为

uik1(12ra)uikruik1uik1fik (2-13)

其中rh2，古典显格式的差分格式的截断误差是o(h2)。

2.1.2 古典显格式稳定性分析

古典显格式(2-13)写成矩阵形式为

uk12raIraCukkh1hfh 其中rh2,ukkkkkh(u1,u2,......,uN2,uN1)。

010C10100 101010(N1)(N1)上面的C矩阵的特征值是：C2cos(jh)j1,2,......,N1 H12raIraC

Hj12raraC=12ra2racos(jh)

12ra1cos(jh) 14rasin2jh2j1,2,......,N1使(H)1，即

114rasin2jh21 0ra12

结论：当0ra12时，所以古典显格式是稳定的。

2.2 古典隐格式

(2-14)

(2-15)

2.2.1 古典隐格式的推导 将ttk1代入式 (2-3)得

u)j,k(tk1tk)o((tk1tk)2) (2-16) tuu(xj,tk1)u(xj,tk)()j,ko(2) (2-17)

t得到对时间的一阶偏导数

u(xj,tk1)u(xj,tk)(u(xj,tk)u(xj,tk1)u()j,k=o() (2-18) t将式(2-9)、(2-18)原方程得到

u(xj,tk)u(xj,tk1)u(xj1,tk)2u(xj,tk)u(xj1,tk)2f(x,t)o(h) jk2h(2-19)

为了方便把(2-19)写成

k1ukjujk1ukjujkkukj12ujuj1kfj (2-20) 2hh2ukj1kk2ukjuj1fj (2-21)

最后得到古典隐格式的差分格式为

kkk1k (12ra)ukruuufjj1j1jj (2-22)

其中rh2，古典隐格式的差分格式的截断误差是o(h2)。

2.2.2 古典隐格式稳定性分析

将古典隐格式(2-22)写成矩阵形式如下

k1kk 12raIraCuufhhh(rh2) (2-23)

误差传播方程

k1k (2-24) 12raIraCvhvhA12raIraC,BI所以误差方程的系数矩阵为

1HA112raIraC

jH使(H)1，显然

112ra2racosjhj1,2,......,N1

jH112ra2racos(jh)112ra(1cos(jh))114rasin2jh2

结论：对于r0，即任意网格比下，古典隐格式是绝对稳定的。

jH1恒成立。

2.3 Richardson格式

2.3.1 Richardson格式的推导 将ttk1和ttk1，代入式(2-3)得

uu(x,t)u(x,t)()i,k(tk1tk)o((tk1tk)2)ik1ikt (2-25) uu(x,t)u(x,t)()(tt)o((tt)2)iki,kk1kk1kik1t即

u2u(x,t)u(x,t)()o()ik1iki,kt  (2-26)

uu(x,t)u(x,t)()o(2)iki,kik1t由此得到可得 (u(xi,tk1)u(xi,tk1)u)i,ko(2) (2-27) t2将式(2-9) 、(2-27)代入原方程得到下式

u(xi,tk1)u(xi,tk1)u(x,t)2u(xi,tk)u(xi1,tk)22i1kf(x,t)o(h) ik22h (2-28)

为了方便可以把式(2-28)写成

uik12uikuik1uik1uik1kfi (2-29) 22h即

uik1uik1h2uki12uikuik1fik (2-30)

最后得到Richardson显格式的差分格式为

uik12ruik12uikuik1uik12fik (2-31)

其中rh2，古典显格式的差分格式的截断误差是o(2h2)。

2.3.2 Richardson稳定性分析

将Richardson显格式(2-31)写成如下矩阵形式

k1kk1uh2rC2Iuhuh2fhk (2-32)

误差传播方程矩阵形式

k1kk1vh2rC2Ivhvh (2-33) kkvhvh再将上面的方程组写成矩阵形式

wk1hvk12ra(C2I)Ikkw (2-34) I0v系数矩阵的特征值是

4racos(jh)4ra1 j1028rasin2解得特征值为

8rasin2jh10 (2-35) 21,2jhjhr2a2sin4422 (2-36)

2jhjh+1+16r2a2sin41 (恒成立) (2-37) 22Max1,2=4rasin2结论：上式对任意的网比都恒成立，即Richardson格式是绝对不稳定的。

4. Crank-Nicholson格式

3.4.1 Crank-Nicholson格式的推导 将ttk1代入式(2-9)得

2u2u(x,t)u(x,t)()(tt)o((tt))111jkjkj,kkkkk222t (2-40) 2u(x,t1)u(x,t)(u)jkjk+1j,k1(tk1tk+1)o((tk1tk+1))222t即

uu(x,t)u(x,t)()o(2)1jkjkj,k2t2 (2-41) u(x,t1)u(x,t)(u)o(2)jkjk+1j,k12t2得到如下方程

2u(x,t1)u(x,t)jkjku2()=o(2)tj,k (2-42) 2u(x,t1)u(x,t)ujkjk12o(2)()j,k1=t所以ttk1处的一阶偏导数可以用下式表示：

2)u(xj,tk)2u(xj,tk1)u(xj,tk1)1uu12u(xj,tk122()()j,k1j,k2tt2u(xj,tk1)u(xj,tk)o(2)o()2 (2-43) 将ttk,xxj1和xxj1代入式(2-6)可以得到式(2-9)； 同理ttk1,xxj1和xxj1代入式(2-6)可以得到

u(xj1,tk1)2u(xj,tk1)u(xj1,tk1)2u2(2)j,k1o(h) (2-44) 2xh

所以ttk1，xj处的二阶偏导数用式(2-6)、(2-44)表示：

2

12u2u()()2j,k2j,k12xx1u(xj1,tk1)2u(xj,tk1)u(xj1,tk1)u(xj1,tk)2u(xj,tk)u(xj1,tk)2o(h)222hh

(2-45)

所以ttk1，xj处的函数值可用下式表示：

21  (2-46) f(xj,tk1)f(xj,tk)2原方程变为：

1k11uu2u2uk2()()()()ffo(h) (2-47) j,k1j,kj,kj,k1jj222tt2xx2将差分格式代入上式得：

u(xj,tk1)u(xj,tk)u(xj1,tk1)2u(xj,tk1)u(xj1,tk1)u(xj1,tk)2u(xj,tk)u(xj1,tk)2h2h2122f(x,t)f(x,t)o(h)jk1jk2 (2-48)

为了方便写成：

1k11k11kkk ukj1ukjrukukfjk (2-49) j12ujj1uj12ujuj1fj22

最后得到Crank-Nicholson格式的差分格式为

1rukj1其中rrk1rk1k11kkkuj1uk=1ruuu+ffj1jj1j1jj (2-50) 222h2，Crank-Nicholson格式的差分格式的截断误差是o(2h2)。

3.4.1 Crank-Nicholson稳定性分析

Crank-Nicholson格式写成矩阵形式如下：

rk1rk1k1(1r)ICu=(1r)ICuh+fhfhk (2-51) h222误差传播方程是：

rk1rkCvh=(1r)ICvh (2-52) (1r)I22可以得到：

rr H(1r)IC(1r)IC (2-53)

221jH

(1r)racos(jh)(1r)racos(jh)jh (2-) 12rasin22jh12rasin22使(H)1 即

jh21 (2-55) 1jh12rasin22jhjhjh 12rasin2 (2-56) 12rasin212rasin222212rasin22jh2jh12rasin12rasin22  (2-57)

jhjh12rasin212rasin2222jh22rasin2rasin2 rasin2jh1rasin22jh2 (2-58) jh2上式恒成立。

结论：Crank-Nicholson格式对任意网格比也是绝对稳定的。

5. Du Fort Frankle格式（Richardson格式的改进）

1将uik(uik1uik1)代入式(2-31)并化简得到Du Fort Frankle：

2 uik12ruik1uik1uik1uik1uik12fik (2-59) (12ra)uik12ruik1uik1(12ra)uik12fik (2-60) 可以证明Du Fort Frankle格式是绝对稳定的。因为此格式是Richardson格式的改进格式，因此截断误差还是o(2h2)。

3.6 总结

(1) 古典显格式

古典显格式的差分格式为

uik1(12ra)uikruik1uik1fik

截断误差：o(h2)。 稳定性：当0ra(2) 古典隐格式

古典隐格式的差分格式为

kkk1(12ra)ukfjk jruj1uj1uj1时，古典显格式稳定。 2截断误差：o(h2)。

稳定性：任意网格比古典隐格式绝对稳定。 (3) Richardson显格式

Richardson显格式的差分格式为

uik12ruik12uikuik1uik12fik

截断误差：o(2h2)。

稳定性：任意网格比Richardson格式绝对不稳定。

(4) Crank-Nicholson格式

Crank-Nicholson格式的差分格式为

1rukj1rk1rk1k11kkkuj1uk=1ruuu+ffj1jj1j1jj 222截断误差：o(2h2)。

稳定性：Crank-Nicholson格式对任意网格比绝对稳定。

(5) Du Fort Frankle格式

(12ra)uik12ruik1uik1(12ra)uik12fik 截断误差：o(2h2)。

稳定性：Du Fort Frankle格式对任意网格比绝对稳定。

三、MATLAB实现五种差分格式对偏微分方程的求解及误差分析

3.1 精确数值解

u2utx20u(x,0)sin(x)u(0,t)u(1,t)00x1

t0上述偏微分方程的精确解是

u(x,t)e2tsin(x)(0x1,t0)

区域取值范围：0x1,0t0.2。

用MATLAB对精确解进行编程画出三维图像精确解程序如下：

close all clc

[x,t]=meshgrid(0:0.01:1,0:0.001:0.2)

u=exp(-pi*pi*t).*sin(pi*x) mesh(x,t,u); surf(x,t,u);

xlabel('x'),ylabel('t'),zlabel('u'); title('精确数值解'); shading interp

图3-1 精确数值解的三维图

(a) 精确数值解X-Y平面 (b) 精确数值解X-Z平面

(c) 精确数值解Y-Z平面

图3-2 精确数值解的三个平面图

3.2 五种差分格式MATLAB程序

3.2.1 古典显格式

close all clc T=0.2 X=1.0 M=41 N=11

u=zeros(M,N); %构造一个M行N列的矩阵用于存放时间t和变量x ra=(T/(M-1))/((X/(N-1))^2); %网格比 fprintf('稳定性系数 S=ra 为:\\n');

disp(ra); % 显示网格比 for i=2:N-1

xx=(i-1)*(X/(N-1));

u(1,i)=sin(pi*xx);

end; %即t=0时刻赋值，边界条件 for k=1:M u(k,1)=0; u(k,N)=0;

end; % x=0,x=1处的边界条件

for k=1:M-1 %矩阵是从y轴表示行k，x轴表示列i。由行开始 for i=2:N-1

u(k+1,i)=((1-2*ra)*u(k,i)+ra*u(k,i+1)+ra*u(k,i-1)); %此处为古典显格式 end end

disp(u); % 显示差分法求得的值

[x,t]=meshgrid(1:N,1:M); %将区域划分成网格对每个点赋值再画图 surf(x,t,u); xlabel('x'),ylabel('t'),zlabel('u');

title(' 古典显格式'); %此程序得到的是图3-3

图3-3古典显格式程序结果图（r=0.5）

图3-4精确数值解、古典显格式程序结果的Y-Z平面图（r=0.5）

图3-5古典显格式在取不同网格比时的误差传播结果图

图3-6 不同时间取值时精确解、与古典显格式的值对比图（网格比r=0.5）红线表示精确解、

蓝色线表示差分格式的解

图3-1、图3-3对比可以看出，精确解和古典显格式（网格比r=0.5）的图形是一致的。图3-4精确数值解、古典显格式的Y-Z平面图结果可以看出古典显格式的边界值和精确解一样。图3-5是r分别取0.245、0.5、0.72、1.125时的误差传播图像，边界位置网格数为5处的误差为0.01得到的，可以看出r小于等于0.5是稳定的；但是r大于0.5出现不稳定现象。很好的验证了古典显格式稳定性。

3.2.2 古典隐格式

close all clc T=0.2 X=1.0 M=41 N=21

ra=(T/(M-1))/((X/(N-1))^2); %网格比 fprintf('稳定性系数 S=ra 为:\\n');

disp(ra); %显示网格比

u=zeros(M,N); %构造一个M行N列的矩阵用于存放时间t和变量x for i=2:N-1

xx=(i-1)*(X/(N-1));

u(1,i)=sin(pi*xx); %t=0时刻的赋值，边界条件 end;

for k=1:M u(k,1)=0; u(k,N)=0;

end; % x=0,x=1处的边界条件

A=zeros(N-1,N-1); %隐格式的矩阵形式中的A矩阵赋值 A(1,1)=1+2*ra;

A(N-1,N-1)=1+2*ra;

A(1,2)=-ra; A(N-1,N-2)=-ra; for m=2:N-2

A(m,m-1)=-ra; A(m,m)=1+2*ra; A(m,m+1)=-ra; end;

%以下是——追赶法求u值

d=zeros(N-1,1); %隐格式右边初始矩阵 n=length(d); U=zeros(n); L=eye(n); y=zeros(n,1); x=zeros(n,1); for i=1:N-1

d(i,1)=sin(pi*(i-1)*(1.0/(N-1)));

end %隐格式右边矩阵赋值 %以下循环对矩阵A进行LU分解 U(1,1)=A(1,1); for i=2:n

L(i,i-1)=A(i,i-1)/U(i-1,i-1);

U(i-1,i)=A(i-1,i);%U的上次对角线即为A的上次对角线 U(i,i)=A(i,i)-L(i,i-1)*U(i-1,i); end

for k=1:M-1 %外层大循环是求时间网格2,3……，M的求解u %以下是追赶法的求解过程

%---------追的过程--------即Ly=d的求解y y(1,1)=d(1,1); for i=2:n

y(i,1)=d(i,1)-L(i,i-1)*y(i-1,1);

end

%------赶的过程---------即Ux=y的求解x x(n,1)=y(n,1)/U(n,n); for i=n-1:-1:1

x(i,1)=(y(i,1)-U(i,i+1)*x(i+1,1))/U(i,i);

end %追赶法结束

for i=1:n

u(k+1,i)=x(i,1) end

d=zeros(N-1,1); %更新右边矩阵

d=x %每次外循环更换右边矩阵 end

for k=1:M u(k,1)=0; end;

disp(u); % 显示差分法求得的值

[x,t]=meshgrid(1:N,1:M); %将区域划分成网格对每个点赋值再画图 surf(x,t,u);

xlabel('x'),ylabel('t'),zlabel('u');

title('古典隐格式'); %此程序得到图是3-7

图3-7古典隐格式程序结果图（r=2）

图3-8精确数值解、古典隐格式程序结果的Y-Z平面图（r=2）

图3-9古典隐格式在取不同网格比时的误差传播结果图

图3-10 不同时间取值时精确解、与古典隐格式的值对比图（网格比r=2）红线表示精确解、

蓝色线表示差分格式的解

图3-7古典隐格式在r=2的图形与精确解是一致的。图3-8精确数值解、古典隐格式的Y-Z平面图结果可以看出古典隐格式在t=0.2处的值的边界值和精确解还是有误差的。图3-5是r分别取0.245、0.5、0.72、1.125时的误差传播图像，边界位置网格数为5处的误差为0.01得到的，可以看出r取不同的值时都是稳定的；即古典隐格式对任意的网格比稳定性。从图3-10可以看出隐格式随着时间

的增加，差分格式计算的结果和精确结果越来越大；隐格式虽然对任意网格比都是稳定的，但是计算的精确度是它的不足。

3.2.3 Richardson显格式

close all clc T=0.2 X=1.0000 M=41 N=11

ra=(T/(M-1))/((X/(N-1))^2);

fprintf('稳定性系数 S=ra 为:\\n'); disp(ra);

u=zeros(M,N); %构造一个M行N列的矩阵 for i=2:N-1

xx=(i-1)*(X/(N-1));

u(1,i)=sin(pi*xx); %边界条件 end;

for k=1:M u(k,1)=0;

u(k,N)=0; % 边界条件 end;

%因为Richadson格式需要知道前两行的值，第二行值我们采用隐格式求得

%下面是隐格式求解第二行，和3.2.2隐格式的程序一样，只是求一行，此处不再做注释

A=zeros(N-1,N-1); A(1,1)=1+2*ra;

A(N-1,N-1)=1+2*ra;

A(1,2)=-ra; A(N-1,N-2)=-ra; for m=2:N-2

A(m,m-1)=-ra; A(m,m)=1+2*ra; A(m,m+1)=-ra; end; n=N-1; U=zeros(n); L=eye(n); y=zeros(n,1); x=zeros(n,1); U(1,1)=A(1,1); for i=2:n

L(i,i-1)=A(i,i-1)/U(i-1,i-1);

U(i-1,i)=A(i-1,i); %U的上次对角线即为A的上次对角线 U(i,i)=A(i,i)-L(i,i-1)*U(i-1,i); end

y(1,1)=u(1,1); for i=2:n

y(i,1)=u(1,i)-L(i,i-1)*y(i-1,1); end

x(n,1)=y(n,1)/U(n,n); for i=n-1:-1:1

x(i,1)=(y(i,1)-U(i,i+1)*x(i+1,1))/U(i,i); end

u(2,1)=0; for i=2:N-1

u(2,i)=x(i,1) end

%通过隐格式已求得第二行的值u(2,i)，下面是Richardson格式的过程 for k=2:M-1 %时间轴第3行开始到第M行 for i=2:N-1 %位置网格点

u(k+1,i)=2*ra*(u(k,i-1)-2*u(k,i)+u(k,i+1))+u(k-1,i); %Richardson格式 end end

disp(u); %显示求解的值

[x,t]=meshgrid(1:N,1:M); %区域划分进行赋值画图 surf(x,t,u);

xlabel('x'),ylabel('t'),zlabel('u');

title('Richardson格式'); %此程序得到的图形是图3-11

图3-11 Richardson显格式程序结果图（r=0.5）

图3-12精确数值解、Richardson显格式程序结果的Y-Z平面图（r=0.5）

图3-13 Richardson显格式在取不同网格比时都不稳定的结果图

图3-11是 Richardson显格式在r=0.5时的程序结果图，可以看出不稳定。

图3-12是精确数值解、Richardson显格式程序结果的Y-Z平面图。从图3-13可以看出Richardson显格式在取不同网格比时都出现不稳定现象，和理论结果相一致。所以说Richardson显格式绝对不稳定，这种差分格式不能使用。后面有改进的格式D-F格式。

3.2.4 Crank-Nicholson格式

close all clc T=0.2 X=1.0000 M=41 N=21

ra=(T/(M-1))/((X/(N-1))^2);

fprintf('稳定性系数 S=ra 为:\\n'); disp(ra);

u=zeros(M,N); %构造一个M行N列的矩阵用于存放时间t和变量x%disp(u); for i=2:N-1

xx=(i-1)*(X/(N-1));

u(1,i)=sin(pi*xx); % 边界条件 end;

for k=1:M u(k,1)=0;

u(k,N)=0; % 边界条件 end;

%Crank-Nicholson格式需要两个矩阵，下面的A、B，参照Crank-Nicholson格式的矩阵形式

A=zeros(N-1,N-1); %下面对A进行填充赋值 A(1,1)=1+ra;

A(N-1,N-1)=1+ra; A(1,2)=-ra/2; A(N-1,N-2)=-ra/2; for m=2:N-2

A(m,m-1)=-ra/2; A(m,m)=1+ra; A(m,m+1)=-ra/2; end;

B=zeros(N-1,N-1); %下面对B矩阵进行填充赋值 B(1,1)=1-ra;

B(N-1,N-1)=1-ra; B(1,2)=ra/2; B(N-1,N-2)=ra/2; for m=2:N-2

B(m,m-1)=ra/2; B(m,m)=1-ra; B(m,m+1)=ra/2; end;

%以下是——追赶法求u值

d=zeros(N-1,1); %首先填充右边向量，然后进行L、U分解 n=length(d); U=zeros(n); L=eye(n); y=zeros(n,1); x=zeros(n,1); U(1,1)=A(1,1); for i=2:n

L(i,i-1)=A(i,i-1)/U(i-1,i-1);

U(i-1,i)=A(i-1,i); %U的上次对角线即为A的上次对角线 U(i,i)=A(i,i)-L(i,i-1)*U(i-1,i); end

for i=1:N-1

d(i,1)=sin(pi*(i-1)*(1.0/(N-1))); end

for k=1:M-1 %按层外围大循环，即时间步长 m=zeros(n,1);

m(1,1)=B(1,1)*d(1,1)+B(1,2)*d(2,1);

m(N-1,1)=B(N-1,N-2)*d(N-2,1)+B(N-1,N-1)*d(N-1,1); for i=2:N-2

m(i,1)=B(i,i-1)*d(i-1,1)+B(i,i)*d(i,1)+B(i,i+1)*d(i+1,1);

end %以上是右边矩阵的填充更新

%---------追--------求解Ly=b，其中b是原来的列向量矩阵乘上B系数矩阵得到的 y(1,1)=m(1,1); for i=2:n

y(i,1)=m(i,1)-L(i,i-1)*y(i-1,1); end

%---------赶--------求解Ux=y x(n,1)=y(n,1)/U(n,n); for i=n-1:-1:1

x(i,1)=(y(i,1)-U(i,i+1)*x(i+1,1))/U(i,i); end

for i=1:n

u(k+1,i)=x(i,1) end

d=zeros(N-1,1); %右边向量

d=x end

for k=1:M u(k,1)=0; end;

disp(u); % u的值全部求出，以下画图 [x,t]=meshgrid(1:N,1:M); surf(x,t,u);

xlabel('x'),ylabel('t'),zlabel('u');

title(' Crank-Nicholson格式'); %此程序得到的图像是图3-14

图3-14 Crank-Nicholson格式程序结果图（r=2）

图3-15精确数值解、Crank-Nicholson格式程序结果的Y-Z平面图（r=2）

图3-16 Crank-Nicholson格式在取不同网格比时的误差传播结果图

图3-17 不同时间取值时精确解、与C-N格式的解对比图（网格比r=2）红线表示精确解、

蓝色线表示差分格式的解

图3-14是程序运行得到的Crank-Nicholson格式在网格比取r=2的结果，和精确解图像一致。在t=0.2时从图3-15 的Y-Z面可以看出结果还是有一定的误差。理论上Crank-Nicholson格式对任意的网格比也是稳定的，从图3-16可以看出在r取0.245、0.5、0.72、1.125误差传播图像可以看出误差不会扩撒。图3-17是不同时间取值时精确解、与C-N格式r=2时的解对比图，计算还是具有误差。

3.2.5 Du Fort Frankle格式

close all clc T=0.2 X=1.0000 M=41 N=21

ra=(T/(M-1))/((X/(N-1))^2);

fprintf('稳定性系数 S=ra 为:\\n'); disp(ra);

u=zeros(M,N); %构造一个M行N列的矩阵 for i=2:N-1

xx=(i-1)*(X/(N-1));

u(1,i)=sin(pi*xx); %i表示x轴，k表示y轴(即t) end;

for k=1:M %其实初始矩阵已经将i=1和i=N列的初值赋零了 u(k,1)=0; u(k,N)=0; end;

%第二行用隐格式求得

A=zeros(N-1,N-1); %下面对A进行填充赋值 A(1,1)=1+2*ra;

A(N-1,N-1)=1+2*ra;

A(1,2)=-ra; %第一行第二个和最后一行倒数第二个一个赋值 A(N-1,N-2)=-ra; for m=2:N-2

A(m,m-1)=-ra; A(m,m)=1+2*ra; A(m,m+1)=-ra; end; n=N-1; U=zeros(n); L=eye(n); y=zeros(n,1); x=zeros(n,1); U(1,1)=A(1,1); for i=2:n

L(i,i-1)=A(i,i-1)/U(i-1,i-1);

U(i-1,i)=A(i-1,i); %U的上次对角线即为A的上次对角线 U(i,i)=A(i,i)-L(i,i-1)*U(i-1,i); end

%---------追-------- y(1,1)=u(1,1); for i=2:n

y(i,1)=u(1,i)-L(i,i-1)*y(i-1,1); end

%------赶--------

x(n,1)=y(n,1)/U(n,n); for i=n-1:-1:1

x(i,1)=(y(i,1)-U(i,i+1)*x(i+1,1))/U(i,i); end

u(2,1)=0; for i=2:N-1 u(2,i)=x(i,1)

end %第二行求出，下面用D-F差分格式 for k=2:M-1 for i=2:N-1

u(k+1,i)=(2*ra*u(k,i-1)+2*ra*u(k,i+1)+(1-2*ra)*u(k-1,i))/(1+2*ra); end end disp(u);

[x,t]=meshgrid(1:N,1:M); surf(x,t,u);

xlabel('x'),ylabel('t'),zlabel('u');

title('Du Fort Frankle格式'); %此程序为Richardson格式的改进，得到图3-18

图3-18 Du Fort Frankle格式程序结果图（r=2）

图3-19精确数值解、Du Fort Frankle格式程序结果的Y-Z平面图（r=2）

图3-20 Du Fort Frankle格式在取不同网格比时的误差传播结果图

图3-21 不同时间取值时精确解、与D-F格式的解对比图（网格比r=2）红线表示精确解、

蓝色线表示差分格式的解

D-F格式是Richardson格式（绝对不稳定）的改进，从图3-18可以看出当r时D-F格式是稳定的；图3-20表示D-F格式网格比r为0.245、0.5、0.72、1.125时误差传播图像，不同的网格比误差都不会扩撒。说明这种格式是稳定的，和理论上的结果相一致。图3-21是不同时间取值时精确解、与D-F格式的解对比，随时间的增加计算的值和差分得到的值有误差。此种格式虽然是绝对稳定的，但是计算的精度还是有待提升。

四、结论及感想

从程序得出的结果与精确解的对比来看和理论是上的结果基本一致。比如古典显格式网格比r小于等于0.5（偏微分方程的系数a取1）才稳定，从MATLAB编程运行的结果也可以看出r小于等于0.5是稳定的，r大于0.5不稳定。又如Richardson格式理论上是绝对不稳定的，从编程的结果在取不同的网格比Richardon格式都是不稳定的，理论和结果一致。理论上对不稳定格式进行改进使其稳定，比如得到的D-F格式，取不同的格式网格比都是稳定的，很好的验证了改进的格式的稳定性。

本次报告首先推导了一维抛物线型偏微分方程的一般差分格式。分别是古典显格式、古典隐格式、Richardson显格式、C-N格式进版的Richardson格式D-F格式。推导中得到各种格式的截断误差，从理论上分析了各种格式的稳定性。对于不稳定的格式进行修改得到稳定的格式，即Richardson格式的修改。

通过MATLAB编程实现了各种格式的程序实现。用实验的方法来验证理论结果，能更好的理解各种差分格式的稳定性及操作过程。报告中的程序都是基本程序，误差图与二维x-u的图像（见附录）都是在基本程序的基础上对参数的修改得到的图像。

通过这次报告的完成学到了很多的东西。首先，对编程有了进一步的了解；尤其是使用MATLAB编程。在这个过程中也遇到了很多的问题，通过查阅资料并利用网络资源寻求解决办法。其次，对于差分格式在程序上的实现并不是简单的书写，需要步步衔接好；比如好几种格式都用到差分格式的矩阵形式，尤其是隐式格式不能直接求解，需要应用追赶法进行求解；编写过程中通过直接求解得出的结果都不正确，通过追赶法才能得到正确的结果。最后，我们专业是电磁场与微波技术且偏计算，经常遇到偏微分方程的求解；所以这次试验毫无疑问的对我们理解有限差分法和求解方程提供了很大的帮助。

最后，感谢张老师在课堂上的知识的讲解；同时也感谢在完成报告期间对我提供帮助的同学！

附 录

误差传播三维图（以显格式（图3-5）为例）：

close all clc T=0.2 X=1.0 M=41 N=8

u=zeros(M,N); 阵用于存放时间 %t构造一个和变量x

M行N列的矩ra=(T/(M-1))/((X/(N-1))^2); fprintf(' %网格比 disp(ra); 稳定性系数u(1,5)=0.01;

% 显示网格比 S=ra 为:\\n'); %for 即 t=0k=1:M 时刻赋值，误差 0.01 end; u(k,1)=0; for k=1:M-1 u(k,N)=0;

% x=0,x=1处的边界条件 %行开始矩阵是从

y

轴表示行k，x轴表示列i。由 u(k+1,i)=((1-2*ra)*u(k,i)+ra*u(k,i+1)+ra*u(k, for i=2:N-1

i-1)); %此处为古典显格式end

end disp(u); [x,t]=meshgrid(1:N,1:M); % 显示差分法求得的值 %subplot(2,2,1)

将区域划分成网格对每个点赋值再画图

surf(x,t,u); xlabel('x'),ylabel('t'),zlabel('u');

title(' 古典显格式-误差传播 r=0.245'); clc T=0.2 X=1.0 M=41 N=11

u=zeros(M,N);

ra=(T/(M-1))/((X/(N-1))^2);

fprintf('disp(ra);

稳定性系数 S=ra 为:\\n'); u(1,5)=0.01; % for k=1:M 即 t=0 时刻赋值，误差 0.01 end;

u(k,1)=0; u(k,N)=0; for k=1:M-1 u(k+1,i)=((1-2*ra)*u(k,i)+ra*u(k,i+1)+ra*u(k, for i=2:N-1

i-1)); end

end disp(u); [x,t]=meshgrid(1:N,1:M);

subplot(2,2,2)

surf(x,t,u);

xlabel('x'),ylabel('t'),zlabel('u');

title('clc 古典显格式-误差传播 r=0.5'); T=0.2 X=1.0 M=41 N=13

u=zeros(M,N);

ra=(T/(M-1))/((X/(N-1))^2);

fprintf('disp(ra);

稳定性系数 S=ra 为:\\n'); u(1,5)=0.01; %for 即t=0时刻赋值，误差0.01 k=1:M end;

u(k,1)=0; u(k,N)=0; for k=1:M-1 u(k+1,i)=((1-2*ra)*u(k,i)+ra*u(k,i+1)+ra*u(k, for i=2:N-1

i-1)); end

end disp(u); [x,t]=meshgrid(1:N,1:M);

subplot(2,2,3) surf(x,t,u); xlabel('x'),ylabel('t'),zlabel('u');

title(' 古典显格式-误差传播 r=0.72'); clc T=0.2 X=1.0 M=41 N=16

u=zeros(M,N);

ra=(T/(M-1))/((X/(N-1))^2); fprintf('disp(ra);

稳定性系数 S=ra 为 :\\n'); u(1,5)=0.01; %for k=1:M 即 t=0 时刻赋值，误差 0.01 end;

u(k,1)=0; u(k,N)=0; for k=1:M-1 u(k+1,i)=((1-2*ra)*u(k,i)+ra*u(k,i+1)+ra*u(k, for i=2:N-1

i-1)); end end disp(u);

[x,t]=meshgrid(1:N,1:M); subplot(2,2,4) surf(x,t,u); xlabel('x'),ylabel('t'),zlabel('u');

title(' 古典显格式-误差传播 r=1.125');

不同时刻的二维图的程序（基础程序运行后运行如下程序得到二维图）： 以显格式（图3-6）为例：

k=11; k=26; xx=(0:X/(N-1):1) xx=(0:X/(N-1):1) a=(xx-(X/(N-1)))/(X/(N-1))+2 a=(xx-(X/(N-1)))/(X/(N-1))+2 i=int8(a) i=int8(a) y=u(k,i) y=u(k,i) x=(0:0.01:1) x=(0:0.01:1) z=exp(-pi*pi*((k-1)*(T/(M-1)))).*sin(pi*x) z=exp(-pi*pi*((k-1)*(T/(M-1)))).*sin(pi*x) t=((k-1)*(T/(M-1))) t=((k-1)*(T/(M-1))) subplot(2,2,1) subplot(2,2,3) plot(xx,y,'b',x,z,'r') xlabel('x'),ylabel('u'); title('t=0.05 线对比'); 时，精确解、古典显格式的解曲 k=20; xx=(0:X/(N-1):1) a=(xx-(X/(N-1)))/(X/(N-1))+2 i=int8(a) y=u(k,i) x=(0:0.01:1) z=exp(-pi*pi*((k-1)*(T/(M-1)))).*sin(pi*x) t=((k-1)*(T/(M-1))) subplot(2,2,2) plot(xx,y,'b',x,z,'r') xlabel('x'),ylabel('u'); title('t=0.095 曲线对比'); 时，精确解、古典显格式的解

plot(xx,y,'b',x,z,'r') xlabel('x'),ylabel('u'); title('t=0.125 曲线对比'); 时，精确解、古典显格式的解k=36; xx=(0:X/(N-1):1) a=(xx-(X/(N-1)))/(X/(N-1))+2 i=int8(a) y=u(k,i) x=(0:0.01:1) z=exp(-pi*pi*((k-1)*(T/(M-1)))).*sin(pi*x) t=((k-1)*(T/(M-1))) subplot(2,2,4) plot(xx,y,'b',x,z,'r') xlabel('x'),ylabel('u'); title('t=0.175 曲线对比');时，精确解、古典显格式的解