边值问题的研究

西北农林科技大学理学院应用数学系《微分方程数值解》结课论文论文题目边值问题的研究2016年1月14日一·问题重述对于下列边值问题：*(0)(1)0uAxBuu其中A为学号的倒数第2位，B为学号的倒数第1位。(1)差分：截断误差、稳定性、收敛半径、递推（隐式）或方程组（显式）(2)有限元：刚度矩阵、算法步骤及代码二·问题分析题目明确指出使用差分方法和有限元解法。什么都不管先构造一种差分格式，然后对求解区域做划分将问题离散化，从微分方程的定解问题转化为求线性代数方程的解，以便于能够使用计算机进行计算。在这里选用的是中心差分法，同时将边界进行处理，同时用Ritz有限元法和Galerkin法有限元法尝试去得到结果，最后再去比较两种解法所得到结果的精确性，分析相容性和截断误差等等。三·解题过程1·首先建立差分格式，考虑两点的边值问题，由题目知道1,0,9*7pqfxx建立中心差分格式如下()(),1.1(a),(b)dduLupqufxaxbdxdxuu对求解区间做网格划分，在a到b之间取N+1个节点，定义为xi（i取1到N）即将区间I=[a,b]分为N个小区间1:,1,2,.iiiIxxxiN由此得到区间I的一个网格剖分。记1iiihxx。用hI表示网格内点1x，2x，…，1Nx的集合，hI表示内点和界点0,Nxaxb的集合。取相邻节点1,iixx的中点1/211()(1,2,,)2iiixxxiN，称为半整数点。由节点01/23/2-1/21/2iNNaxxxxxxb又构成[,]ab的一个对偶剖分。用差商代替微商，将方程（1.1）在内点ix离散化.逼近边值问题（1.1）（1.2）的差分方程为：111/21/2110()()()()2()[]()1.31,2,,1,iiiihiiiiiiiiiiNuxuxuxuxLuxppquxfhhhhiNuu，，当网格均匀，即(1,2,,)ihhiN时差分方程简化为1/211/21/21/2121[()]1.4hiiiiiiiiiiiLupuppupuqufh这相当于用一阶中心差商，二阶中心差商依次代替（1.1）的一阶微商和二阶微商的结果。这个方程就是中心差分格式。式（1.4）用方程组展开：1/203/21/2113/2212221/211/21/21/212223/221/23/2111/21222111[()]111[()]111[()]kkkkkkkkkNNNNNNNNNpuppqupufhhhpuppqupufhhhpuppqupufhhh这是一个以121,,Nuuu……，为未知量的线性方程组。到此为止，中心差分格式展开完毕，接下来处理方程(1.1)将方程1.1在节点离散化，由泰勒公式展开得：242211224()2()()12iiiiiduxduxuxuxuxhOhhdxdx所以截断误差为422412iiduxhRuOhdx下一步是分析差分格式的稳定性差分格式的截断误差：2()oh，而边界条件的截断误差为()oh收敛性和稳定性是从不同角度讨论差分法的精确情况，稳定性主要是讨论初值的误差和计算中的舍入误差对计算结果的影响，收敛性则主要讨论推算公式引入的截断误差对计算结果的影响.使用既收敛有稳定的差分格式才有比较可靠的计算结果，这也是讨论收敛性和稳定性的重要意义.截断误差：()()[]ihiiRuLuxLu，即2411(()),,12mhmmmmmmhRLuxuuxx。差分方程组的解mu满足：111max,()()max,1,2,,12mmmmmNuxabxfmN其中a、b代表边界点，、代表边界点的取值。42max,1,2,,196mmaxbbauxuhuxmN上式给出了差分方程的解的误差估计，而且表明当0h差分解收敛到原边值问题的解，收敛速度为2h。2·接下来是有限元的解法从Ritz法出发，单元刚度矩阵为：()()1,11,()()(),1Kiiiiiiiiiiiiiaaaa1()121,11101()121101()()11,,1110[()()(1)],[()()],[()()(1)].iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiahpxhhqxhdahpxhhqxhdaahpxhhqxhd按规则组装成总刚度矩阵()1niiKK。令()()()1(,),iiiTiifff其中1()1101()10()(1),(),iiiiiiiiiifhfxhdfhfxhd以及12(,,,),Tnbbbb(1)(2)(2)(3)()111222,,,,nnnbffbffbf则有限元方程为.Kub从Galerkin有限元法出发，Galerkin有限元方程为：1(,),1,2,,.nbijijaiaufdxjn系数矩阵第j行只有三个非零元素，即111110(,)[()()(1)],jjjjjjjjahpxhhqxhd11211011211110(,)[()()][()()(1)],jjjjjjjjjjjjjjahpxhhqxhdhpxhhqxhd11111110(,)[()()(1)],jjjjjjjjahpxhhqxhd这里2,3,,1.jn第一行只有两个非零元素：1112(,),(,).aa第n行只有两个非零元素：1(,)nna和112110(,)[()()]nnnnnnnnahpxhhqxhd方程的右端项1111100(,)()()(1)jjjjjjjfhfxhdhfxhd四·求解过程*(0)(1)0uAxBuu其精确解为32()622ABAuxxBx。算例中()1,()0,pxqx()*fxAxB。(1)从Ritz法出发以将积分区间等分为10份为例，则步长1,1,2,ihinN，记为h。.Kub为：20100.83501020100.92501020101.01501020101.19501020101.1050*u1020101.28501020101.37501020101.46501020101.555010100.7850以步长取为h=1/10为例，从Ritz法出发的有限元法得到的数值解与精确解为Ritz数值解精确解001.15401.11352.22452.14803.20253.09454.07903.94404.84504.68755.49155.31606.00955.82056.39006.19206.62406.42156.702506.5000图像为分析：最大误差为0.202500，Ritz有限元法求解两点边值问题很接近精确解。00.10.20.30.40.50.60.70.80.9101234567xuRitz有限元法解与精确解对比图Ritz数值解精确解以步长取为h=1/50为例，从Ritz法出发的有限元法得到的数值解与精确解图像为最大误差为0.002025步长1/101/501/1001/500Ritz最大误差0.20250.0441000.020250.004491分析：最大误差为0.02025，Ritz有限元法求解两点边值问题很接近精确解，且步长越大，误差越小。（2）从Galerkin法出发以将积分区间等分为10份为例，则步长1,1,2,ihinN，记为h。.Kub为：00.10.20.30.40.50.60.70.80.9101234567xuRitz有限元法解与精确解对比图Ritz数值解精确解20100.791020100.881020100.971020101.061020101.15*u1020101.241020101.331020101.421020101.5110101.6Galerkin有限元法最大误差:0.815000，图像为：以将积分区间等分为100份为例，图像为：00.10.20.30.40.50.60.70.80.91012345678xuGalerkin有限元法解与精确解对比图Galerkin数值解精确解分析：最大误差为0.080150，Galerkin有限元法求解两点边值问题很接近精确解，且步长越大，误差越小。步长1/101/501/1001/500Calerkin最大误差0.801500.160060..0801500.016006最后收敛性和误差分析：令u和hu分别表示精确解和有限元解在剖分区间节点处的值，收敛性表示为01max()()khihiiNMhuuuxux记最大误差为err,则问题转化为在方程中，kMherr已知h和err，求解M和k的拟合问题。在Matlab中拟合采用最小二乘法实现。对_errRitz和h进行最小二乘幂函数拟合，求得从Ritz法出发的误差阶为k=0.9612,M=0.4115.00.10.20.30.40.50.60.70.80.9101234567xuGalerkin有限元法解与精确解对比图Galerkin数值解精确解对_errGalerkin和h进行最小二乘幂函数拟合，求得从Galerkin法出发的误差阶为k=1.004,M=3.061.五·操作代码主程序：functionRitz(a,b,N)%-D2u=a*x+b，u(0)=0，Du(1)=0%a=9;b=7;%a为学号倒数第二位，b为学号倒数第一位，%N为剖分份分数%调用格式：Ritz(9,7,10);%N为剖分份分数symsx;ua=0;%区间界点ub=1;%区间界点u_a=0;%左边界条件du_b=0;%右边界条件p=1;q=0;f=a*x+b;%f=9x+7h=1/N;X=0:h:1;K=Stiff_matrix(p,q,h,N,X);%得到总刚度矩阵KB=vector(p,q,X,h,N,f);%得到BU=K\B;%差分解%处理右边值条件u_b=(2*h*du_b-U(end-1)+4*U(end))/3;U=[0;U;u_b];%精确解y0=dsolve('-D2y=9*X+7','y(0)=0','Dy(1)=0','X');precise_value=double(subs(y0));deta=U-precise_value';deta_max=max(abs(deta));%最大误差fprintf('最大的误差是%f\n',deta_max)%差分解与精确解对比表figure;subplot(1,2,1);plot(X,U,'b*',X,precise_value,'r--');xlabel('x');