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;
  댓글 수: 1
Timothy Cerrie
Timothy Cerrie 2018년 3월 12일
Thank you Walter this answer really helped.

댓글을 달려면 로그인하십시오.

카테고리

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