main.m

用matlab便的支持向量机预测方面的界面

%fxx program
clear; 
% DEFINE BOUNDARIES/PARAMETERS 
fid=fopen('input.data','r');
Lb=fscanf(fid,'%i',1);
D=fscanf(fid,'%i',1);
young=fscanf(fid,'%E',1);
nu=fscanf(fid,'%f',1);
fc=fscanf(fid,'%i',1);
ft=fscanf(fid,'%i',1);
P=0;
Dmat0 = (young/(1-nu^2))*[1 nu 0;nu 1 0;0 0 (1-nu)/2];

% SET UP NODAL COORDINATES
ndivl=fscanf(fid,'%i',1);
ndivw=fscanf(fid,'%i',1);
[x,conn1,conn2,numnod] = mesh2(Lb,D,ndivl,ndivw);
% DETERMINE DOMAINS OF INFLUENCE - UNIFORM NODAL SPACING
dmax=fscanf(fid,'%f',1);
xspac = Lb/ndivl;
yspac = D/ndivw;
dm(1,1:numnod)=dmax*xspac*ones(1,numnod);
dm(2,1:numnod)=dmax*yspac*ones(1,numnod);

% SET UP QUADRATURE 求积分 CELLS
ndivlq = 10;ndivwq = 4;
[xc,conn,numcell,numq] = mesh2(Lb,D,ndivlq,ndivwq);
stress(1:3,:)=0;
stress1(1:3,:)=0;
for ix=1:10
    P=P+1;
% SET UP GAUSS POINTS, WEIGHTS, AND JACOBIAN FOR EACH CELL
quado = 4;
[gauss] = gauss2(quado);
numq2 = numcell*quado^2;
gs = zeros(4,numq2);
[gs] = egauss(xc,conn,gauss,numcell);
L=length(gs);
stress(1:3,:)=stress;
gs(5,:)=stress(1,:);
gs(6,:)=stress(2,:);
gs(7,:)=stress(3,:);
% LOOP OVER GAUSS POINTS TO ASSEMBLE DISCRETE EQUATIONS
k = sparse(numnod*2,numnod*2);
for gg=gs
gpos = gg(1:2);
weight = gg(3);
jac = gg(4);
stressx(1,1)=gs(5);
stressx(2,1)=gs(6);
stressx(3,1)=gs(7);
[young,nu]=fxx(stressx,young,nu,fc,ft);
% PLANE STRESS DMATRIX 
Dmat = (young/(1-nu^2))*[1 nu 0;nu 1 0;0 0 (1-nu)/2];
% DETERMINE NODES IN NEIGHBORHOOD OF GAUSS POINT
v = domain(gpos,x,dm,numnod);
L = length(v);
en = zeros(1,2*L);
[phi,dphix,dphiy] = shape(gpos,dmax,x,v,dm);
Bmat=zeros(3,2*L);
for j=1:L
Bmat(1:3,(2*j-1):2*j) = [dphix(j) 0;0 dphiy(j);dphiy(j) dphix(j)];
end
for i=1:L
en(2*i-1) = 2*v(i)-1;
en(2*i) = 2*v(i);
end
k(en,en) = k(en,en)+sparse((weight*jac)*Bmat'*Dmat*Bmat);
end



% DETERMINE NODES ON BOUNDARY, SET UP BC'S
ind1 = 0;ind2 = 0;ind3=0;ind4=0;
for j=1:numnod
if(x(1,j)==0.0)
ind1=ind1+1;
nnu(1,ind1) = x(1,j);
nnu(2,ind1) = x(2,j);
no1(ind1)=j;
end
if(x(1,j)==Lb)
ind2=ind2+1;
nt(1,ind2) = x(1,j);
nt(2,ind2) = x(2,j);
no2(ind2)=j;
end
if(x(2,j)==D/2)
ind3=ind3+1;
nt2(1,ind3) = x(1,j);
nt2(2,ind3) = x(2,j);
no3(ind3)=j;
end
if(x(2,j)==-D/2)
ind4=ind4+1;
nt1(1,ind4) = x(1,j);
nt1(2,ind4) = x(2,j);
no4(ind4)=j;
end
end
lthu = length(nnu);
ltht = length(nt);
ubar = zeros(lthu*2,1);
f = zeros(numnod*2,1);

%SET UP GAUSS POINTS ALONG TRACTION BOUNDARY
ind=0;
gauss=gauss2(quado);
for i=1:(ltht-1)
ycen=(nt(2,i+1)+nt(2,i))/2;
jcob = abs((nt(2,i+1)-nt(2,i))/2);
for j=1:quado
mark(j) = ycen-gauss(1,j)*jcob;
ind = ind+1;
gst(1,ind)=nt(1,i);
gst(2,ind)=mark(j);
gst(3,ind)=gauss(2,j);
gst(4,ind)=jcob;
end
end


%SET UP GAUSS POINTS ALONG DISPLACEMENT BOUNDARY
gsu=gst;
gsu(1,1:ind)=zeros(1,ind);
qk = zeros(1,2*lthu);

%INTEGRATE FORCES ALONG BOUNDARY
Imo = (1/12)*D^3;
for gt=gst
gpos = gt(1:2);
weight=gt(3);
jac = gt(4);
v = domain(gpos,x,dm,numnod);
L = length(v);
en = zeros(1,2*L);
force=zeros(1,2*L);
[phi,dphix,dphiy] = shape(gpos,dmax,x,v,dm);
tx=0;
ty= -(P/(2*Imo))*((D^2)/4-gpos(2,1)^2);
for i=1:L
en(2*i-1) = 2*v(i)-1;
en(2*i) = 2*v(i);
force(2*i-1)=tx*phi(i);
force(2*i) = ty*phi(i);
end
f(en) = f(en) + jac*weight*force';
end

% INTEGRATE G MATRIX AND Q VECTOR ALONG DISPLACEMENT BOUNDARY
GG = zeros(numnod*2,lthu*2);
ind1=0;ind2=0;
for i=1:(lthu-1)
ind1=ind1+1;
m1 = ind1; m2 = m1+1;
y1 = nnu(2,m1); y2 = nnu(2,m2);
len = y1-y2;
for j=1:quado
ind2=ind2+1;
gpos = gsu(1:2,ind2);
weight = gsu(3,j);
jac = gsu(4,j);
fac11 = -(-P*nnu(2,m1))/(6*young*Imo);       %有p出现
fac2 = -P/(6*young*Imo);
xp1 = gpos(1,1);
yp1 = gpos(2,1);
uxex1 = (6*Lb-3*xp1)*xp1 + (2+nu)*(yp1^2 - (D/2)^2);
uxex1 = uxex1*fac11;
uyex1 = 3*nu*yp1^2*(Lb-xp1)+0.25*(4+5*nu)*xp1*D^2+(3*Lb-xp1)*xp1^2;
uyex1 = uyex1*fac2; 
N1 = (gpos(2,1)-y2)/len; N2 = 1-N1;
qk(2*m1-1) = qk(2*m1-1)-weight*jac*N1*uxex1;
qk(2*m1) = qk(2*m1) - weight*jac*N1*uyex1;
qk(2*m2-1) = qk(2*m2-1) -weight*jac*N2*uxex1;
qk(2*m2) = qk(2*m2) - weight*jac*N2*uyex1;
v = domain(gpos,x,dm,numnod);
[phi,dphix,dphiy] = shape(gpos,dmax,x,v,dm);
L = length(v);
for n=1:L
G1 = -weight*jac*phi(n)*[N1 0;0 N1];
G2 = -weight*jac*phi(n)*[N2 0;0 N2];
c1=2*v(n)-1;c2=2*v(n);c3=2*m1-1;c4=2*m1;
c5=2*m2-1;c6=2*m2;
GG(c1:c2,c3:c4)=GG(c1:c2,c3:c4)+ G1;
GG(c1:c2,c5:c6)=GG(c1:c2,c5:c6)+G2;
end
end
end

% ENFORCE BC'S USING LAGRANGE MULTIPLIERS
f = [f;zeros(lthu*2,1)];
f(numnod*2+1:numnod*2+2*lthu,1) = qk';
m = sparse([k GG;GG' zeros(lthu*2)]);
d=m\f;
u=d(1:2*numnod);
for i=1:numnod
u2(1,i) = u(2*i-1);
u2(2,i) = u(2*i);
end
% SOLVE FOR OUTPUT VARIABLES - DISPLACEMENTS
displ=zeros(1,2*numnod);
ind=0;
for gg=x
ind = ind+1;
gpos = gg(1:2);
v = domain(gpos,x,dm,numnod);
[phi,dphix,dphiy] = shape(gpos,dmax,x,v,dm);
displ(2*ind-1) = phi*u2(1,v)';
displ(2*ind) = phi*u2(2,v)';
end
for i=1:numnod
    disp2(1,i)=displ(2*i-1);
    disp2(2,i)=displ(2*i);
end
weiyi(1,ix)=P;
weiyi(2,ix)=100*abs(disp2(2,55));

% SOLVE FOR STRESSES AT GAUSS POINTS
ind = 0;
enorm=0;
L=length(gs);
stress(1:3,:)=stress;
gs(5,:)=stress(1,:);
gs(6,:)=stress(2,:);
gs(7,:)=stress(3,:);
for gg=gs
ind = ind+1;
gpos = gg(1:2);
weight = gg(3);
jac = gg(4);
stressx(1,1)=gs(5);
stressx(2,1)=gs(6);
stressx(3,1)=gs(7);
[young,nu]=fxx(stressx,young,nu,fc,ft);
% PLANE STRESS DMATRIX 
Dmat = (young/(1-nu^2))*[1 nu 0;nu 1 0;0 0 (1-nu)/2];
v = domain(gpos,x,dm,numnod);
L = length(v);
en = zeros(1,2*L);
[phi,dphix,dphiy] = shape(gpos,dmax,x,v,dm);
Bmat=zeros(3,2*L);
for j=1:L
Bmat(1:3,(2*j-1):2*j) = [dphix(j) 0;0 dphiy(j);dphiy(j) dphix(j)];
end
for i=1:L
en(2*i-1) = 2*v(i)-1;
en(2*i) = 2*v(i);
end
stress(1:3,ind) = Dmat*Bmat*u(en);
stressex(1,ind) = (1/Imo)*P*(Lb-gpos(1,1))*gpos(2,1);
stressex(2,ind) = 0;
stressex(3,ind) = -0.5*(P/Imo)*(0.25*D^2 - gpos(2,1)^2);
err = stress(1:3,ind)-stressex(1:3,ind);
err2 = weight*jac*(0.5*(inv(Dmat)*err)'*(err));
enorm = enorm + err2;
end
enorm = sqrt(enorm);
%求解节点应力
ind=0;
enorm=0;
sx=x;
L1=length(x);
stress1(1:3,:)=stress1;
sx(3,:)=stress1(1,:);
sx(4,:)=stress1(2,:);
sx(5,:)=stress1(3,:);
for gg=sx
ind=ind+1;
gpos=gg(1:2);
stressx(1,1)=sx(5);
stressx(2,1)=sx(6);
stressx(3,1)=sx(7);
[young,nu]=fxx(stressx,young,nu,fc,ft);
% PLANE STRESS DMATRIX 
Dmat = (young/(1-nu^2))*[1 nu 0;nu 1 0;0 0 (1-nu)/2];
v = domain(gpos,x,dm,numnod);
L=length(v);
en=zeros(1,2);
[phi,dphix,dphiy]=shape(gpos,dmax,x,v,dm);
Bmat=zeros(3,2);
for j=1:L
Bmat(1:3,(2*j-1):2*j)=[dphix(j) 0;0 dphiy(j);dphiy(j) dphix(j)];
en(2*j-1)=2*v(j)-1;
en(2*j)=2*v(j);
stress1(1:3,ind)=Dmat*Bmat*u(en);
end
stressex1(1,ind) = (1/Imo)*P*(Lb-gpos(1,1))*gpos(2,1);
stressex1(2,ind) = 0;
stressex1(3,ind) = -0.5*(P/Imo)*(0.25*D^2 - gpos(2,1)^2);
end
err1 = (stressex1(1:3,1:numnod)-stress1(1:3,1:numnod));
%for i=1:3
%    for j=1:numnod
%    wucha(i,j)=abs(err1(i,j))/abs((stressex(i,j)));
%end
%end
end
x;
stress1
stressex1
stressex;
disp2
figure(8);
for i=1:ix
weiyi1(1,i)=0+i*P/ix;
weiyi1(2,i)=i*weiyi(2,1);
end
plot(weiyi(2,:),weiyi(1,:),'-o',weiyi1(2,:),weiyi1(1,:),'-o');
xlabel('挠度');
ylabel('荷载');
title('荷载挠度图');
legend('非线性挠度','弹性挠度');
%wucha;
s1=x(1,1:numnod);
s2=x(2,1:numnod);
sx=x+40*disp2;
s3=sx(1,1:numnod);
s4=sx(2,1:numnod);
figure(1);
plot(s1,s2,'o',s3,s4,'*');
xlabel('悬臂梁各点X坐标');
ylabel('悬臂梁各点Y坐标');
title('悬臂梁各点位移图');
legend('变形前','变形后');
axis equal;
%plot the stress
ind=fscanf(fid,'%i',1);
if(ind==1)
figure(2);
t1=stress1(1,no1);
t11=stressex1(1,no1);
xt=x(2,no1);
plot(xt,t1,'-o');
xlabel('悬臂梁左端各点Y坐标');
ylabel('悬臂梁左端各点σx值');
title('悬臂梁左端各点计算解');
legend('EFG计算解');
figure(3);
t1=stress1(1,no2);
t11=stressex1(1,no2);
xt=x(2,no2);
plot(xt,t1,'-o');
xlabel('悬臂梁右端各点Y坐标');
ylabel('悬臂梁右端各点σx值');
title('悬臂梁右端各点计算解');
legend('EFG计算解');
figure(4);
t1=stress1(1,no3);
t11=stressex1(1,no3);
xt=x(1,no3);
plot(xt,t1,'-o');
xlabel('悬臂梁上端各点X坐标');
ylabel('悬臂梁上端各点σx值');
title('悬臂梁上端各点计算解');
legend('EFG计算解');
figure(5);
t1=stress1(1,no4);
t11=stressex1(1,no4);
xt=x(1,no4);
plot(xt,t1,'-o');
xlabel('悬臂梁下端各点X坐标');
ylabel('悬臂梁下端各点σx值');
title('悬臂梁下端各点计算解');
legend('EFG计算解');
%绘制等值线和应力云图
x1=linspace(0,Lb,ndivl+1);
x2=linspace(D/2,-D/2,ndivw+1);
z=stress1(1,1:numnod);
z1=reshape(z,ndivw+1,ndivl+1);
figure(6);
s=contour(x1,x2,z1,10);
clabel(s);
xlabel('悬臂梁各点X坐标');
ylabel('悬臂梁各点Y坐标');
title('悬臂梁各点σx等值线图');
axis equal;
figure(7)
pcolor(x1,x2,z1);
shading interp;
xlabel('悬臂梁各点X坐标');
ylabel('悬臂梁各点Y坐标');
title('悬臂梁各点σx应力云图');
colorbar;
axis equal;
grid on;
end
if(ind==2)
figure(2);
subplot(2,2,1);
t1=stress1(2,no1);
t11=stressex1(2,no1);
xt=x(2,no1);
plot(xt,t1,'-o',xt,t11,'-*');
xlabel('悬臂梁左端各点Y坐标');
ylabel('悬臂梁左端各点σy值');
title('悬臂梁左端各点计算解与理论解的拟合');
legend('EFG计算解','理论解');
subplot(2,2,2);
t1=stress1(2,no2);
t11=stressex1(2,no2);
xt=x(2,no2);
plot(xt,t1,'-o',xt,t11,'-*');
xlabel('悬臂梁右端各点Y坐标');
ylabel('悬臂梁右端各点σy值');
title('悬臂梁右端各点计算解与理论解的拟合');
legend('EFG计算解','理论解');
subplot(2,2,3);
t1=stress1(2,no3);
t11=stressex1(2,no3);
xt=x(1,no3);
plot(xt,t1,'-o',xt,t11,'-*');
xlabel('悬臂梁上端各点X坐标');
ylabel('悬臂梁上端各点σy值');
title('悬臂梁上端各点计算解与理论解的拟合');
legend('EFG计算解','理论解');
subplot(2,2,4);
t1=stress1(2,no4);
t11=stressex1(2,no4);
xt=x(1,no4);
plot(xt,t1,'-o',xt,t11,'-*');
xlabel('悬臂梁下端各点X坐标');
ylabel('悬臂梁下端各点σy值');
title('悬臂梁下端各点计算解与理论解的拟合');
legend('EFG计算解','理论解');
%绘制等值线和应力云图
x1=linspace(0,Lb,ndivl+1);
x2=linspace(D/2,-D/2,ndivw+1);
z=stress1(2,1:numnod);
z1=reshape(z,ndivw+1,ndivl+1);
figure(3);
s=contour(x1,x2,z1,10);
clabel(s);
xlabel('悬臂梁各点X坐标');
ylabel('悬臂梁各点Y坐标');
title('悬臂梁各点σy等值线图');
figure(4)
pcolor(x1,x2,z1);
shading interp;
xlabel('悬臂梁各点X坐标');
ylabel('悬臂梁各点Y坐标');
title('悬臂梁各点σy应力云图');
colorbar;
end
if(ind==3)
figure(2);
subplot(2,2,1);
t1=stress1(3,no1);
t11=stressex1(3,no1);
xt=x(2,no1);
plot(xt,t1,'-o',xt,t11,'-*');
xlabel('悬臂梁左端各点Y坐标');
ylabel('悬臂梁左端各点τxy值');
title('悬臂梁左端各点计算解与理论解的拟合');
legend('EFG计算解','理论解');
subplot(2,2,2);
t1=stress1(3,no2);
t11=stressex1(3,no2);
xt=x(2,no2);
plot(xt,t1,'-o',xt,t11,'-*');
xlabel('悬臂梁右端各点Y坐标');
ylabel('悬臂梁右端各点τxy值');
title('悬臂梁右端各点计算解与理论解的拟合');
legend('EFG计算解','理论解');
subplot(2,2,3);
t1=stress1(3,no3);
t11=stressex1(3,no3);
xt=x(1,no3);
plot(xt,t1,'-o',xt,t11,'-*');
xlabel('悬臂梁上端各点X坐标');
ylabel('悬臂梁上端各点τxy值');
title('悬臂梁上端各点计算解与理论解的拟合');
legend('EFG计算解','理论解');
subplot(2,2,4);
t1=stress1(3,no4);
t11=stressex1(3,no4);
xt=x(1,no4);
plot(xt,t1,'-o',xt,t11,'-*');
xlabel('悬臂梁下端各点X坐标');
ylabel('悬臂梁下端各点τxy值');
title('悬臂梁下端各点计算解与理论解的拟合');
legend('EFG计算解','理论解');
%绘制等值线和应力云图
x1=linspace(0,Lb,ndivl+1);
x2=linspace(D/2,-D/2,ndivw+1);
z=stress1(3,1:numnod);
z1=reshape(z,ndivw+1,ndivl+1);
figure(3);
s=contour(x1,x2,z1,10);
clabel(s);
xlabel('悬臂梁各点X坐标');
ylabel('悬臂梁各点Y坐标');
title('悬臂梁各点τxy等值线图');
figure(4)
pcolor(x1,x2,z1);
shading interp;
xlabel('悬臂梁各点X坐标');
ylabel('悬臂梁各点Y坐标');
title('悬臂梁各点τxy应力云图');
colorbar;
end
fclose(fid);