一维有限元法

实习三、一维问题的有限元方法一）实习问题：设''14(0,1)(0)0,(1)xuuxexuuee，~1()uxeeu令将原问题的边界条件齐次化''~~1~~4()(0,1)(0)0,(1)0xxexeexuuuu，二）算法描述：1()21,101101()1,011101(),101110(),111[()()()]1[()()()()]1[()()()()]1[()(iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiaphhqhNdhaphhqhNNdhaphhqhNNdhaphhqhhxxxxxxxx12010)()()]NNd1()21,101101[()()()]iiiiiiiiiaphhqhNdhxx1()10101()110()()()()iiiiiiiiiibhfhNdbhfhNdxx1，单元剖分(1,2,,)iineL2，i=1~~00Ab3,计算数值积分：()()()()()()1,11,,1,1,,,,,iiiiiiiiiiiiiiiiaaaabb即得单元上的iiAb4，将iiAb迭加到总的~~Ab中5，若i=n,则i=i+1并转到底三步；否则继续下一步6，根据边界条件调整~~Ab（掐头去尾），即得A和b7，解线性方程组Au=b,得u从而的hu三）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));fori=1:n-1forj=(i+1):n-1ifl(i)l(j)l2=l(i);l(i)=l(j);l(j)=l2;endendendfori=1:n-1x(i+1)=l(i);endx(1)=0;x(n+1)=1;fori=1:nh(i)=x(i+1)-x(i);end%一次区间元法%A的求解fori=1:na(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);endfori=1:n-1A(i,i)=a(i,4)+a(i+1,1);endfori=1:n-2A(i,i+1)=a(i+1,2);endfori=2:n-1A(i,i-1)=a(i,3);end%b的求解fori=2:n+1b1(i,i-1)=f1(x(i-1),h(i-1));b1(i,i)=f2(x(i-1),h(i-1));endfori=2:nb2(i)=b1(i,i)+b1(i+1,i);endfori=1:n-1b(i)=b2(i+1);endu=inv(A)*b';fori=2:nun(i)=u(i-1);endun(1)=0;un(n+1)=0;%还原原始的u值uz=un'+x'*(exp(1)-exp(-1));uz%真解的求解fori=1:n+1u1(i)=exp(x(i))-exp(-x(i))+x(i)*x(i)*exp(x(i))-x(i)*exp(x(i));endu1'%误差的计算fori=1:n+1e(i)=abs((uz(i)-u1(i))/u1(i)*100);ende'%作图subplot(1,2,1)plot(x,u1)xlabel('自变量x的范围');ylabel('函数值u的取值');title('真解的图象');gridsubplot(1,2,2)plot(x,uz)xlabel('自变量x的范围');ylabel('函数值u的取值');title('有限元法算得的近似解的图象');grid四）图形显示的计算结果：将区间随机分（利用rand函数）为20份的计算结果：00.5100.511.522.5自变量x的范围函数值u的取值真解的图象00.5100.511.522.5自变量x的范围函数值u的取值有限元法算得的近似解的图象图1：预测值与真实值的作图比较表1：预测值与真实值的数值比较真实值预测值误差（%）000.00000.01350.01350.07260.03360.03360.07250.18160.18150.06570.34920.34900.05150.41170.41150.04650.54600.54580.03770.66460.66440.03160.73150.73130.02870.82680.82660.02540.87530.87510.02400.87820.87800.02391.09411.09390.01811.32561.32540.01281.35511.35490.01221.52691.52680.00921.82531.82520.00501.99631.99630.00302.02042.02030.00272.15292.15290.00152.35042.35040.0000分析：从图1与表1中可以看出，预测值与真实值极其的接近，其误差最大值是0.0726%，所以，该有限元方法对此题的近似进度相当的高，可以作为求值的近似。