Need to remake this optimize function to be efficient takes way to long currently

%% other given parameters
clear; clc;
E = 70*10^9
G = 30*10^9
L = 1
data = 1
%%Simple Script 
[h b] = optimize(data,E,G,L) %% optimized hieght and base to get highest Pcr



%% Problem function 

function [h b] = optimize(data,E,G,L) 
    t = .002 %% constant 
    A1 = 6*10^-4 %% area has to be smaller than this 
    i = 1 %counter 

    %%for loop that find all pvals for combination of h and bs 
    if data  == 1
    for h = .001:.001:.5
        for b = .001:.001:.5
            [A,Ixx,Iyy,Xs,Ys,J,C_w,I0] = Ibeam(t,h,b) 
            if A < A1
                [Pcritical_matrix] = Pcrit(E,G,L,J,Iyy,Ixx,A,I0,C_w,Xs,Ys);
                pval = min(Pcritical_matrix) %%finds the min of this matrix but want to find the max possible of this min
                m(:,i) = [pval h b]
                i = i+1
            end 
        end 
    end 
    end 

   [M I] =  max(m(1,:))
   h = m(2,I)
   b = m(3,I)
end 

%% other functions 

function [A,Ixx,Iyy,Xs,Ys,J,C_w,I0] = Ibeam(t,h,b)
Cy = h/2;
Cx = b/2;
Ixx = 2*((b*t^3/12)+(t*b*((h/2)+(t/2))^2))+(t*((h))^3/12)
Iyy = 2*(t*b^3/12)+((h*t^3/12))
I0 = Ixx + Iyy;
J = (2*b*t^3+h*t^3)/3
C_w = t*h^2*b^3/24
A = 2*(b*t)+((h)*t)
Xs = 0
Ys = 0
end 
        function [Pcritical_matrix] = Pcrit(E,G,L,J,Iyy,Ixx,A,I0,C_w,Xs,Ys)
i = 1
Pcritical_matrix = zeros(1,3);
Py = pi^2*E*Iyy/(L^2)
Px = pi^2*E*Ixx/(L^2)
P_alpha = (A /I0)*(G*J+(E*C_w*pi^2/L^2))


if Xs == 0 && Ys == 0 
    Pcr1 = Px;
    Pcr2 = Py;
    Pcr3 = P_alpha;
    Pcritical_matrix = [Pcr1 Pcr2 Pcr3];
else 

syms pCr;
m = [0 (pCr-Px) (pCr*.076);(pCr-Py) 0 0; 0 (.076*pCr) (.0089*(pCr-P_alpha))];
pCrit = solve(det(m));
Pcritical_matrix = vpa(pCrit);

 
end 
end 

Answer 1

[Torsten]

You can easily solve when det(m) = 0 in advance and insert the formula in your code. This will enhance speed enormously:

syms pCr Px Py c1 c2 P_alpha
m = [0 (pCr-Px) (pCr*c1);(pCr-Py) 0 0; 0 (c1*pCr) (c2*(pCr-P_alpha))];
pCrit = solve(det(m))

pCrit = 

with c1 = 0.076 and c2 = 0.0089.

Answer 2

[the cyclist]

The only thing that is taking a long time is that your code is displaying almost every output to the command window.

Put semicolons at the ends of the lines, and that output will be suppressed. Your code ran in less than one second for me.

Answer 3

[Walter Roberson]

Q = @(v) sym(v);

E = Q(70)*Q(10)^9;
G = Q(30)*Q(10)^9;
L = Q(1);

t = Q(2)/1000; %% constant 
A1 = Q(6)/10^4; %% area has to be smaller than this 

syms h b positive

Cy = h/2;
Cx = b/2;
Ixx = 2*((b*t^3/12)+(t*b*((h/2)+(t/2))^2))+(t*((h))^3/12);
Iyy = 2*(t*b^3/12)+((h*t^3/12));
I0 = Ixx + Iyy;
J = (2*b*t^3+h*t^3)/3;
C_w = t*h^2*b^3/24;
A = 2*(b*t)+((h)*t);
Xs = 0;
Ys = 0;

Py = Q(pi)^2*E*Iyy/(L^2);
Px = Q(pi)^2*E*Ixx/(L^2);
P_alpha = (A /I0)*(G*J+(E*C_w*Q(pi)^2/L^2))

P_alpha = 

syms pCr
m = [0 (pCr-Px) (pCr*Q(76)/1000);(pCr-Py) 0 0; 0 (Q(76)/1000*pCr) (Q(89)/10000*(pCr-P_alpha))];
pCrit = solve(det(m));

pCrit

pCrit = 

simplify(expand(pCrit(2) - pCrit(3)), 'steps', 20)

ans = 

Each of the terms for pCrit(2) - pCrit(3) is positive for positive real b, h. Therefore pCrit(2) > pCrit(3) under those conditions so min() between those two is pCrit(3) . We have already seen that the polynomial term is much larger than the other values; if that continues to hold then it follows that pCrit(3) is the minimum of the three terms. The task would then be to maximize pCrit(3) over those ranges.

However... this analysis really only follows if (Xs == 0 & Ys == 0) is false because you use a different calculation in that particular case. But you do set Xs = 0 and Ys = 0, so really you should do a separate analysis for the alternate flow, the Px Py Palpha flow.