第七章 偏微分方程数值解法——Crank-Nicolson 格式****（学号） *****（姓名）上机题目要求见教材P346,10题。

一、算法原理本文研究下列定解问题(抛物型方程)22(,) (0,0)(,0)() (0)(0,)(), (1,)() (0)u ua f x t x l t T t x u x x x l u t t u t t t T ϕαβ⎧∂∂-=<<≤≤⎪∂∂⎪=≤≤⎨⎪==<≤⎪⎩(1)的有限差分法，其中a 为正常数，,,,f ϕαβ为已知函数，且满足边界条件和初始条件。

关于式(1)的求解，采用离散化方法，剖分网格，构造差分格式。

其中，网格剖分是将区域{}0,0D x l t T =≤≤≤≤用两簇平行直线(0)(0)i k x x ih i M t t k k N τ==≤≤⎧⎨==≤≤⎩ 分割成矩形网格，其中,l Th M Nτ==分别为空间步长和时间步长。

将式(1)中的偏导数使用不同的差商代替，将得到不同的差分格式，如古典显格式、古典隐格式、Crank-Nicolson 格式等。

其中，Crank-Nicolson 格式具有更高的收敛阶数，应用更广泛，故本文采用Crank-Nicolson 格式求解抛物型方程。

Crank-Nicolson 格式推导：在节点(,)2i k x t τ+处考虑式(1)，有22(,)(,)(,)222i k i k i k u u x t a x t f x t t x τττ∂∂+-+=+∂∂ (2)对偏导数(,)2i k u x t t τ∂+∂用中心差分展开 []2311+131(,)(,)(,)(,) ()224k k i k i k i k i i k i k u ux t u x t u x t x t t t t ττηητ++∂∂+=--<<∂∂ (3) 将22(,)2i k u x t x τ∂+∂在节点(,)i k x t 和1(,)i k x t +表示为222+122224+1221(,)=(,)+(,)22 (,) ()8i k i k i k k k i i k i k u u ux t x t x t x x x ux t t x tττηη⎡⎤∂∂∂+⎢⎥∂∂∂⎣⎦∂-<<∂∂ (4)对以上两个偏导数用二阶差分展开[]2112224i+141(,)(,)2(,)(,) (,) ()12i k i k i k i k kk i k i i u x t u x t u x t u x t x hh u t x x x ξξ+-∂=-+∂∂-<<∂ (5)[]211111122241i+141(,)(,)2(,)(,) (,) ()12i k i k i k i k k ki k i i u x t u x t u x t u x t x hh u t x x xξξ++++-++∂=-+∂∂-<<∂ (6)将式(4)(5)(6)分别代入式(3)，略去高阶小量，用k i u 代替(,)i k u x t 并化简得()()()2111111112122,22k k k k k k k k i i i i i i i i i k a u u u u u u u u f x t h ττ+++++-+-⎛⎫⎡⎤---++-+=+ ⎪⎣⎦⎝⎭ (7) 令2/r a h τ=，将式(7)联合式(1)初始条件和边界条件，用矩阵的形式表示为：11112212211111221122221122221122 k k k k k kM M k k M M r r r r u u r r r r r ru u r r r r u u r r u u r r r r +++--+--⎡⎤⎡⎤+---⎢⎥⎢⎥⎢⎥⎢⎥⎡⎤⎡⎤⎢⎥⎢⎥-+--⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-+--⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎢⎥⎢⎥⎢⎥⎢⎥-+-⎢⎥⎢⎥⎣⎦⎣⎦()()()()()()112211,22,2 ,2,22k k k k M k M k k k r f x t t t f x t f x t r f x t t t ττααττττττββ+--+⎡⎤⎛⎫+++ ⎪⎢⎥⎝⎭⎢⎥⎢⎥⎛⎫+ ⎪⎢⎥⎝⎭⎢⎥+⎢⎥⎢⎥⎛⎫⎢⎥+ ⎪⎢⎥⎝⎭⎢⎥⎛⎫⎢⎥+++ ⎪⎢⎥⎝⎭⎣⎦(8) Crank-Nicolson 格式的截断误差为22()R O h τ=+，具有较高的精度。

二、计算代码Crank_Nicolson 格式完整代码function U=Crank_Nicolson(f,a,x0,xn,dx,t0,tm,dt,g,s0,sn) %采用Crank_Nicolson 格式求解抛物线型偏微分方程 % du/dt-a*d2u/dx2=f(x,t) %Input - f 抛物方程右端函数 % - a 为二阶导系数 % - x0 and xn 空间域端点 % - t0 and tm 时间域端点 % - dx 为空间步长，dt 为时间步长 % - g 为初始条件函数 % - s0,sn 为边界条件函数%Output - U 输出,横坐标为空间，纵坐标为时间 M=(tm-t0)/dt;N=(xn-x0)/dx;%网格数 x=x0+dx:dx:xn-dx; t=t0:dt:tm; u=feval(g,x);u=u'; r=a*dt/dx/dx;%步长比%Crank-Nicolson 格式：A*u_(k+1)=B*u_k+CA=-r/2*[zeros(1,N-1);eye(N-2,N-2),zeros(N-2,1)]-r/2* ...[zeros(N-2,1),eye(N-2,N-2);zeros(1,N-1)]+(1+r)*eye(N-1,N-1); A=inv(A);B=r/2*[zeros(1,N-1);eye(N-2,N-2),zeros(N-2,1)]+r/2* ...[zeros(N-2,1),eye(N-2,N-2);zeros(1,N-1)]+(1-r)*eye(N-1,N-1); U=[]; for k=1:MC=dt*feval(f,x,t(k)+0.5*dt); C=C';C(1)=C(1)+r/2*(feval(s0,t(k))+feval(s0,t(k+1))); C(N-1)=C(N-1)+r/2*(feval(sn,t(k))+feval(sn,t(k+1))); u=A*(B*u+C); U=[U;u']; end三、计算结果及分析对于定解问题2210 (01,01)(,0) (01)(0,), (1,) (01)x t t u ux t t x u x e x u t e u t e t +⎧∂∂-=<<≤≤⎪∂∂⎪=≤≤⎨⎪==<≤⎪⎩(9)取空间步长1/x M ∆=，时间步长1/t N ∆=，采用Crank_Nicolson 格式计算，并将计算结果与精确值()x t y x e +=作比较，计算调用程序Crank_Nicolson.m ，如下clcformat longf=inline('0*x+0*t');%抛物方程右端函数 g=inline('exp(x)');%初始条件函数s0=inline('exp(t)');sn=inline('exp(t+1)');%边界条件函数 a=1;x0=0;xn=1;t0=0;tm=1; M=40;N=40; err=[]; kk=4; for ii=1:kkdx=1/N;dt=(tm-t0)/M;%空间步长，时间步长U=Crank_Nicolson(f,a,x0,xn,dx,t0,tm,dt,g,s0,sn);%调用Crank_Nicolson 函数 U_CN=U(length(U),N/5:N/5:4*N/5); %在四个目标点处的数值解 fu=inline('exp(x+t)'); %精确方程 xx=0.2:0.2:0.8;tt=1.0;U_real=feval(fu,xx,tt); %在四个目标点处的精确解 %abs(U_CN-U_real)err=[err,sqrt((U_CN-U_real)*(U_CN-U_real)')]; %平方误差和开方 M=M*2;N=N*2;%步长反复二分 end for ii=2:kkerr(ii)/err(ii-1)%误差观察 endfprintf('%16.15f ', err)在MATLAB 中运行以上程序，计算结果如表1所示表1 Crank_Nicolson 格式计算结果及其误差(40M N ==)从表1可以看出，在空间步长和时间步长均为0.25时，采用Crank_Nicolson 格式计算抛物线方程(5)的数值解与精确解之间的误差达到了410-，故很好的实现偏微分方程的数值解法，且具有较高的精度和较低的运算复杂度。

同时，为了比较不同的步长所带来的计算精度，本文通过将步长反复二分，观察四个输出点的误差均方和减小的规律，如表2所示。

表2 不同步长下的误差及其规律从表2可以看出，随着空间步长和时间步长的减小，数值计算的误差呈下降趋势，说明，随着时间步长和空间步长的减小，计算结果可以收敛到精确解。

同时，通过观察表2最后一列，可见，时间步长和空间步长同时减半，计算误差降低为原来的四分之一，这与Crank_Nicolson 格式的截断误差22()R O h τ=+（其中为τ时间步长，h 为空间步长）是符合一致的。

四、结论本文介绍了偏微分方程求解的Crank_Nicolson 方法，其主要思想是用离散的、只含有限个未知数的差分方程去近似代替连续变量的微分方程和边界条件，并将相应的差分方程的解作为原问题的近似解。

Crank_Nicolson 格式是一种隐格式，是绝对稳定的。

同时，本文考察了其收敛性，实验证明，Crank_Nicolson 格式是收敛的，且具有二阶收敛性。