How to speed up the iteration of the for loop and make some further improvement on the below code?

Question

Said Bolluk 2020-9-10

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/591814-how-to-speed-up-the-iteration-of-the-for-loop-and-make-some-further-improvement-on-the-below-code

评论： Said Bolluk 2020-9-13

Hello,

The moment_curvature function is to analyze the ductile behavior of a reinforced concrete member. There are some significant strain values, given in ecu_list, affecting the performance of the section. To calculate the curvature and moment values specific to these strains, I need to calculate the strains and stresses by finding the neutral axis depth, c. After that, I can calculate the forces. The summation of the internal reaction forces, which are provided by the steel and concrete layers, must be equal to the total sum of the external forces -zero (0) in this case.

So, my approach is to iterate over the c value and find the exact depth. As soon as the F_total closes to zero, I accepted that value of c and performed accordingly. However, these operations take too long and are a bit away from reality. I need the higher precision, but more importantly, I should speed up the running time of this code. Thus, I am open to any suggestions, and I will appreciate if I can get some help.

Thank you.

function moment_curvature
    %%% section property
    h = 600; %height
    b = 300; %width
    %%% concrete properties:
    fcu = 30; %MPa
    %%% reinforcing steel properties:
    fy = 420; %MPa
    Es = 200e3; %MPa
    %%% layers input = [area, effective depth, moment arm] all in mm:
    steel1 = [2000, 560, h/2 - 570];
    steel2 = [1000, 40, h/2 - 30];
    %%% strain values
    ecu_list = [2.5e-4, 5e-4, 1e-3, 1.5e-3, 2e-3, 2.5e-3, 3e-3, 3.5e-3, 4e-3];
    %%% alpha&beta are the cooefficients specific to the given strain values
    alpha = [0.178, 0.336, 0.595, 0.779, 0.889, 0.929, 0.934, 0.925, 0.910];
    beta = [0.674, 0.682, 0.700, 0.722, 0.750, 0.785, 0.818, 0.849, 0.874];
    N = length(ecu_list);
    curvature_list = zeros(1, N);
    moment_list = zeros(1, N);
    for ecu_idx = 1 : N
        ec_top = ecu_list(ecu_idx);
        for c = 0 : 1e-7 : 2*h %the neutral axis depth where there is no stress
            %%% steel layer 1
            es1 = (c - steel1(2)) * ec_top / c; %mm/mm
            fs1 = Es * es1; %MPa
            if abs(fs1) > fy
                if fs1 < 0
                    fs1 = -fy;
                else
                    fs1 = fy;
                end
            end
            Fs1 = fs1 * steel1(1) * 1e-3; %kN
            m1 = Fs1 * steel1(3) * 1e-3; %kN.m
            %%% steel layer 2
            es2 = (c - steel2(2)) * ec_top / c; %mm/mm
            fs2 = Es * es2; %MPa
            if abs(fs2) > fy
                if fs2 < 0
                    fs2 = -fy;
                else
                    fs2 = fy;
                end
            end
            Fs2 = fs2 * steel2(1) * 1e-3; %kN
            m2 = Fs2 * steel2(3) * 1e-3; %kN.m
            %%% concrete layer
            Fc = alpha(ecu_idx) * fcu * beta(ecu_idx) * c * b * 1e-3;
            mc = Fc * (h/2 - beta(ecu_idx)*c/2) * 1e-3; %kN.m
            %%% F_total must be equal to the sum of external forces
            F_total = Fs1 + Fs2 + Fc;
            %%% setting F_total to zero gives the neutral axis depth, c
            if -0.0001 < F_total && 0.0001 > F_total
                break;
            end
            moment = m1 + m2 + mc; %kN.m
            curvature = ec_top / c; %1/mm
            moment_list(ecu_idx) = moment;
            curvature_list(ecu_idx) = curvature;
        end
    end
    plot(curvature_list, moment_list);
    drawnow;
end

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Answer 1

Alan Stevens 2020-9-10

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/591814-how-to-speed-up-the-iteration-of-the-for-loop-and-make-some-further-improvement-on-the-below-code#answer_492763

在 MATLAB Online 中打开

The following, which makes use of fzero to find the value of c that makes F_total equal to zero, runs a lot faster. I've no idea if the results are sensible!

    %%% section property
    h = 600; %height
    b = 300; %width
    %%% concrete properties:
    fcu = 30; %MPa
    %%% reinforcing steel properties:
    fy = 420; %MPa
    Es = 200e3; %MPa
    materialProps = [h,b,fcu,fy,Es];
    %%% layers input = [area, effective depth, moment arm] all in mm:
    steel1 = [2000, 560, h/2 - 570];
    steel2 = [1000, 40, h/2 - 30];
    Layers = [steel1; steel2];
    %%% strain values
    ecu_list = [2.5e-4, 5e-4, 1e-3, 1.5e-3, 2e-3, 2.5e-3, 3e-3, 3.5e-3, 4e-3]; 
    %%% alpha&beta are the cooefficients specific to the given strain values
    alpha = [0.178, 0.336, 0.595, 0.779, 0.889, 0.929, 0.934, 0.925, 0.910];
    beta = [0.674, 0.682, 0.700, 0.722, 0.750, 0.785, 0.818, 0.849, 0.874];
    alphabeta = [alpha; beta];
    
    N = length(ecu_list);
    curvature_list = zeros(1, N);
    moment_list = zeros(1, N);
    c = h; % Initial guess
for ecu_idx = 1 : N
    ec_top = ecu_list(ecu_idx);
    
    c0 = c; %  Use previous converged value as guess for next ecu_idx
    c = fzero(@Ffn, c0,[],ec_top, ecu_idx,materialProps, Layers, alphabeta);
   
       %%% steel layer 1
    es1 = (c - steel1(2)) * ec_top / c; %mm/mm
    fs1 = Es * es1; %MPa
    if abs(fs1) > fy
        if fs1 < 0
            fs1 = -fy;
        else
            fs1 = fy;
        end
    end
    Fs1 = fs1 * steel1(1) * 1e-3; %kN
    m1 = Fs1 * steel1(3) * 1e-3; %kN.m
    %%% steel layer 2
    es2 = (c - steel2(2)) * ec_top / c; %mm/mm
    fs2 = Es * es2; %MPa
    if abs(fs2) > fy
        if fs2 < 0
            fs2 = -fy;
        else
            fs2 = fy;
        end
    end
    Fs2 = fs2 * steel2(1) * 1e-3; %kN
    m2 = Fs2 * steel2(3) * 1e-3; %kN.m
    %%% concrete layer
    Fc = alpha(ecu_idx) * fcu * beta(ecu_idx) * c * b * 1e-3;
    mc = Fc * (h/2 - beta(ecu_idx)*c/2) * 1e-3; %kN.m
    moment = m1 + m2 + mc; %kN.m
    curvature = ec_top / c; %1/mm
    moment_list(ecu_idx) = moment;
    curvature_list(ecu_idx) = curvature;
end
    
    plot(curvature_list, moment_list,'o');
    drawnow;
    function F_total = Ffn(c, ec_top, ecu_idx,materialProps, Layers,alphabeta)   
    
            h = materialProps(1);
            b = materialProps(2);
            fcu = materialProps(3);
            fy = materialProps(4);
            Es = materialProps(5);
            steel1 = Layers(1,:);
            steel2 = Layers(2,:);
            alpha = alphabeta(1,:);
            beta = alphabeta(2,:);
    
            %%% steel layer 1
            es1 = (c - steel1(2)) * ec_top / c; %mm/mm
            fs1 = Es * es1; %MPa
            if abs(fs1) > fy
                if fs1 < 0
                    fs1 = -fy;
                else
                    fs1 = fy;
                end
            end
            Fs1 = fs1 * steel1(1) * 1e-3; %kN
            %%% steel layer 2
            es2 = (c - steel2(2)) * ec_top / c; %mm/mm
            fs2 = Es * es2; %MPa
            if abs(fs2) > fy
                if fs2 < 0
                    fs2 = -fy;
                else
                    fs2 = fy;
                end
            end
            Fs2 = fs2 * steel2(1) * 1e-3; %kN
            %%% concrete layer
            Fc = alpha(ecu_idx) * fcu * beta(ecu_idx) * c * b * 1e-3;
            %%% F_total must be equal to the sum of external forces
            F_total = Fs1 + Fs2 + Fc;
            %%% setting F_total to zero gives the neutral axis depth, c
            
    end

10 个评论
显示 8更早的评论隐藏 8更早的评论

Alan Stevens 2020-9-10

在 MATLAB Online 中打开

% The first part from here
    h = 600; %height
% to here
    moment_list = zeros(1, N);
% is basically just your original code,except that I've added 
    materialProps = [h,b,fcu,fy,Es];
    Layers = [steel1; steel2];
    alphabeta = [alpha; beta];
% to collect together data that will later be passed to a function.
% Then, within your loop
    for ecu_idx = 1 : N
        ...
    end
% I've inserted the fzero function.  This finds the value of a parameter (in this case 
% the parameter c) that makes the output of a function zero. The function here is called Ffn.
% It takes in an initial guess for c and works out the corresponding value of F_total.
% It also takes in several other parameters (materialProps etc) that it needs in order to 
% calculate F_total.  Fzero then adjusts c until it finds a value such that Ffn
% returns zero for F_total (well, zero, to within the built-in tolerance - which
% can be changed if desired).
% See help fzero for details of its syntax.
% All the calculations within the function are basically your original calculations, 
% except that I've inserted calculations to extract material properties etc so the
% subsequent calculations use the same variable names that you created.
% 
% fzero returns the value of c that makes F_total zero.
% Subsequent calculations are again yours.
%
% I also changed your plot command slightly to plot points rather than lines.
%
% I hope this helps.

Said Bolluk 2020-9-12

在 MATLAB Online 中打开

Well, wise move. I studied a bit more about this fzero function, and understood what you did in the initial code. But, what is the [ ] when you call for the Ffn with the other parameters?

c = fzero(@Ffn, c0,[],ec_top, ecu_idx,materialProps, Layers, alphabeta);

Also, this approach is not valid for every strain value between zero and 4%. It works only for some particular strains -the ones in ecu_list- with their peculiar coefficients (alpha and beta). The most proper way is to integrate over the c for finding the actual strain value ec_i. Consequently, the after product Fc will be a function of c. However, there are two different equations for fc_i, being a function of ec_i. But, there might be two solutions since the results are too close. For example, let say that ec_i = 0.003. The result of the first equation is about 27 MPa, while that of the second is 24 MPa. So, how can this code decide for the right c value? I will upload a photo of these equations. Is there any way to make the code possible on finding the strain ec_i out of the c value, and then the real stress fc_i?

syms y;
ec_i = ec_top * (c - y) / c;
fc_i = fcu * ((2 * ec_i / ec0) - (ec_i / ec0) ^ 2); %for 0<ec_i<0.002
%%% or
fc_i = fcu * (1 - 0.15 * ((ec_i - ec0) / (ecu - ec0))); %for 0.002<ec_i<ecu
%%% then calculate the force
Fc = int(fc_i * b, y, 0, c) * 1e-3;

Thank you so much again for all the help. I hope I'm not asking too much.

Alan Stevens 2020-9-13

I guess so. I'm unfamiliar with the physics of this, so you will need to decide yourself!

Said Bolluk 2020-9-13

Thank you so much again for these great efforts. Your suggestions will answer the purpose for sure. I appreciate that.

请先登录，再进行评论。

How to speed up the iteration of the for loop and make some further improvement on the below code?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

10 个评论
显示 8更早的评论隐藏 8更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

How to speed up the iteration of the for loop and make some further improvement on the below code?

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

10 个评论 显示 8更早的评论隐藏 8更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

10 个评论
显示 8更早的评论隐藏 8更早的评论