guaslinear.m

有限元分析作业源码, piecewise单元,积分 高斯近似算法两种

clear;

syms a x s

N=10; % number of elements
L=1/N;
n=2; % integral points
[y,W]=gauss(n);

K=zeros(N+1,N+1);
k=zeros(2,2);
f=zeros(N+1,1);
A=zeros(N+1,1);
X=zeros(N+1,1);


for i=1:N+1
    X(i)=(i-1)/N
end

fi(1)=(1-s)/2;
fi(2)=(1+s)/2;

for j=1:2
   for i=1:2
      for m=1:n
          s=y(m);
          a=eval(100*fi(i)*fi(j));
          k(i,j)=k(i,j)+W(m)*eval((4/L^2*diff(fi(i),s)*diff(fi(j),s)+eval(100*fi(i)*fi(j))))*L/2;
      end
   end
end

for i=1:N
    K(i,i)=K(i,i)+k(1,1);
    K(i,i+1)=k(1,2);
    K(i+1,i)=k(2,1);
    K(i+1,i+1)=k(2,2);
end
K(N+1,N+1)=K(N+1,N+1)+1;

for m=1:n
    s=y(m);
    f(1)=f(1)+W(m)*eval(100*fi(1)*L/2);          %f(1)
end

for i=2:N
    for m=1:n
        s=y(m);
        f(i)=f(i)+W(m)*eval(100*(fi(1)+fi(2))*L/2);  %f(2)-f(N)
    end
    
end

for m=1:n
    s=y(m);
    f(N+1)=f(N+1)+W(m)*eval(100*fi(2)*L/2);   % f(N+1)
end

KK=K; ff=f;

f(1)=5/3;
f(2)=f(2)-5/3*K(2,1);
f(N+1)=f(N+1)+1;

K(1,:)=0.0;
K(:,1)=0.0;
K(1,1)=1.0;

A=inv(K)*f;
PI=1/2*A'*KK*A-A'*ff


% =====Plot==========
X1=[0:0.001:1];
for n=1:1:1001
   T(n)=2/3*exp(-10*X1(n))+1; 
end

hold on
plot(X,A,'-.')
plot(X1,T)
legend('FEM-Gauss','Exact')';
hold off