抛物型方程隐格式的程序实现

1淮海工学院实验报告书课程名称：微分方程数值解法实验名称：抛物型方程隐格式的程序实现班级信息101姓名：学号：日期：地点数学实验室指导教师：成绩：数理科学系21.实验目的：按照实验内容对-方法导出的几种隐式差分格式编程实现求出数值解；估计数值解的误差界；大致比较几种不同格式的优劣；结合格式的相容性、稳定性和收敛性条件简单分析计算结果。2.实验内容：用-方法导出的隐式格式求如下初边值问题,10,0),1(,0),0(,10,sin)0,(10,10,022ttutuxxxutxxutu（2）的数值解，观察取不同值时，数值解的误差情况,大致比较至少三种不同格式的优劣，分析解释计算结果。3.实验步骤：第一步:对求解区域作网格剖分。我们把求解区域}10,10|),{(txtxG分割为一些均匀的矩形，其边与坐标轴平行，取空间步长为Jx/1和时间步长为Mt/1，其中J，M为正整数，用两族平行直线：Jjxjxxj...,,2,1,平行于t轴作竖直线和Mntnttn...,,2,1,平行于x轴作水平线,将矩形区域G分割成矩形网格，得网格剖分图。第二步:差分法的目的是：求初边值问题（2）的解),(txu在节点),(njtx的近似值),...,2,1;1,,2,1(MnJjUnj。为此，需要构造逼近微分方程定解问题的差分格式。采用-方法，构造如下差分格式0,0),()(2)1()(2002112111111nJnjjjnjnjnjnjnjnjnjnjUUxUxUUUxUUUtUU以2)/(xtr表示网比。由于本是实验中xxsin)(,所以上述差分格式简化为0,0,sin)1())1(21()1()21(001111111nJnjjnjnjnjnjnjnjUUxUrUUrrUrUUrrU利用Taylor展式，可得差分格式的截断误差为：3))()((])(81)(241[])(121)21[()(43212221221xtOututuxtuuTnjxxtttttnjxxxxxxtnj利用Fourier方法判断差分格式的稳定性条件是：当210，当且仅当1)21(21r时稳定，当121，对所有的r格式均稳定。利用最大值原理可以证明差分格式的收敛条件是：1)1(21r。第三步:有限差分方程的解法。根据问题的规模和计算机的容量和速度，选取适当的解法。这是数值计算的关键一步，有限差分法解方程的主要计算都集中在这里。由于本实验问题的有限差分方程为-方法导出的隐式差分格式，可以用Thomas方法求解。对每个时间层1n需要求解的方程为：1...,,2,1)1())1(21()1()21(1111111JjrUUrrUrUUrrUnjnjnjnjnjnj初边值条件为0,0,sin00nJnjjUUxU。第四步:根据前面的分析，编制程序，上机计算，求出数值解，必要时画出解的几何图形，分析观察解的性质。为了方便比较可以用傅立叶级数方法求出定解问题的解析解为：xetxutsin),(2。另外，上述差分程组为三对角方程组，可用Thomas法（追赶法）求解。若有三对角方程组：njdxcxbxajjjjjjj...,,2,1,11，则Thomas法（追赶法）公式为niaebafdfniaebcebdfbceiiiiiiiiiiii,,3,2,1,,3,2,,111111111，41,,2,1,1nnixefxfxiiiinn其中nxxx...,,,21为所求三对角方程组的解。四．实验数据记录及分析（或程序及运行结果）：程序如下：clearcita=input('cita=');M=input('M=');J=input('J=');u=zeros(J+1,M+1);d=zeros(1,J);e=zeros(1,J-1);f=zeros(1,J);x=zeros(1,J);uu=zeros(1,J);forn=1:Mu(n,1)=0;u(n,M)=0;endforj=1:Ju(1,j)=sin(j*pi/J);endtt=1/M;xx=1/J;r=tt/(xx*xx);forn=1:Mforj=2:J-1;d(1,j)=(1-cita)*r*u(n,j+1)+(1-2*(1-cita)*r)*u(n,j)+(1-cita)*r*u(n,j-1);a=cita*r;b=1+2*cita*r;c=cita*r;ende(1,2)=c/b;f(1,2)=d(1,2)/b;fori=3:J-2e(1,i)=c/(b-e(1,i-1)*a);endfori=3:J-15f(1,i)=(d(1,i)+f(1,i-1)*a)/(b-e(1,i-1)*a);endx(1,J-1)=f(1,J-1);fori=J-2:-1:2x(1,i)=f(1,i)+e(1,i)*x(1,i+1);endfori=2:J-1u(n+1,i)=x(1,i);endendforj=1:J-1uu(1,j)=exp(-pi*pi)*sin(pi*j/J);h(j)=uu(1,j)-u(M,j);endxuhmesh(u);运行结果如下：cita=0M=10J=10x=1.0e+010*Columns1through90-2.49993.6685-3.25752.0465-0.93830.3102-0.07060.0099Column100u=1.0e+010*Columns1through960.00000.00000.00000.00000.00000.00000.00000.00000.000000.00000.00000.00000.00000.00000.00000.00000.00000-0.00000.00000.00000.00000.00000.00000.00000.000000.0000-0.00000.00000.00000.00000.00000.00000.00000-0.00000.0000-0.00000.00000.00000.00000.00000.000000.0000-0.00000.0000-0.00000.00000.0000-0.00000.00000-0.00000.0000-0.00000.0000-0.0000-0.00000.0000-0.000000.0001-0.00010.0001-0.00000.0000-0.0000-0.00000.00000-0.00230.0031-0.00240.0012-0.00040.0001-0.0000-0.000000.0754-0.10670.0887-0.05050.0200-0.00530.0009-0.00010-2.49993.6685-3.25752.0465-0.93830.3102-0.07060.0099Columns10through110.0000000000000000000000000h=1.0e+009*0.0000-0.75431.0667-0.88740.5047-0.20010.0534-0.00870.00077cita=0.5M=10J=10x=1.0e-003*Columns1through90-0.92150.20090.29070.0655-0.0568-0.0596-0.0252-0.0037Column1008u=Columns1through90.30900.58780.80900.95111.00000.95110.80900.58780.309000.13350.23380.29810.32500.31480.27040.19750.104100.01630.04250.06640.08130.08460.07600.05710.030500.01390.01740.01830.01850.01780.01560.01170.00630-0.00290.00120.00460.00610.00590.00490.00340.001700.00370.00170.00050.00040.00080.00100.00090.00060-0.0023-0.00020.00090.00090.00050.00020.0000-0.000000.00190.0001-0.0005-0.00030.00000.00020.00020.00010-0.00140.00010.00050.00020.0000-0.0001-0.0001-0.000000.0011-0.0002-0.0004-0.00010.00010.00010.00010.00000-0.00090.00020.00030.0001-0.0001-0.0001-0.0000-0.0000Columns10through110.0000000000000000000000000h=90.0000-0.00110.00020.00040.0002-0.0000-0.0000-0.0000-0.0000cita=1M=10J=10x=1.0e-003*Columns1through900.12940.24320.32770.37260.37260.32770.24320.1294Column10010u=Columns1through90.30900.58780.80900.95111.00000.95110.80900.58780.309000.18340.32640.42110.46280.45080.38870.28470.150300.07570.14060.18690.20980.20750.18080.13330.070600.03340.06250.08390.09500.09450.08280.06130.032500.01500.02820.03790.04300.04290.03770.02790.014900.00680.01270.01710.01950.01950.01710.01270.006800.00310.00580.00780.00880.00880.00780.00580.003100.00140.00260.00350.00400.00400.00350.00260.001400.00060.00120.00160.00180.00180.00160.00120.000600.00030.00050.00070.00080.00080.00070.00050.000300.00010.00020.00030.00040.00040.00030.00020.0001Columns10through110.000000000000000000000000011h=1.0e-003*0.0160-0.2551-0.4947-0.6737-0.7703-0.7729-0.6810-0.5061-0.26956、思考题1.什么是隐式格式?写出隐格式的截断误差.2.什么是-方法?写出-方法构造的差分格式及其稳定性条件?3.为提高-方法构造的差分格式截断误差的精度阶，应该如何选择参数?此时截断误差的主项为什么？解答如下：1、隐式格式：由下一行的值推出前一行的值截断误差：o（k^2+h^2）2、-方法：取【0,1】中的任何值时求差分方程的解12格式：0,0),()(2)1()(2002112111111nJnjjjnjnjnjnjnjnjnjnjUUxUxUUUxUUUtUU