边值问题的研究

西北农林科技大学理学院应用数学系微分方程数值解结课论文论文题目边值问题的研究2016年1月14日一问题重述对于下列边值问题01UAXB其中A为学号的倒数第2位，B为学号的倒数第1位。1差分截断误差、稳定性、收敛半径、递推（隐式）或方程组（显式）2有限元刚度矩阵、算法步骤及代码二问题分析题目明确指出使用差分方法和有限元解法。什么都不管先构造一种差分格式，然后对求解区域做划分将问题离散化，从微分方程的定解问题转化为求线性代数方程的解，以便于能够使用计算机进行计算。在这里选用的是中心差分法，同时将边界进行处理，同时用RITZ有限元法和GALERKIN法有限元法尝试去得到结果，最后再去比较两种解法所得到结果的精确性，分析相容性和截断误差等等。三解题过程1首先建立差分格式，考虑两点的边值问题，由题目知道1,097PQFX建立中心差分格式如下,1A,BDULPQFXABXU对求解区间做网格划分，在A到B之间取N1个节点，定义为XI（I取1到N）即将区间IA,B分为N个小区间1,2,IIIX由此得到区间的一个网格剖分。记。用表示网格内点，I1IIHHI1X2的集合，表示内点和界点的集合。1NXH0,NAB取相邻节点的中点，称为半整数点。由1,IX1/21,2IIIXXN节点又构成的一个对偶剖分。0/23/2INAB,A用差商代替微商，将方程（11）在内点离散化逼近边值问题（11）（12）I的差分方程为11/21/210131,2,IIIIHIIIIIINUXUXLUXPPQUXFHHH，当网格均匀，即时差分方程简化为,2I1/21/2/1/214HIIIIIIIIIILUPUPUQUF这相当于用一阶中心差商，二阶中心差商依次代替（11）的一阶微商和二阶微商的结果。这个方程就是中心差分格式。式（14）用方程组展开1/203/21/13/211/21/2/1/23/21/23/211/21KKKKKKNNNNPUPQUPFHHUFPUPQPFHHH这是一个以为未知量的线性方程组。121,，到此为止，中心差分格式展开完毕，接下来处理方程11将方程在节点离散化，由泰勒公式展开得24221121IIIIIDUXDUXUXXHOHH所以截断误差为4221IIDUXHROH下一步是分析差分格式的稳定性差分格式的截断误差，而边界条件的截断误差为2OHOH收敛性和稳定性是从不同角度讨论差分法的精确情况，稳定性主要是讨论初值的误差和计算中的舍入误差对计算结果的影响，收敛性则主要讨论推算公式引入的截断误差对计算结果的影响使用既收敛有稳定的差分格式才有比较可靠的计算结果，这也是讨论收敛性和稳定性的重要意义截断误差，即IHIIRULXU。241,1MHMMMHX差分方程组的解满足1AX,A,2,12MMMNUXBFN其中A、B代表边界点，代表边界点的取值。、42AX,2,196MBUXHU上式给出了差分方程的解的误差估计，而且表明当差分解收敛到原边值问0H题的解，收敛速度为。2H2接下来是有限元的解法从RITZ法出发，单元刚度矩阵为1,1,KIIIIIA12,11011,110,IIIIIIIIIIIIIIIIIIIAHPXHQXD按规则组装成总刚度矩阵。1NIIK令1,IIITFF其中101,IIIIIIIIFHFXDH以及12,TNBB1223,NFBFBF则有限元方程为KU从GALERKIN有限元法出发，GALERKIN有限元方程为1,12,NBIJIJAIAFDXN系数矩阵第J行只有三个非零元素，即11110,JJJJJJJJAHPXHQXD211110,JJJJJJJJJJJJJXX110,JJJJJJJJAHPHQD这里第一行只有两个非零元素第N行只有23JN2,A两个非零元素和1,1110,NNNNAXXH方程的右端项100,JJJJJJJFHFXDHFD四求解过程01UAXB其精确解为。算例中326ABUX1,0,PQX。FX1从RITZ法出发以将积分区间等分为10份为例，则步长，记为。1,2,IHNNH为KUB201083592012105U02821137460508以步长取为H1/10为例，从RITZ法出发的有限元法得到的数值解与精确解为RITZ数值解精确解0011540111352224521480320253094540790394404845046875549155316060095582056390061920662406421567025065000图像为0010203040506070809101234567XURITZ与与与与与与与与与与与与RITZ与与与与与与分析最大误差为0202500，RITZ有限元法求解两点边值问题很接近精确解。以步长取为H1/50为例，从RITZ法出发的有限元法得到的数值解与精确解图像为0010203040506070809101234567XURITZ与与与与与与与与与与与与RITZ与与与与与与最大误差为0002025步长1/101/501/1001/500RITZ最大误差0202500441000020250004491分析最大误差为002025，RITZ有限元法求解两点边值问题很接近精确解，且步长越大，误差越小。（2）从GALERKIN法出发以将积分区间等分为10份为例，则步长，记为。1,2,IHNNH为KUB2010798021165U0242131056GALERKIN有限元法最大误差0815000，图像为00102030405060708091012345678XUGALERKIN与与与与与与与与与与与与GALERKIN与与与与与与以将积分区间等分为100份为例，图像为0010203040506070809101234567XUGALERKIN与与与与与与与与与与与与GALERKIN与与与与与与分析最大误差为0080150，GALERKIN有限元法求解两点边值问题很接近精确解，且步长越大，误差越小。步长1/101/501/1001/500CALERKIN最大误差08015001600600801500016006最后收敛性和误差分析令和分别表示精确解和有限元解在剖分区间节点处的值，收敛性表示为UH01MAXKHIHIINMUUX记最大误差为ERR,则问题转化为在方程中，KER已知H和ERR，求解M和K的拟合问题。在MATLAB中拟合采用最小二乘法实现。对和进行最小二乘幂函数拟合，求得从RITZ法出发的误差阶为_ERRITZK09612,M04115对和进行最小二乘幂函数拟合，求得从GALERKIN法出发的误差_ERGALKINH阶为K1004,M3061五操作代码主程序FUNCTIONRITZA,B,ND2UAXB，U00，DU10A9B7A为学号倒数第二位，B为学号倒数第一位，N为剖分份分数调用格式RITZ9,7,10N为剖分份分数SYMSXUA0区间界点UB1区间界点U_A0左边界条件DU_B0右边界条件P1Q0FAXBF9X7H1/NX0H1KSTIFF_MATRIXP,Q,H,N,X得到总刚度矩阵KBVECTORP,Q,X,H,N,F得到BUKB差分解处理右边值条件U_B2HDU_BUEND14UEND/3U0UU_B精确解Y0DSOLVED2Y9X7,Y00,DY10,XPRECISE_VALUEDOUBLESUBSY0DETAUPRECISE_VALUEDETA_MAXMAXABSDETA最大误差FPRINTF最大的误差是FN,DETA_MAX差分解与精确解对比表FIGURESUBPLOT1,2,1PLOTX,U,B,X,PRECISE_VALUE,RXLABELXYLABELUTITLE差分解与精确解对比图LEGEND差分解,精确解,LOCATION,BEST误差图SUBPLOT1,2,2PLOTX,DETAXLABELXYLABELUEND子程序得到刚度矩阵KFUNCTIONKSTIFF_MATRIXP,Q,H,N,XSYMSXKZEROSN1,N1DIAG_0ZEROSN1,1确定K的对角元DIAG_1ZEROSN2,1确定K的偏离对角的上对角元DIAG_2ZEROSN2,1确定K的偏离对角的下对角元获取对角元素FORJ1N1DIAG_0JSUBSP,X,XJXJ1/2SUBSP,X,XJXJ1/2HSUBSQ,X,XJ1H2END获取A的第三条对角FORJ1N2DIAG_2JSUBSP,X,XJ1XJ2/2SUBS0,X,XJ2H/2END获取A的第二条对角FORJ1N2DIAG_1JSUBSP,X,XJ1XJ2/2SUBS0,X,XJ1H/2ENDKDIAGDIAG_0DIAGDIAG_1,1DIAGDIAG_2,1KN1,N12KN1,N22得到BFUNCTIONBVECTORP,Q,X,H,N,FBZEROSN1,1SYMSXX0H1FORJ1N1BJH2SUBSF,X,XJ1END系数矩阵B1B10SUBSP,X,X2X1/2SUBS0,X,X2H/2BN13BN120SUBSP,X,XNXN1/2SUBS0,X,XNH/2主程序P1Q0A9B7N100剖分份数H1/NX0H1A_GALERKINMATIRIXA,B,P,Q,H,NVALUES_F_GALERKINVECTOR1A,B,X,H,NU_GALERKINA_GALERKINVALUES_F_GALERKINY0DSOLVED2YAXB,Y00,DY10,XPRECISE_VALUEDOUBLESUBSY0GALERKIN有限元法解与精确解对比图FIGURESUBPLOT1,2,1PLOTX,0U_GALERKIN,B,X,PRECISE_VALUE,RXLABELXYLABELUTITLEGALERKIN有限元法解与精确解对比图LEGENDGALERKIN数值解,精确解,LOCATION,BESTERR_GALERKINMAXABS0U_GALERKINPRECISE_VALUEFPRINTFSPRINTFGALERKIN有限元法最大误差FN,ERR_GALERKIN子程序GALERKIN有限元法求解一维问题构造系数矩阵FUNCTIONA_GALERKINMATIRIXA,B,P,Q,H,N,A_GALERKINZEROSNFORI3NFUN1_GALERKINKSP/HHQKS1KSFUN2_GALERKINKSP/HHQKS2P/HHQ1KS2FUN3_GALERKINKSP/HHQKS1KSA_GALERKINI1,I2INTEGRALFUN1_GALERKIN,0,1A_GALERKINI1,I1INTEGRALFUN2_GALERKIN,0,1A_GALERKINI1,IINTEGRALFUN3_GALERKIN,0,1ENDFUN2_GALERKINKSP/HHQKS2P/HHQ1KS2A_GALERKIN1,1INTEGRALFUN2_GALERKIN,0,1FUN3_GALERKINKSP/HHQKS1KSA_GALERKIN1,2INTEGRALFUN3_GALERK