东南大学_数值分析_第七章_偏微分方程数值解法

东南大学《数值分析》上机练习——算法与程序设计实验报告第七章偏微分方程数值解法——Crank-Nicolson格式****（学号）*****（姓名）上机题目要求见教材P346,10题。一、算法原理本文研究下列定解问题(抛物型方程)22(,)(0,0)(,0)()(0)(0,)(),(1,)()(0)uuafxtxltTtxuxxxluttutttT(1)的有限差分法，其中a为正常数，,,,f为已知函数，且满足边界条件和初始条件。关于式(1)的求解，采用离散化方法，剖分网格，构造差分格式。其中，网格剖分是将区域0,0DxltT用两簇平行直线(0)(0)ikxxihiMttkkN分割成矩形网格，其中,lThMN分别为空间步长和时间步长。将式(1)中的偏导数使用不同的差商代替，将得到不同的差分格式，如古典显格式、古典隐格式、Crank-Nicolson格式等。其中，Crank-Nicolson格式具有更高的收敛阶数，应用更广泛，故本文采用Crank-Nicolson格式求解抛物型方程。Crank-Nicolson格式推导：在节点(,)2ikxt处考虑式(1)，有22(,)(,)(,)222ikikikuuxtaxtfxttx(2)对偏导数(,)2ikuxtt用中心差分展开2311+131(,)(,)(,)(,)()224kkikikikiikikuuxtuxtuxtxtttt(3)将22(,)2ikuxtx在节点(,)ikxt和1(,)ikxt表示为第七章偏微分方程数值解法222+122224+1221(,)=(,)+(,)22(,)()8ikikikkkiikikuuuxtxtxtxxxuxttxt(4)对以上两个偏导数用二阶差分展开2112224i+141(,)(,)2(,)(,)(,)()12ikikikikkkikiiuxtuxtuxtuxtxhhutxxx(5)211111122241i+141(,)(,)2(,)(,)(,)()12ikikikikkkikiiuxtuxtuxtuxtxhhutxxx(6)将式(4)(5)(6)分别代入式(3)，略去高阶小量，用kiu代替(,)ikuxt并化简得2111111112122,22kkkkkkkkiiiiiiiiikauuuuuuuufxth(7)令2/rah，将式(7)联合式(1)初始条件和边界条件，用矩阵的形式表示为：11112212211111221122221122221122kkkkkkMMkkMMrrrruurrrrrruurrrruurruurrrr112211,22,2,2,22kkkkMkMkkkrfxtttfxtfxtrfxttt(8)Crank-Nicolson格式的截断误差为22()ROh，具有较高的精度。东南大学《数值分析》上机练习——算法与程序设计实验报告二、计算代码Crank_Nicolson格式完整代码functionU=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为二阶导系数%-x0andxn空间域端点%-t0andtm时间域端点%-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=[];fork=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)xttuuxttxuxexuteutet(9)第七章偏微分方程数值解法取空间步长1/xM，时间步长1/tN，采用Crank_Nicolson格式计算，并将计算结果与精确值()xtyxe作比较，计算调用程序Crank_Nicolson.m，如下clcformatlongf=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;forii=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;%步长反复二分endforii=2:kkerr(ii)/err(ii-1)%误差观察endfprintf('%16.15f',err)在MATLAB中运行以上程序，计算结果如表1所示表1Crank_Nicolson格式计算结果及其误差(40MN)输出点(,xt)数值解精确解误差(0.2,1.0)3.3201482661014403.3201169227365470.000031343364893(0.4,1.0)4.0552499870950874.0551999668446750.000050020250413(0.6,1.0)4.9530861218311254.9530324243951150.000053697436010(0.8,1.0)6.0496862214381036.0496474644129470.000038757025156从表1可以看出，在空间步长和时间步长均为0.25时，采用Crank_Nicolson格式计算抛物线方程(5)的数值解与精确解之间的误差达到了410，故很好的实现偏微分方程的数值解法，且具有较高的精度和较低的运算复杂度。同时，为了比较不同的步长所带来的计算精度，本文通过将步长反复二分，观察四个输出点的误差均方和减小的规律，如表2所示。表2不同步长下的误差及其规律东南大学《数值分析》上机练习——算法与程序设计实验报告k步长误差均方和kerr误差减少倍数1/kkerrerr10.250.00008871275898720.1250.0000221800894380.25002141395377530.06250.0000055451415520.25000537386696040.031250.0000013862901650.25000086154510850.0156250.0000003465612770.24999187453977460.00781250.0000000865617080.249773167584743从表2可以看出，随着空间步长和时间步长的减小，数值计算的误差呈下降趋势，说明，随着时间步长和空间步长的减小，计算结果可以收敛到精确解。同时，通过观察表2最后一列，可见，时间步长和空间步长同时减半，计算误差降低为原来的四分之一，这与Crank_Nicolson格式的截断误差22()ROh（其中为时间步长，h为空间步长）是符合一致的。四、结论本文介绍了偏微分方程求解的Crank_Nicolson方法，其主要思想是用离散的、只含有限个未知数的差分方程去近似代替连续变量的微分方程和边界条件，并将相应的差分方程的解作为原问题的近似解。Crank_Nicolson格式是一种隐格式，是绝对稳定的。同时，本文考察了其收敛性，实验证明，Crank_Nicolson格式是收敛的，且具有二阶收敛性。