椭圆方程的有限元法.docVIP

椭圆形方程的有限元 ——上机实习报告 两点边值问题有限元法（必做） 从Galerkin原理出发用线性元解两点边值问题： 精确解：。 1.1变分形式 从Galerkin原理出发推导出两点边值问题的变分形式，将积分区间等分为N份，则步长，记为。写出有限元方程及系数矩阵元素。 解： 由题可知p(x)=1,q(x)=1,f(x)= 所以有限元方程为,j=1,2,...,n 其中， ， ， ， 计算有： 1.2利用MATLAB求解问题的 h=1/N; x=0:h:1; A=zeros(N-1); for i=2:N-1 a3=@(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); end A(1,1)=quad(a2,0,1); f=zeros(N-1,1); for i=2:N f1=@(t)(x(i-1)+h.*t).^2.*t+(x(i)+h.*t).^2.*(1-t); f(i-1)=h.*quad(f1,0,1); end u=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 N=4 N=8 N=16 N=32 N=64 N=128 N=256 不同N对应的误差表格如下： N 4 8 16 32 64 128 256 err 2.46e-04 8.56e-05 3.01e-05 1.06e-05 3.76e-06 1.33e-06 4.70e-07 1.3 方法总结及分析 由上图和表格可以看出，从Galerkin原理出发推导的两点边值问题的解和真实解的误差还是很小的，尤其，随着步长h的减小，所求的近似解与真实解更为接近。 讨论组：庞瑞 王丹 刘锡兰 李笑鹏 数值线性代数上机实习报告 第1页