piecewise.m

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

clear;

syms x
N=20;
L=1/N;
K=zeros(N+1,N+1);
f=zeros(N+1,1);
A=zeros(N+1,1);
fa=(1-x)/2;
fi=(1+x)/2;


K(1,1)=int((4/L^2*diff(fa,x)^2+100*fa*fa)*L/2,-1,1);
for i=2:(N+1)
    K(i-1,i)=int((4/L^2*diff(fa,x)*diff(fi,x)+100*fa*fi)*L/2,-1,1);
    K(i,i-1)=int((4/L^2*diff(fi,x)*diff(fa,x)+100*fi*fa)*L/2,-1,1);
    K(i,i)=int((4/L^2*diff(fi,x)^2+100*fi*fi)*L/2,-1,1)+K(1,1);
end
K(N+1,N+1)=int((4/L^2*diff(fi,x)^2+100*fi*fi)*L/2,-1,1)+1;

f(1)=int(100*fa*L/2,-1,1);
for i=2:N
    f(i)=int(100*(fa+fi)*L/2,-1,1);
end
f(N+1)=int(100*fi*L/2,-1,1)+1;
ff=f;

ff(1)=5.0/3.0;
ff(2)=f(2)-5.0/3.0*K(2,1);

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

A=inv(KK)*ff;

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

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

hold on
plot(X,A);
plot(X1,T);
hold off;

PI=1/2*A'*K*A-A'*f;
PIex=eval(1/2*int((diff(2/3*exp(-10*x)+1,x))^2+100*(2/3*exp(-10*x)+1)^2+1,0,1)-int(100*(2/3*exp(-10*x)+1)+1,0,1));