如何解决矩阵运算内存不足的问题？
本人在用无网格方法做粘弹性分析，由于利用弹性-粘弹性相应原理，在laplace域里进行求解，会产生字符变量，导致矩阵运算内存不足，请问各位大侠如何解决？跪求！clear;
% DEFINE BOUNDARIES/PARAMETERS
tic;
Lb = 48;
D = 12;
%young = 30e6;nu=0.3;
P=1000;
% PLANE STRESS DMATRIX
syms t Gt;
n1=2.824e12;
n2=5.01e11;
E1=4.701e10;
E2=2.309e10;
P1=n1/E1+(n1+n2)/E2;
P2=(n1*n2)/(E1*E2);
a=(P1+sqrt(P1^2-4*P2))/(2*P2);
b=(P1-sqrt(P1^2-4*P2))/(2*P2);
g1=E1/(2+2*n2);
g2=E1/(2+2*n1);
a1=(g2/n2-a)/(b-a);
a2=(g2/n2-b)/(a-b);
T1=1/a;
T2=1/b;
Gt=(g1+g2)*(a1*exp(-t/T1)+a2*exp(-t/T2));
Gs=laplace(Gt);
GG1=2*Gs;
syms K kk;
K=2.1e7*kk;
Ks=laplace(K);
GG2=3*Ks;
Ds=vpa((GG1/(2*GG1+GG2)*),5);
Dmat=Ds;
young=9*GG1*GG2/(3*GG2+GG1);
nu=(3*GG2-2*GG1)/(2*(3*GG2+GG1));
% SET UP NODAL COORDINATES
ndivl = 12;ndivw = 6;
= mesh2(Lb,D,ndivl,ndivw);
% DETERMINE DOMAINS OF INFLUENCE - UNIFORM NODAL SPACING
dmax=3.5;
xspac = Lb/ndivl;
yspac = D/ndivw;
dm = zeros(2,numnod);
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;
= mesh2(Lb,D,ndivlq,ndivwq);
% SET UP GAUSS POINTS, WEIGHTS, AND JACOBIAN FOR EACH CELL
quado = 4;
= pgauss(quado);
numq2 = numcell*quado^2;
gs = zeros(4,numq2);
= egauss(xc,conn,gauss,numcell);
% LOOP OVER GAUSS POINTS TO ASSEMBLE DISCRETE EQUATIONS
k =vpa(sym(zeros(numnod*2,numnod*2)),5);
for gg=gs
gpos = gg(1:2);
weight = gg(3);
jac = gg(4);
% DETERMINE NODES IN NEIGHBORHOOD OF GAUSS POINT
v = domain(gpos,x,dm,numnod);
L = length(v);
en = zeros(1,2*L);
= shape(gpos,dmax,x,v,dm);
Bmat=zeros(3,2*L);
for j=1:L
Bmat(1:3,(2*j-1):2*j) = ;
end
for i=1:L
en(2*i-1) = 2*v(i)-1;
en(2*i) = 2*v(i);
end
A1=vpa((sym(k(en,en))),5);
%A=str2mat((weight*jac)*Bmat'*Dmat*Bmat);
A=vpa(((weight*jac)*Bmat'*Dmat*Bmat),5);
% for h=1,length(en)
% for n=1,length(en)
%for g=1,32;
% for e=1,32;
% A11=A(g,e);
%k(en(h),en(n)) = k(en(h),en(n))+A11;
% end
% end
%end
% end
k(en,en) = A1+A;
end
k=vpa(k,5);
toc;
tt=toc
% DETERMINE NODES ON BOUNDARY, SET UP BC'S
ind1 = 0;ind2 = 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);
end
if(x(1,j)==Lb)
ind2=ind2+1;
nt(1,ind2) = x(1,j);
nt(2,ind2) = x(2,j);
end
end
lthu = length(nnu);
ltht = length(nt);
ubar = zeros(lthu*2,1);
f = vpa(zeros(numnod*2,1),5);
%SET UP GAUSS POINTS ALONG TRACTION BOUNDARY
ind=0;
gauss=pgauss(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 = vpa(zeros(1,2*lthu),5);
%INTEGRATE FORCES ALONG BOUNDAR
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);
= 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);
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(1,2*m1-1) = qk(1,2*m1-1)-weight*jac*N1*uxex1;
qk(1,2*m1) = qk(1,2*m1) - weight*jac*N1*uyex1;
qk(1,2*m2-1) = qk(1,2*m2-1) -weight*jac*N2*uxex1;
qk(1,2*m2) = qk(1,2*m2) - weight*jac*N2*uyex1;
v = domain(gpos,x,dm,numnod);
= shape(gpos,dmax,x,v,dm);
L = length(v);
for n=1:L
G1 = -weight*jac*phi(n)*;
G2 = -weight*jac*phi(n)*;
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(numnod*2+1:numnod*2+2*lthu,1) = -qk';
m =;
d=vpa(m\f,10);
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);
= shape(gpos,dmax,x,v,dm);
displ(2*ind-1) = phi*u2(1,v)';
displ(2*ind) = phi*u2(2,v)';
end
%SOLVE FOR STRESSES AT GAUSS POINTS
ind = 0;
enorm = 0;
for gg=gs
ind = ind+1;
gpos = gg(1:2);
weight = gg(3);
jac = gg(4);
v = domain(gpos,x,dm,numnod);
L = length(v);
en = zeros(1,2*L);
= shape(gpos,dmax,x,v,dm);
Bmat=zeros(3,2*L);
for j=1:L
Bmat(1:3,(2*j-1):2*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);
stres**(1,ind) = (1/Imo)*P*(Lb-gpos(1,1))*gpos(2,1);
stres**(2,ind) = 0;
stres**(3,ind) = -0.5*(P/Imo)*(0.25*D^2 - gpos(2,1)^2);
err = ilaplace(stress(1:3,ind)-stres**(1:3,ind));
err2 = weight*jac*(0.5*(inv(Dmat)*err)'*(err));
enorm = enorm + err2;
end
enorm = sqrt(enorm)
换个64位系统试一试 qibbxxt 发表于 2012-6-11 20:11 static/image/common/back.gif
换个64位系统试一试
就是在64位的服务器上计算的，还是不行。 ccz19860422 发表于 2012-6-21 10:21 static/image/common/back.gif
就是在64位的服务器上计算的，还是不行。
正常来说，应该是i的程序的问题 qibbxxt 发表于 2012-6-25 14:56 static/image/common/back.gif
正常来说，应该是i的程序的问题
能不能请高手指点一二啊:P 本帖最后由 pasuka 于 2012-6-28 21:04 编辑
ccz19860422 发表于 2012-6-28 09:09 static/image/common/back.gif
能不能请高手指点一二啊
符号运算的表达式推导出来后，简化后对照着再编程用matlab计算，不要直接vpa或者subs
如果这样还不行的话，根据符号推导的表达式用C、C++、Fortran编程计算
如果符号运算就报内存不足的话，换mathematics推导公式，其余同上 十分感谢！
页:
[1]