一维有限元法..doc

实习三、一维问题的有限元方法 一）实习问题： 设 将原问题的边界条件齐次化 二）算法描述： 1，单元剖分 2，i=1 3,计算数值积分：即得单元上的 4，将迭加到总的中 5，若i=n,则i=i+1并转到底三步；否则继续下一步 6，根据边界条件调整（掐头去尾），即得 A 和b 7，解线性方程组Au=b,得u从而的 三）matlab程序： 问题M文件: function [y1]=f1(x,h) y1=h*(1/6*(-24*exp(x+h)*h-3*h^2*exp(1)*x+48*exp(x+h)-24*exp(x+h)*x-h^3*exp(1)+3*h^2*exp(-1)*x+h^3*exp(-1)+24*exp(x)*h*x-48*exp(x)-24*exp(x)*h+24*x*exp(x))/h^2); function [y2]=f2(x,h) y2=h*(1/6*(48*exp(x+h)*h-24*exp(x+h)*h^2-24*exp(x+h)*x*h-3*h^2*exp(1)*x+3*h^2*exp(-1)*x-2*h^3*exp(1)+2*h^3*exp(-1)+24*exp(x+h)*x-48*exp(x+h)-24*x*exp(x)+48*exp(x))/h^2); 主M 文件： function []=yiweiyouxianyuan(n) %对x的随机剖分及区间长度的计算 l=abs(rand(1,n-1)); for i=1:n-1 for j=(i+1):n-1 if l(i)l(j) l2=l(i); l(i)=l(j); l(j)=l2; end end end for i=1:n-1 x(i+1)=l(i); end x(1)=0;x(n+1)=1; for i=1:n h(i)=x(i+1)-x(i); end %一次区间元法 %A的求解 for i=1:n a(i,1)=1/h(i)+h(i)/3; a(i,2)=-1/h(i)+h(i)/6; a(i,3)=a(i,2); a(i,4)=a(i,1); end for i=1:n-1 A(i,i)=a(i,4)+a(i+1,1); end for i=1:n-2 A(i,i+1)=a(i+1,2); end for i=2:n-1 A(i,i-1)=a(i,3); end %b的求解 for i=2:n+1 b1(i,i-1)=f1(x(i-1),h(i-1)); b1(i,i)=f2(x(i-1),h(i-1)); end for i=2:n b2(i)=b1(i,i)+b1(i+1,i); end for i=1:n-1 b(i)=b2(i+1); end u=inv(A)*b ; for i=2:n un(i)=u(i-1); end un(1)=0;un(n+1)=0; %还原原始的u值 uz=un+x*(exp(1)-exp(-1)); uz %真解的求解 for i=1:n+1 u1(i)=exp(x(i))-exp(-x(i))+x(i)*x(i)*exp(x(i))-x(i)*exp(x(i)); end u1 %误差的计算 for i=1:n+1 e(i)=abs((uz(i)-u1(i))/u1(i)*100); end e %作图 subplot(1,2,1) plot(x,u1) xlabel(自变量x的范围);ylabel(函数值u的取值);title(真解的图象);grid subplot(1,2,2) plot(x,uz) xlabel(自变量x的范围);ylabel(函数值u的取值);title(有限元法算得的近似解的图象);grid 四）图形显示的计算结果： 将区间随机分（利用rand函数）为20份的计算结果： 图1：预测值与真实值的作图比较 表1：预测值与真实值的数值比较 真实值 预测值 误差（%） 0 0 0.0000 0.0135 0.0135 0.0726 0.0336 0.0336 0.0725 0.1816 0.1815 0.0657 0.3492 0.3490 0.0515 0.4117 0.4115 0.0465 0.5460 0.5458 0.0377 0.6646 0.66