Is there any way to speed this up?

1 次查看（过去 30 天）
Timothy Cerrie
Timothy Cerrie 2018-3-10
评论： Timothy Cerrie 2018-3-12
Is there any why this code can be optimized currently takes about 45 min to run through everything.
clc, clear
syms x y
a=0.4826; %length in x
b=0.19; %length in y
v=0.33; %poisson's ratio
h=0.00127; %thickness
Xl=a/2;%location of x for analysis
Yl=b/2;%location of y analysis
Zl=h/2;%location of z from nuetral plane for analysis
rho = 2700.*h; %volumetric density
E= 68900000000; %Youngs modulas
FT= 1398294.769; %force per unit in z-direction
D = (E.*h.^3)./(12.*(1-v.^2));
gamma = [0 4.73 7.85 9.42];
X=sym(zeros(1,4));
Y=sym(zeros(1,4));
for n = 1:length(gamma)
if n==2,4;
X(n) = cos(gamma(n)*(x/a-1/2))+(sin(gamma(n)/2)/sinh(gamma(n)/2))*cosh(gamma(n)*(x/a-1/2));% even
Y(n) = cos(gamma(n)*(y/b-1/2))+(sin(gamma(n)/2)/sinh(gamma(n)/2))*cosh(gamma(n)*(y/b-1/2));% even
else
X(n) = sin(gamma(n)*(x/a-1/2))-(sin(gamma(n)/2)/sinh(gamma(n)/2))*sinh(gamma(n)*(x/a-1/2));% odd
Y(n) = sin(gamma(n)*(y/b-1/2))-(sin(gamma(n)/2)/sinh(gamma(n)/2))*sinh(gamma(n)*(y/b-1/2));% odd
end
end
count=1;
mode_shapes=sym(zeros(9,1));
freq=zeros(9,1);
for i=2:4
if i==2
G_x = 1.506;
H_x = 1.248;
J_x = 1.248;
else
G_x = i-0.5;
H_x = (i-0.5).^2.*(1 - (2./(i-0.5).*pi));
J_x = (i-0.5).^2.*(1 - (2./(i-0.5).*pi));
end
for j=2:4
mode_shapes(count,1) = X(i)*Y(j);
if j==2
G_y = 1.506;
H_y = 1.248;
J_y = 1.248;
freq(count,1) = sqrt(((pi.^4.*D)./(a.^4.*rho)).*((G_x.^4+(G_y.^4.*(a./b).^4)) +((2.*(a./b).^2).*((v.*H_x.*H_y)+((1-v).*J_x.*J_y))))); % frequency9
else
G_y = j-0.5;
H_y = (j-0.5).^2.*(1 - (2./(j-0.5).*pi));
J_y = (j-0.5).^2.*(1 - (2./(j-0.5).*pi));
freq(count,1) = sqrt(((pi.^4.*D)./(a.^4.*rho)).*((G_x.^4+(G_y.^4.*(a./b).^4)) +((2.*(a./b).^2).*((v.*H_x.*H_y)+((1-v).*J_x.*J_y))))); % frequency
end
count=count+1;
end
end
[~, C]=sort(freq); %sorting frequencies
new_mode_shapes=mode_shapes(C, :); % sorting mode shapes in order to frequencies
MO=new_mode_shapes;
syms z w(x,y)
G=E/(2*(1+v));
U=-z.*diff(w,x);
V=-z.*diff(w,y);
ex=diff(U,x);
ey=diff(V,y);
gxy=diff(V,x)+diff(U,y);
Sx=(1./(1-v.^2)).*(E.*ex+v.*E.*ey);
Sy=(1./(1-v.^2)).*(E.*ey+v.*E.*ex);
Txy=G.*gxy;
Nx=int(Sx,z,-h/2,h/2);
Ny=int(Sy,z,-h/2,h/2);
Nxy=int(Txy,z,-h/2,h/2);
Mx=int(Sx*z,z,-h/2,h/2);
My=int(Sy*z,z,-h/2,h/2);
Mxy=int(Txy*z,z,-h/2,h/2);
Qx=diff(Mx,x)+diff(Mxy,y);
Qy=diff(My,y)+diff(Mxy,x);
qz=diff(Qx,x)+diff(Qy,y)+diff(Nx*diff(w,x),x)+diff(Ny*diff(w,y),y)+diff(Nxy*diff(w,y),x)+diff(Nxy*diff(w,x),y);
qx=diff(Nx,x)+diff(Nxy,y)-diff(Qx*diff(w,x),x)-diff(Qy*diff(w,x),y);
qy=diff(Ny,y)+diff(Nxy,x)-diff(Qx*diff(w,y),x)-diff(Qy*diff(w,y),y);
P=zeros(9,9);
qrz=zeros(9,1);
%%Glerkin method
for ii=1:9
Wo=MO(ii);
for jj=1:9
Wi=MO(jj);
qZ=subs(qz,w,Wi);
P(ii,jj)=int(int(Wo.*(qZ),x,0,b),y,0,a);
end
qrz(ii)=int(int(Wo.*FT,x,0,b),y,0,a);
end
CC=inv(P)*qrz;
WXY=sym(zeros(9,1));
for ij=1:9
WXY(ij)=CC(ij).*MO(ij);
end
W=cumsum(WXY);
Deflection=double(subs(W(3),[x y],[Xl,Yl]));
sx=subs(Sx,[z,w],[Zl,W(3)]);
sy=subs(Sy,[z,w],[Zl,W(3)]);
txy=subs(Txy,[z,w],[Zl,W(3)]);
sxx=subs(sx,[x,y],[Xl,Yl]);
syy=subs(sy,[x,y],[Xl,Yl]);
txxyy=subs(txy,[x,y],[Xl,Yl]);
%%von misses
Vm=double(sqrt(0.5.*((sxx-syy).^2+syy.^2+sxx.^2+(6.*txxyy.^2))));%our stess
Yeild=3.1e+08;%actual yeild stress
K=double(Yeild/Vm);%conentration factor
fprintf( 'Total deflection = %f m \n',Deflection)
fprintf('Total stress = %f Pa \n',Vm)
fprintf('Stress concentration from dent = %f \n',K)

请先登录，再回答此问题。

回答（2 个）

Arthur Vidal
Arthur Vidal 2018-3-11
Work with syms is usually slow. Try to do not use it and make the calculations directly with numbers instead.
Walter Roberson
Walter Roberson 2018-3-11
You have
if n==2,4;
The meaning of that is the same as
if n==2
4;
which executes 4 and throws the result away because of the semi-colon.
The closest MATLAB syntaxes to what you appear to be trying to do are:
if any(n==[2,4])
or
if n == 2 || n == 4
or
if ismember(n, [2, 4])
In your section
if j==2
G_y = 1.506;
H_y = 1.248;
J_y = 1.248;
freq(count,1) = sqrt(((pi.^4.*D)./(a.^4.*rho)).*((G_x.^4+(G_y.^4.*(a./b).^4)) +((2.*(a./b).^2).*((v.*H_x.*H_y)+((1-v).*J_x.*J_y))))); % frequency9
else
G_y = j-0.5;
H_y = (j-0.5).^2.*(1 - (2./(j-0.5).*pi));
J_y = (j-0.5).^2.*(1 - (2./(j-0.5).*pi));
freq(count,1) = sqrt(((pi.^4.*D)./(a.^4.*rho)).*((G_x.^4+(G_y.^4.*(a./b).^4)) +((2.*(a./b).^2).*((v.*H_x.*H_y)+((1-v).*J_x.*J_y))))); % frequency
end
your J_y appears to be identical to your H_y in both cases. Why not just calculate H_y and then assign J_y = H_y ?
The purpose of using int() is to find closed form solutions to the integral if a closed form is available. When I test your code, it appears that closed form solutions are available, but for example the exact solution to P(9,9) is
(25303565968580561233*cos(12771/25400)*sin(12771/25400)*((19*cos(473/200)*sin(473/200))/946 + (19*cos(24123/2500)*sin(24123/2500))/946 + 2413/10000))/73344658800000 - (25303565968580561233*cos(24123/2500)*sin(24123/2500))/1961004561600000 - (64150623052873913699*cos(473/200)*sin(473/200)*((19*cos(473/200)*sin(473/200))/946 + (19*cos(24123/2500)*sin(24123/2500))/946 + 2413/10000))/73344658800000 - (25303565968580561233*cos(473/200)*sin(473/200))/1961004561600000 - (305051014177721440629*sin(473/200)^2*((19*sinh(473/100))/1892 + (19*sinh(24123/1250))/1892 + 2413/10000))/(887121111200000*sinh(473/200)^2) - (11968586703138605463209*sin(473/200)^4*((19*sinh(473/100))/1892 + (19*sinh(24123/1250))/1892 + 2413/10000))/(18629543335200000*sinh(473/200)^4) - (9187329000435377858209*sin(473/200)^3*((19*cos(473/200)*sinh(473/200))/946 + (19*cosh(473/200)*sin(473/200))/946 + (19*cos(24123/2500)*sinh(24123/2500))/946 + (19*cosh(24123/2500)*sin(24123/2500))/946))/(9314771667600000*sinh(473/200)^3) - (305051014177721440629*sin(473/200)^2*((19*cos(473/200)*sin(473/200))/946 + (19*cos(24123/2500)*sin(24123/2500))/946 + 2413/10000))/(887121111200000*sinh(473/200)^2) - (9187329000435377858209*sin(473/200)*((19*cos(473/200)*sinh(473/200))/946 + (19*cosh(473/200)*sin(473/200))/946 + (19*cos(24123/2500)*sinh(24123/2500))/946 + (19*cosh(24123/2500)*sin(24123/2500))/946))/(9314771667600000*sinh(473/200)) - (19423528542146676233*cosh(473/200)*sin(473/200)^2*((19*cos(473/200)*sin(473/200))/946 + (19*cos(24123/2500)*sin(24123/2500))/946 + 2413/10000))/(36672329400000*sinh(473/200)) - (644928148367275773*sin(473/200)^2*sinh(473/100)*((19*cos(473/200)*sin(473/200))/946 + (19*cos(24123/2500)*sin(24123/2500))/946 + 2413/10000))/(6985205600000*sinh(473/200)^2) + (644928148367275773*sin(473/200)^2*sinh(12771/12700)*((19*cos(473/200)*sin(473/200))/946 + (19*cos(24123/2500)*sin(24123/2500))/946 + 2413/10000))/(6985205600000*sinh(473/200)^2) - (52390548200006143699*cos(473/200)*sin(473/200)^3*((19*sinh(473/100))/1892 + (19*sinh(24123/1250))/1892 + 2413/10000))/(73344658800000*sinh(473/200)^2) - (19423528542146676233*cosh(473/200)*sin(473/200)^4*((19*sinh(473/100))/1892 + (19*sinh(24123/1250))/1892 + 2413/10000))/(36672329400000*sinh(473/200)^3) - (25303565968580561233*sin(473/200)^4*sinh(473/100)*((19*sinh(473/100))/1892 + (19*sinh(24123/1250))/1892 + 2413/10000))/(146689317600000*sinh(473/200)^4) + (25303565968580561233*sin(473/200)^4*sinh(12771/12700)*((19*sinh(473/100))/1892 + (19*sinh(24123/1250))/1892 + 2413/10000))/(146689317600000*sinh(473/200)^4) - (19423528542146676233*cos(473/200)*sin(473/200)^2*((19*cos(473/200)*sinh(473/200))/946 + (19*cosh(473/200)*sin(473/200))/946 + (19*cos(24123/2500)*sinh(24123/2500))/946 + (19*cosh(24123/2500)*sin(24123/2500))/946))/(12224109800000*sinh(473/200)) - (19423528542146676233*cosh(473/200)*sin(473/200)^3*((19*cos(473/200)*sinh(473/200))/946 + (19*cosh(473/200)*sin(473/200))/946 + (19*cos(24123/2500)*sinh(24123/2500))/946 + (19*cosh(24123/2500)*sin(24123/2500))/946))/(18336164700000*sinh(473/200)^2) - (19423528542146676233*sin(473/200)^3*sinh(473/100)*((19*cos(473/200)*sinh(473/200))/946 + (19*cosh(473/200)*sin(473/200))/946 + (19*cos(24123/2500)*sinh(24123/2500))/946 + (19*cosh(24123/2500)*sin(24123/2500))/946))/(73344658800000*sinh(473/200)^3) + (19423528542146676233*sin(473/200)^3*sinh(12771/12700)*((19*cos(473/200)*sinh(473/200))/946 + (19*cosh(473/200)*sin(473/200))/946 + (19*cos(24123/2500)*sinh(24123/2500))/946 + (19*cosh(24123/2500)*sin(24123/2500))/946))/(73344658800000*sinh(473/200)^3) + (19423528542146676233*cos(12771/25400)*sin(473/200)^2*sinh(12771/25400)*((19*cos(473/200)*sinh(473/200))/946 + (19*cosh(473/200)*sin(473/200))/946 + (19*cos(24123/2500)*sinh(24123/2500))/946 + (19*cosh(24123/2500)*sin(24123/2500))/946))/(18336164700000*sinh(473/200)^2) + (19423528542146676233*cosh(12771/25400)*sin(473/200)^2*sin(12771/25400)*((19*cos(473/200)*sinh(473/200))/946 + (19*cosh(473/200)*sin(473/200))/946 + (19*cos(24123/2500)*sinh(24123/2500))/946 + (19*cosh(24123/2500)*sin(24123/2500))/946))/(18336164700000*sinh(473/200)^2) + (19423528542146676233*cos(12771/25400)*sin(473/200)*sinh(12771/25400)*((19*cos(473/200)*sin(473/200))/946 + (19*cos(24123/2500)*sin(24123/2500))/946 + 2413/10000))/(36672329400000*sinh(473/200)) + (19423528542146676233*cosh(12771/25400)*sin(473/200)*sin(12771/25400)*((19*cos(473/200)*sin(473/200))/946 + (19*cos(24123/2500)*sin(24123/2500))/946 + 2413/10000))/(36672329400000*sinh(473/200)) + (19423528542146676233*cos(12771/25400)*sin(473/200)*sin(12771/25400)*((19*cos(473/200)*sinh(473/200))/946 + (19*cosh(473/200)*sin(473/200))/946 + (19*cos(24123/2500)*sinh(24123/2500))/946 + (19*cosh(24123/2500)*sin(24123/2500))/946))/(36672329400000*sinh(473/200)) + (644928148367275773*cos(12771/25400)*sin(473/200)^2*sin(12771/25400)*((19*sinh(473/100))/1892 + (19*sinh(24123/1250))/1892 + 2413/10000))/(3492602800000*sinh(473/200)^2) + (19423528542146676233*cos(12771/25400)*sin(473/200)^3*sinh(12771/25400)*((19*sinh(473/100))/1892 + (19*sinh(24123/1250))/1892 + 2413/10000))/(36672329400000*sinh(473/200)^3) + (19423528542146676233*cosh(12771/25400)*sin(473/200)^3*sin(12771/25400)*((19*sinh(473/100))/1892 + (19*sinh(24123/1250))/1892 + 2413/10000))/(36672329400000*sinh(473/200)^3) - 11968586703138605463209/77204904000000000
I have to ask: do you really need that kind of accuracy? Because you assign to P which you initialized as double precision, so there is no point in holding that much accuracy. You might as well do a much faster numeric integration.
And then you do
CC=inv(P)*qrz;
which risks numeric problems. Avoid using inv(). Your expression would be better written
CC = P\qrz;

类别

在 Help Center 和 File Exchange 中查找有关 Formula Manipulation and Simplification 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by