《2022年东南大学数值分析第七章偏微分方程数值解法.pdf》由会员分享，可在线阅读，更多相关《2022年东南大学数值分析第七章偏微分方程数值解法.pdf（5页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、第七章 偏微分方程数值解法 Crank-Nicolson格式* （学号） *（姓名）上机题目要求见教材P346,10 题。一、算法原理本文研究下列定解问题 (抛物型方程 ) 22( , ) (0,0)( ,0)( ) (0)(0, )( ),(1, )( ) (0)uuaf x txltTtxu xxxluttutttT(1) 的有限差分法，其中a为正常数，, ,f为已知函数，且满足边界条件和初始条件。关于式 (1)的求解，采用离散化方法，剖分网格，构造差分格式。其中，网格剖分是将区域0,0DxltT 用两簇平行直线 (0) (0)ikxxihiMttkkN分割成矩形网格，其中,lThMN分别
2、为空间步长和时间步长。将式(1)中的偏导数使用不同的差商代替， 将得到不同的差分格式， 如古典显格式、 古典隐格式、Crank-Nicolson 格式等。其中， Crank-Nicolson 格式具有更高的收敛阶数，应用更广泛，故本文采用Crank-Nicolson 格式求解抛物型方程。Crank-Nicolson 格式推导： 在节点(,)2ikx t处考虑式 (1)，有22(,)(,)(,)222ikikikuux tax tf x ttx(2) 对偏导数(,)2ikux tt用中心差分展开2311+131(,)(,)(,)(,) ()224kkikikikiikikuux tu x tu
3、x txtttt(3) 将22(,)2ikux tx在节点(,)ikx t和1(,)ikx t表示为精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 1 页，共 5 页 - - - - - - - - - - 222+122224+1221(,)=(,)+(,)22(,) ()8ikikikkkiikikuuux tx tx txxxuxttxt(4) 对以上两个偏导数用二阶差分展开2112224i+141(,)(,)2 (,)(,)(,) ()12ikikikikkkikiiux tu xtu x tu
4、xtxhhutxxx(5) 211111122241i+141(,)(,)2 (,)(,)(,) ()12ikikikikkkikiiux tu xtu x tu xtxhhutxxx(6) 将式(4)(5)(6)分别代入式 (3)，略去高阶小量，用kiu代替(,)iku x t并化简得2111111112122,22kkkkkkkkiiiiiiiiikauuuuuuuufx th (7) 令2/rah，将式 (7)联合式 (1)初始条件和边界条件，用矩阵的形式表示为：11112212211111221122221122221122kkkkkkMMkkMMrrrruurrrrrruurrrru
5、urruurrrr112211,22,2,2,22kkkkMkMkkkrfx tttfx tfxtrfxttt(8) Crank-Nicolson 格式的截断误差为22()ROh，具有较高的精度。精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 2 页，共 5 页 - - - - - - - - - - 二、计算代码Crank_Nicolson 格式完整代码functionU=Crank_Nicolson(f,a,x0,xn,dx,t0,tm,dt,g,s0,sn) %采用 Crank_Nicolson格式
6、求解抛物线型偏微分方程% 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
7、*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:M C=dt*feval(f,x,t(k)+0.5*dt); C=C; C(1)=
8、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)xttuuxttxu xexuteutet(9) 精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 3 页，共 5 页 - - - - - - - - - - 取空间步长1/xM
9、 ，时间步长1/tN ，采用 Crank_Nicolson格式计算，并将计算结果与精确值( )x ty xe作比较，计算调用程序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,
10、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;%步长反复二分endfor ii=2:kkerr(ii)/err(ii-1)%误差观察endf
11、printf(%16.15f , err)在 MATLAB 中运行以上程序，计算结果如表1 所示表 1 Crank_Nicolson格式计算结果及其误差(40MN) 输出点 (, x t) 数值解精确解误差(0.2,1.0) 3.320148266101440 3.320116922736547 0.000031343364893 (0.4,1.0) 4.055249987095087 4.055199966844675 0.000050020250413 (0.6,1.0) 4.953086121831125 4.953032424395115 0.000053697436010 (0.8,
12、1.0) 6.049686221438103 6.049647464412947 0.000038757025156 从表 1 可以看出，在空间步长和时间步长均为0.25 时， 采用 Crank_Nicolson格式计算抛物线方程 (5)的数值解与精确解之间的误差达到了410， 故很好的实现偏微分方程的数值解法， 且具有较高的精度和较低的运算复杂度。同时，为了比较不同的步长所带来的计算精度，本文通过将步长反复二分， 观察四个输出点的误差均方和减小的规律，如表2 所示。表 2 不同步长下的误差及其规律精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳
13、- - - - - - - - - -第 4 页，共 5 页 - - - - - - - - - - k 步长误差均方和kerr误差减少倍数1/kkerrerr1 0.25 0.000088712758987 2 0.125 0.000022180089438 0.250021413953775 3 0.0625 0.000005545141552 0.250005373866960 4 0.03125 0.000001386290165 0.250000861545108 5 0.015625 0.000000346561277 0.249991874539774 6 0.0078125 0
14、.000000086561708 0.249773167584743 从表 2 可以看出，随着空间步长和时间步长的减小，数值计算的误差呈下降趋势，说明，随着时间步长和空间步长的减小，计算结果可以收敛到精确解。同时，通过观察表 2 最后一列，可见，时间步长和空间步长同时减半，计算误差降低为原来的四分之一，这与Crank_Nicolson 格式的截断误差22()ROh（其中为时间步长， h为空间步长）是符合一致的。四、结论本文介绍了偏微分方程求解的Crank_Nicolson 方法，其主要思想是用离散的、只含有限个未知数的差分方程去近似代替连续变量的微分方程和边界条件，并将相应的差分方程的解作为原问题的近似解。Crank_Nicolson 格式是一种隐格式，是绝对稳定的。同时，本文考察了其收敛性，实验证明，Crank_Nicolson格式是收敛的，且具有二阶收敛性。精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 5 页，共 5 页 - - - - - - - - - -