Is there any way to speed this up?

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)

Answer 1

[Arthur Vidal]

Work with syms is usually slow. Try to do not use it and make the calculations directly with numbers instead.

Answer 2

[Walter Roberson]

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;

[Timothy Cerrie]

Thank you Walter this answer really helped.