"Error using symengine Unable to convert expression into double array" when using solve function

Hi,

I'm trying to set up a symbolic solver to help solve an equation I can't solve analytically. I need the output of the code in a usable format (double) so that I can use it in another equation. However, when I attempt to run the code, it throws the following error:

Error using symengine

Unable to convert expression into double array.

Error in sym/double (line 698)

Xstr = mupadmex('symobj::double', S.s, 0);

Error in script (line 8)

f_x(counter) = double(solve(eqn, x));

Does anyone have ideas for how to remedy this?

Code for reference:

%% Constants
g = 1.2; y = 0.1; counterlim = 50
while counter < counterlim
    syms x
    eqn = y ==  c * sqrt((x ^ (2 / g) - x ^ ((g + 1) / g)));
    f_x(counter) = double(solve(eqn, x));
    c = c + 0.01
    counter = counter + 1;
end

Answer 1

[Walter Roberson]

%% Constants
C_v = 0.5;
A_o = 1;
P_upstream = 1;
MW = 0.004;
Z = 0.95;
R = 8.314;
T1 = 100;
gamma = 1.660;
mdot_desired = 0.1;
counterlim = 50;

f_Pratio = zeros(1,counterlim);
for counter = 1 : counterlim
    syms Pratio
    eqn = mdot_desired == C_v * A_o * P_upstream * sqrt((2 * MW / ...
        (Z * R * T1)) * (gamma / (gamma - 1)) * (Pratio ^ (2 / gamma) - ...
        Pratio ^ ((gamma + 1) / gamma)));
    sol = vpasolve(eqn, Pratio);
    f_Pratio(counter) = double(sol);
    A_o = A_o + 0.01;
    counter = counter + 1;
end

Note: the answers are all complex valued. In each case, the imaginary component is roughly 2*counter .

There are two solutions for each equation; the solutions are complex conjugates of each other.

[Walter Roberson]

For the revised question, and selecting down to real values:

g = 1.2; y = 0.1; counterlim = 50;
c0 = 1;    %CHECK THIS
delta_c = 0.01;

syms c x real
eqn = y ==  c * sqrt((x ^ (2 / g) - x ^ ((g + 1) / g)));
sol = solve(eqn, x, 'returnconditions', true);
cvals = c0 + (0:counterlim-1) * delta_c;
eqns = subs(sol.conditions, c, cvals);
fx_cell = arrayfun(@(eqn) double(subs(sol.x, sol.parameters, solve(eqn))).', eqns, 'uniform', 0);
fx = vertcat(fx_cell{:});
plot(cvals, fx)

There are two real solutions for each c value, at least in the range being worked over.

Note that your revised question did not initialize c, so I had to speculate that your c was taking the place of your A0 of your original question, which you had initialized to 1. If that is not correct, then adjust the c0 variable.