偏微分方程数值解法13

1/16《偏微分方程数值解法》13——差分格式的矩阵分析变系数方程差分格式有限体积法简介差分格式的矩阵MATLAB创建矩阵2/16变系数两点边值问题)(,)0(0),(LuuLxxfcuubujjjjxjjxfucuhbuh2122二阶差分格式jjjjjjjjjfucuuhbuuuh)(2)2(111112jjch22记2/1jjhb2/1jjhbjjfhF2三对角矩阵A=TriDiag{}jjj,,AU=F3/16四阶差分逼近)(121121214222222hOLhLbhbhxxx四阶差分格式jjxjxjxuLhLubhbh)121()12121(222222)(xfcuubuxxbDDL2jjjjfucLu1)3()2(1)1(11jjjjjjjjjjjjfrfrfruuu4/16),(),(]),([]),([yxfuyxcyuyxbyxuyxax差分格式ijijijijyijyijxijxfucubhuah][1][122其中jijiijxuuu,2/1,2/12/1,2/1,jijiijyuuujijiijjijijijiuauaaua,1,2/1,2/1,2/1,1,2/1)()]([][,2/1,2/1jijiijxijxijxuuaua)()(,1,2/1,1,2/1jiijjiijjijiuuauua变系数二维问题5/16同理1,2/1,2/1,2/1,1,2/1,)(][jijiijjijijijiijyijyububbubub计算模板41032其中2/1,2jiajia,2/13jia,2/112/1,4jiaijjijijijicaaaa2/1,2/1,,2/1,2/10ijfh2练习:设2)1,0(),(yx写出差分格式对应的矩阵6/16变量代换方法引入新函数xuvyuw),(),(]),([]),([yxfuyxcyuyxbyxuyxax),(),(]),([]),([yxfuyxcwyxbyvyxax由于)(6121432hODhDhxxx)(12222hODhxx)()611(422hODDhxx)()611(42hODxx7/16原微分方程变换后的差分形式为][21]4[61,1,1,1,1jijijiijjiuuhvvv][21]4[611,1,1,1,jijijiijjiuuhvww])[611(21])[611(211,1,1,1,2,1,1,1,12jijijijixjijijijiywbwbhvavahijyxijijyxfuc)611)(611())(611)(611(2222补充方程如下)4(61)611(,1,12jiijjiijxuuuu8/16有限体积法在求解区域上选择的一个局部区域进行积分.以内网格点(xi，yj)为中心的局部积分区域},|),{(21,21,,21,21jijijijiijyyyxxxyx2/12/12/12/1]),([]),([jjiiijyyxxDdxxuyxaxdydxdyxuyxax2/12/1]),(),([2/12/1jjyyiidyxuyxaxuyxa利用数值积分的中矩形公式，上式近似等于9/16)],(),(),(),([2/12/12/12/1jijijijiyxxuyxayxxuyxah利用差商近似一阶导数，有近似表达式)()(]),([,1,2/1,1,2/1jiijjiijjijiDuuauuadxdyxuyxaxij)()(]),([1,2/1,2/1,2/1,jiijjiijjijiDuubuubdxdyyuyxbyij同理根据数值求积公式ijDfhdxdyyxfij2),(结果与直接差分格式同10/16特殊三对角矩阵0110110S对称三对角矩阵TSI二阶差分算子2112112TSI22x01101100T11/16x01101100Tjjjjxjjxfucuhbuh2122两点边值问题二阶差分格式系数矩阵IDhTDhTAcb202其中Db=Diag{b1,······,bN-1}Dc=Diag{c1,······,cN-1}12/16二维泊松方程ijijyxfuh)(1222系数矩阵)(ITTIA五点差分格式),()(2222yxfyuxu二维Helmholtz方程),(2yxfukuijijijyxfukuh2222)(1系数矩阵IIkhITTIA22五点差分格式13/16二维对流扩散方程),(yxfruqupuuuyxyyxx),(yxfruuPuP=[p,q]ijijyxyxfurqph])([1222系数矩阵IrIIqTTpIITTIA00)(差分格式MATLAB的张量积运算kron(X,Y)istheKroneckertensorproductofXandY14/16三对角矩阵创建方法n=4;E=eye(n)；sub=diag(ones(n-1,1),1)T0=sub-sub‘T=-2*E+sub+sub'T0=0100-10100-10100-10T=-21001-21001-21001-2sub=010000100001000015/16)(ITTIA矩阵创建方法n=3;E=eye(n);sub=diag(ones(n-1,1),1);T=-2*E+sub+sub‘;A=-(kron(E,T)+kron(T,E))A=4-10-100000-14-10-100000-1400-1000-1004-10-1000-10-14-10-1000-10-1400-1000-1004-100000-10-14-100000-10-14spy(A)16/16IrIIqTTpIITTIA00)(创建矩阵n=4;E=eye(n);sub=diag(ones(n-1,1),1);T0=sub-sub;T=-2*E+sub+sub';A=-(kron(E,T)+kron(T,E));p=1;q=1;r=1;A=A+p*kron(E,T0)+q*kron(T0,E);A=A+r*kron(E,E);Spy(A)