24小时热门版块排行榜

返回列表

明明那天

新虫 (初入文坛)

应助: 0 (幼儿园)
金币: 237.4
散金: 1
红花: 2
帖子: 28
在线: 15.1小时
虫号: 1366161
注册: 2011-08-11
专业: 机械动力学

[求助] matlab 数值积分速度上不去

我是新手，在用matlab计算数值积分的程序，可能是存在奇点问题，运行速度非常慢。不知有没有高手可以帮忙指导一下。谢谢了！

clear all
clc
syms x y u1 u2 k1 k2 r1 r2 ss1 ss2 a1 c1 s1 cc1 c_c1
A11=1.58e8; A66=1.05e7;A22=7.99e7;

A12=5.625e5;
D1=A66*A22/(A12^2-A11*A22);
D2=A66*A11/(A12^2-A11*A22);
D3=-2*A66*A12/(A12^2-A11*A22)+1;
u=[1.875;4.694;7.855;10.996;14.173;17.279];%x
q=[4.730;7.853;10.996;14.137;17.279;20.378];%y
FUN5=0;

X=[];Y=[];int1=zeros(36,36);int2=zeros(36,36);int3=zeros(36,36);row=0;
for mode_num2=1:6
   u_2=u(mode_num2);k_2=q(mode_num2);
r_2=(cos(u_2)-cosh(u_2))/(sinh(u_2)-sin(u_2));
   s_2=(sinh(k_2)-sin(k_2))/(cos(k_2)+cosh(k_2));
X=[X,cosh(u_2*x)+cos(u_2*x)+r_2*(sin(u_2*x)+sinh(u_2*x))];
Y=[Y,cosh(k_2*y)-cos(k_2*y)+s_2*(sin(k_2*y)-sinh(k_2*y))];
end
for mode_num=1:6
   W5=0;F5=0;L=0;row=row+1;
   u_1=u(mode_num);k_1=q(mode_num);
   r_1=(cos(u_1)-cosh(u_1))/(sinh(u_1)-sin(u_1));
   s_1=(sinh(k_1)-sin(k_1))/(cos(k_1)+cosh(k_1));
  X5=[sinh(u1*x)^2*u1^2,-2*sinh(u1*x)*u1^2*sin(u1*x),2*sinh(u1*x)*u1^2*r1*cos(u1*x),2*sinh(u1*x)*u1^2*r1*cosh(u1*x),sin(u1*x)^2*u1^2,-2*sin(u1*x)*u1^2*r1*cos(u1*x),-2*sin(u1*x)*u1^2*r1*cosh(u1*x),r1^2*cos(u1*x)^2*u1^2,2*r1^2*cos(u1*x)*u1^2*cosh(u1*x),r1^2*cosh(u1*x)^2*u1^2];
  Y5=[sinh(k1*y)^2*k1^2,2*sinh(k1*y)*k1^2*sin(k1*y),2*sinh(k1*y)*k1^2*ss1*cos(k1*y),-2*sinh(k1*y)*k1^2*ss1*cosh(k1*y),sin(k1*y)^2*k1^2,2*sin(k1*y)*k1^2*ss1*cos(k1*y),-2*sin(k1*y)*k1^2*ss1*cosh(k1*y),ss1^2*cos(k1*y)^2*k1^2,-2*ss1^2*cos(k1*y)*k1^2*cosh(k1*y),ss1^2*cosh(k1*y)^2*k1^2];

  for i=1:10
for  j=1:10
   W5=W5-u1^2*k1^2*X5(i)*Y5(j);
end
  end

  syms c_c1 c_c101 c_c201 c_c301 c_c401 c_c501
c_c=[c_c1, c_c101, c_c201, c_c301, c_c401, c_c501];
if mode_num==1
   str5=c_c1;
      for i=2:100
syms(['c_c',num2str(i)])
str5=[str5,['c_c',num2str(i)]];
      end

else
   i_1=2+(mode_num-1)*100;i_2=100+(mode_num-1)*100;
   str5=c_c(mode_num);
for  i=i_1:i_2
syms(['c_c',num2str(i)])
str5=[str5,['c_c',num2str(i)]];
end
end
%str5
x_5=[cosh(u1*x),cos(u1*x),sin(u1*x),sinh(u1*x)];
  y_5=[cosh(k1*y),cos(k1*y),sin(k1*y),sinh(k1*y)];
  xm1=[cosh(u_1*x),cos(u_1*x),sin(u_1*x),sinh(u_1*x)];
  ym1=[cosh(k_1*y),cos(k_1*y),sin(k_1*y),sinh(k_1*y)];
  intmode1=zeros(100,36);intmode2=zeros(100,36);
  intmode3=zeros(100,36);
  for i=1:4
   for j=i:4
      for r=1:4
            for k=r:4
      L=L+1;
      F5=F5+str5(L)*x_5(i)*x_5(j)*y_5(r)*y_5(k);
   for mu=1:6
      for s=1:6
      m=s+(mu-1)*6;

Ax1=diff(X(mu),x,2);
Ay1=Y(mu);
Fx1=xm1(i)*xm1(j);
Fy1=diff(ym1(r)*ym1(k),y,2);
Sx1=vectorize(Fx1*X(s)*Ax1);
Sy1=vectorize(Fy1*Y(s)*Ay1);
intx1=quadl(Sx1,0,1) ;
inty1=quadl(Sy1,0,1) ;
intmode1(L,m)=intx1*inty1;

Ax2=diff(X(mu),x);
Ay2=diff(Y(mu),y);
Fx2=diff(xm1(i)*xm1(j),x);
Fy2=diff(ym1(r)*ym1(k),y);
Sx2=vectorize(Fx2*X(s)*Ax2);
Sy2=vectorize(Fy2*Y(s)*Ay2);
intx2=quadl(Sx2,0,1) ;
inty2=quadl(Sy2,0,1) ;
intmode2(L,m)=intx2*inty2;

Ax3=X(mu);
Ay3=diff(Y(mu),y,2);
Fx3=diff(xm1(i)*xm1(j),x,2);
Fy3=ym1(r)*ym1(k);
Sx3=vectorize(Fx3*X(s)*Ax3);
Sy3=vectorize(Fy3*Y(s)*Ay3);
intx3=quadl(Sx3,0,1) ;
inty3=quadl(Sy3,0,1) ;
intmode3(L,m)=intx3*inty3;
      end
   end
            end
      end
   end
  end

回复此楼