内容发布更新时间 : 2024/12/27 10:28:11星期一 下面是文章的全部内容请认真阅读。
有限元法课程作业(二):一维问题的有限元方法
一、解题步骤:
’-u’?u’?-4xex?x(e?e?1) 将原问题的边界齐次化:
u(0)?0, u(1)?0 1)单元剖分:ei(i?1,2,?400,n); 2) i=1 A?0b?0;
i)(i)(i)(i)(i)(i)ii 3)计算数值积分:ai(?1,i?1,ai?1,i,ai,i?1,ai,i,bi?1,bi即得单元上的Ab;
~
~
~~ 4)将Ab迭加到总的Ab中;
5)若i<=n,则i=i+1并转到底三步;否则继续下一步; 6)根据边界条件调整Ab(掐头去尾),即得 A 和b; 7)解线性方程组Au=b,得u从而的uh。 二 、程序编写
主程序:
function []=main(n)
n=400;%对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;
~
~
ii 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
两个子程序编写: 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);
三、 MATLAB程序运行结果
图1:理论值与有限元法近似值图
表1:理论值与有限元法近似值表
理论值 0 0.0144 0.0172 0.0303 0.0355 0.0508 0.0523 0.0561 0.0575 0.0592 0.0593 0.0664 0.0684 0.0739 0.0859 0.0952 0.1091 0.1152 有限元法近似值 0 0.0144 0.0172 0.0303 0.0355 0.0508 0.0523 0.0561 0.0575 0.0592 0.0593 0.0664 0.0685 0.0739 0.0859 0.0952 0.1091 0.1152 误差(100%) 0 0 0 0 0 0 0 0 0 0 0 0 0.1462 0 0 0 0 0 理论值 0.668 0.685 0.7045 0.7074 0.7136 0.7296 0.7908 0.7934 0.8357 0.8433 0.8436 0.899 0.906 0.9078 0.9251 0.959 0.961 0.9628 有限元法近似值 0.668 0.685 0.7045 0.7074 0.7136 0.7297 0.7908 0.7935 0.8357 0.8433 0.8436 0.899 0.906 0.9079 0.9251 0.959 0.961 0.9629 误差(100%) 0 0 0 0 0 0.0137 0 0.0126 0 0 0 0 0 0.0110 0 0 0 0.0104