Is there any way to speed this up?
1 view (last 30 days)
Show older comments
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)
0 Comments
Answers (2)
Arthur Vidal
on 11 Mar 2018
Work with syms is usually slow. Try to do not use it and make the calculations directly with numbers instead.
0 Comments
Walter Roberson
on 11 Mar 2018
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;
See Also
Categories
Find more on Optimization Toolbox in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!