切换到宽版

 找回密码
 注册
Simdroid-非首页
Simwe仿真论坛(forum.simwe.com)，CAE/CAD/CAM/，FEA/FEM/有限元分析论坛---(手机验证注册)»论坛 › ..:: 科学软件专区 › K01：MatLab › 如何解决矩阵运算内存不足的问题？
返回列表 发新帖
查看: 92|回复: 6

如何解决矩阵运算内存不足的问题？

[复制链接]
ccz19860422
发表于 2012-6-8 09:23:09 | 显示全部楼层 |阅读模式 来自 湖北武汉
本人在用无网格方法做粘弹性分析，由于利用弹性-粘弹性相应原理，在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)*[GG1+2*GG2 ,GG2-GG1,0;GG2-GG1,GG1+2*GG2,0;0,0,(2*GG1+GG2)/2]),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;
[x,conn1,conn2,numnod] = 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;
[xc,conn,numcell,numq] = mesh2(Lb,D,ndivlq,ndivwq);
% SET UP GAUSS POINTS, WEIGHTS, AND JACOBIAN FOR EACH CELL
quado = 4;
[gauss] = pgauss(quado);
numq2 = numcell*quado^2;
gs = zeros(4,numq2);
[gs] = 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);
   [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
   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);
   [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);
   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);
   [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 =[k GG;GG' zeros(lthu*2)];
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);
   [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
%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);
   [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);
   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)
回复

使用道具 举报

qibbxxt
发表于 2012-6-11 20:11:48 | 显示全部楼层 来自 河北廊坊
Simdroid开发平台
换个64位系统试一试
回复 赞 不支持

使用道具 举报

ccz19860422
 楼主| 发表于 2012-6-21 10:21:14 | 显示全部楼层 来自 湖北武汉
qibbxxt 发表于 2012-6-11 20:11
换个64位系统试一试

就是在64位的服务器上计算的，还是不行。
回复 赞 不支持

使用道具 举报

qibbxxt
发表于 2012-6-25 14:56:31 | 显示全部楼层 来自 河北廊坊
ccz19860422 发表于 2012-6-21 10:21
就是在64位的服务器上计算的，还是不行。

正常来说，应该是i的程序的问题
回复 赞 不支持

使用道具 举报

ccz19860422
 楼主| 发表于 2012-6-28 09:09:14 | 显示全部楼层 来自 湖北武汉
qibbxxt 发表于 2012-6-25 14:56
正常来说，应该是i的程序的问题

能不能请高手指点一二啊:P
回复 赞 不支持

使用道具 举报

pasuka
发表于 2012-6-28 21:03:42 | 显示全部楼层 来自 上海杨浦区
本帖最后由 pasuka 于 2012-6-28 21:04 编辑
ccz19860422 发表于 2012-6-28 09:09
能不能请高手指点一二啊

符号运算的表达式推导出来后，简化后对照着再编程用matlab计算，不要直接vpa或者subs
如果这样还不行的话，根据符号推导的表达式用C、C++、Fortran编程计算
如果符号运算就报内存不足的话，换mathematics推导公式，其余同上

点评

rocwoods
可以用matlabFunction这个函数将符号表达式转化成数值函数句柄来计算，符号计算用来推导公式，计算交给数值。  发表于 2012-6-29 09:29
回复 赞 不支持

使用道具 举报

ccz19860422
 楼主| 发表于 2012-7-2 08:52:56 | 显示全部楼层 来自 湖北武汉
十分感谢！
回复 赞 不支持

使用道具 举报

返回列表 发新帖
高级模式
B Color Image Link Quote Code Smilies
您需要登录后才可以回帖 登录 | 注册

本版积分规则

Archiver|小黑屋|联系我们|仿真互动网 ( 京ICP备15048925号-7 )

GMT+8, 2024-7-3 14:12 , Processed in 0.038048 second(s), 15 queries , Gzip On, MemCache On.

Powered by Discuz! X3.5 Licensed

© 2001-2024 Discuz! Team.

快速回复 返回顶部 返回列表