椭圆方程的有限元法

椭圆形方程的有限元——上机实习报告数值线性代数上机实习报告第1页两点边值问题有限元法（必做）从Galerkin原理出发用线性元解两点边值问题：2,01(0)(1)0uuxxuu精确解：1221()[(23)(23)]21xxuxeeeexe。1.1变分形式从Galerkin原理出发推导出两点边值问题的变分形式，将积分区间等分为N份，则步长1,1,2,ihinN，记为h。写出有限元方程及系数矩阵元素。解：由题可知p(x)=1,q(x)=1,f(x)=2x所以有限元方程为))(,())(),((1xfuxxajiniji,j=1,2,...,n其中101])1(1[),(dhhajj，dhhajj101])1(1[),(，dhhdhhajj102102))1(1()1(),(，dhxhdhxhdxfjjbaj)1()()(2102101，计算有：6/556/256/7),(),(0),(),(),(0,),(33332133322322211211hhhuuuaaaaaaa6/556/256/7)31(216016)31(216016)31(2333321hhhuuuhhhhhhhhhhhhhh数值线性代数上机实习报告第2页1.2利用MATLAB求解问题的过程依次取2,2,3,4,5,6,7,8.nNn用MATLAB求解并图形比较数值解与精确解，用表格列出不同剖分时的2L误差。程序如下：function[u]=bianzhi(p,q,N)h=1/N;x=0:h:1;A=zeros(N-1);fori=2:N-1a3=@(t)-p./h+h.*q.*t.*(1-t);a2=@(t)p./h+h.*q.*(t.^2)+p./h+h.*q.*((1-t).^2);a1=@(t)-p./h+h.*q.*t.*(1-t);A(i,i-1)=quad(a1,0,1);A(i,i)=quad(a2,0,1);A(i-1,i)=quad(a3,0,1);endA(1,1)=quad(a2,0,1);f=zeros(N-1,1);fori=2:Nf1=@(t)(x(i-1)+h.*t).^2.*t+(x(i)+h.*t).^2.*(1-t);f(i-1)=h.*quad(f1,0,1);endu=inv(A)*f;precise_value=((exp(2)-1)^(-1)).*((2-3*exp(1))*exp(x)-(2*exp(1)-3)*exp(1-x))+x.^2+2;plot(x,[0;u;0],'b--',x,precise_value,'r--');legend('数值解','精确解');err=norm([0;u;0]-precise_value')end数值线性代数上机实习报告第3页N=4N=8N=16数值线性代数上机实习报告第4页N=32N=64N=128数值线性代数上机实习报告第5页N=256不同N对应的误差表格如下：N48163264128256err2.46e-048.56e-053.01e-051.06e-053.76e-061.33e-064.70e-071.3方法总结及分析由上图和表格可以看出，从Galerkin原理出发推导的两点边值问题的解和真实解的误差还是很小的，尤其，随着步长h的减小，所求的近似解与真实解更为接近。讨论组：庞瑞王丹刘锡兰李笑鹏