MATLAB: How to update the variables of the function after solving for a symbolic expression

rank-deficient systemsolve

function [curvature, moment] = moment_curvature
    curvature_list = zeros(1, 1000000);
    moment_list = zeros(1, 1000000);
    h = 600;
    d = 570;
    b = 300;
    As = 800;
    fcd = 30 / 1.5;
    Ec = 12680 + 460 * fcd;
    Es = 200 * 1e3;
    ec0 = 2 * fcd / Ec;
    ecu = 0.0038;
    for ec_top = 0:5e-4:ecu
        syms c y;
        %%% tension zone: contribution from the bottom steel layers
        fs = Es * (c - d) * ec_top / c;
        Fs = fs * As * 1e-3;
        %%% compression zone: contribution from the upper concrete layers
        ec_i = ec_top * y / c;
        fc_i = fcd * ((2 * ec_i / ec0) - (ec_i / ec0) ^ 2);
        Fc = int(fc_i * b, y, 0, c) * 1e-3;
        F_total = Fs + Fc == 0;
        %%% the net force must be zero, to find c:
        c = solve(F_total,c,'Real',true); %%% I need the positive real root
        %%% after finding c, reassign the variables
        fs = subs(fs);
        Fs = subs(Fs);
        ec_i = subs(ec_i);
        fc_i = subs(fc_i);
        Fc = subs(Fc);
        %%% calculating the M and K values after finding c
        moment = Fs * (d - h / 2) * 1e-3 + int((fc_i * b * (h/2 - y)), y, 0, c) * 1e-3;
        curvature = ec_top / c;
        %%%
        curvature_list = [curvature_list, curvature];
        moment_list = [moment_list, moment];
        %%% plot the moment-curvature curve of the beam section
        plot(curvature_list, moment_list);
    end
end
Hello,
In this code, I am trying to plot the Moment-Curvature Diagram of a reinforced concrete beam section. The for loop is there to find the values of the diagram for the specific strain values. For each value, there will be different stresses, and consequently forces and moments. All these values are based on c, which is the neutral axis depth of the section. This depth determines the compression (upper) and tension (lower) zones of the beam. After assigning the values with c, I want to solve the final equation of F_total, which is the total force acting on the cross-section. It must be zero if there is no external normal force so that the c value should be defined according to this.
However, I have many errors in the code. To start with, I need to figure out how to find the real positive root, which gives the actual neutral axis depth c, of the function of F_total. When I run the code, I see this warning: Solution is not unique because the system is rank-deficient. I will attach the full error screen. I will be pleased if I could get any help.
Thank you,

Best Answer

Eventually you get two real roots for c. Your line dividing by c then has a problem because you used / rather than ./ and your attempt to append values onto the list has problems because you are appending two values in a column to the end of a list.
Related Question