Error and Warning solving equations.

Hi, I'm trying to solve some equations with Matlab for a project. The next script runs perfectly (I'm only adding the troubuling part of the program, I have defined before the values of all variables);

%WORKS
syms rho_h; 
Pr_oh = rho_h * S_r * (Vtip_r)^(3) * (sigma_r * Cd0)/8;
Pr_ih = Wo * (sqrt(Wo/(2 * rho_h *S_r)));
Pd_h = Pmc * rho_h/rho;
E1 = ((1+kt) * (Pr_ih + Pr_oh)) - Pd_h;
R1 = double(solve(E1, rho_h));
h = (288.15/(-6.5 * 10^(-3))) * ((R1/rho)^(1/4.252) - 1);
disp('flight roof: ');
disp(h);

When I only change the value of Pr_ih adding Wo*0.3

Pr_ih = Wo * (sqrt(Wo/(2 * rho_h *S_r)) + 0.3);

I the following outout. I already tryed everything and don't understand why this keeps happening;

Warning: Solutions are parameterized by the symbols: z1. To include parameters and conditions in the solution,
specify the 'ReturnConditions' value as 'true'. 
> In sym/solve>warnIfParams (line 475)
In sym/solve (line 357)
In HS_VPF (line 31) 
Warning: Solutions are only valid under certain conditions. To include parameters and conditions in the solution,
specify the 'ReturnConditions' value as 'true'. 
> In sym/solve>warnIfParams (line 478)
In sym/solve (line 357)
In HS_VPF (line 31) 
Error using symengine
Unable to convert expression containing symbolic variables into double array. Apply 'subs' function first to
substitute values for variables.

Error in sym/double (line 872)
        Xstr = mupadmex('symobj::double', S.s, 0);

Error in HS_VPF (line 31)
R1 = double(solve(E1, rho_h));

Related documentation

[John D'Errico]

If you really want help, you would need to provide sufficient information that we can actually run that code. There are multiple parameters undefined, so all we can do is look at a complicated expression, and try to guess what is happening.

For example, what are S_r, Vtip_r, Sigma_r, Wo, rho, kt, Pmc? Maybe some others.

Why make it harder for you to obtain help?

[Steven Lord]

What happens when you add 'ReturnConditions', true to your solve call (which will require separating your solve call and the conversion to double into two statements, as you'll need to call solve with multiple outputs) as per the "Use Parameters and Conditions to Refine Solution" example on the solve documentation page?

[Tamara Laranga Barreiro]

I tried what you suggested Steven Lord but it doesnt work either. I'm going to follow John's advice and post the full code.

This first part its only the declaration of variables;

%General Data

Wo = 3584.74 * 9.8;
Vm = 268.76/3.6; 
Rm = 680 * 10^3;
rho = 1.225;
l_tr = 6.71;
Pmc = 958.98 * 10^3;
Tmc = 751.5 * 10^3;

%NACA0012 Polars
Cl = 5.74;
Cd0 = 0.0085;
Cd2 = 0.263;

%Main Rotor
b_r = 4;
D_r = 11.5;
R_r = D_r/2; 
c_r = 0.33;
sigma_r = 0.073309651;
omega_r = 37.01;
S_r = pi * R_r^2;
Vtip_r = omega_r * R_r;
Tu_r = rho * S_r * (Vtip_r)^2;
Pu_r = rho * S_r * (Vtip_r)^3;
Vi0 = sqrt(Wo/(2*rho*S_r));

%Torque rotor
b_t = 2;
D_t = 2.21;
R_t = D_t/2;
c_t = 0.2214;
sigma_t = 0.155749906;
omega_t = 188.8;
S_t = pi * R_t^2;
Vtip_t = omega_t * R_t;
Tu_t = rho * S_t * (Vtip_t)^2;
Pu_t = rho * S_t * (Vtip_t)^3;

This is the one that calculates the flight roof;

%Power VPF sea level
%Main Rotor
Pr_i = Wo*Vi0;
Pr_o = (sigma_r*Cd0/8)*Pu_r;
Pr_VPF_sl = Pr_o + Pr_i;

%Torque Rotor
Tt_sl = Pr_VPF_sl/(l_tr*omega_r);
CTt_sl = Tt_sl/Tu_t;
lambda_t_sl = sqrt(CTt_sl/2);
Pt_i = CTt_sl*lambda_t_sl*Pu_t;
Pt_o = (sigma_t*Cd0/8)*Pu_t;
Pt_VPF_sl = Pt_i + Pt_o;
%Constant ratio between Main and torque rotor 
kt = Pt_VPF_sl/Pr_VPF_sl;

%FLIGHT ROOF

syms rho_h; 
Pr_oh = rho_h * S_r * (Vtip_r)^(3) * (sigma_r * Cd0)/8;
Pr_ih = Wo * (sqrt(Wo/(2 * rho_h *S_r)));
Pd_h = Pmc * rho_h/rho;
E1(rho_h) = ((1+kt) * (Pr_ih + Pr_oh)) - Pd_h;
R1 = solve(E1, rho_h, true);
R1_1 = double(solve(R1, rho_h));
h = (288.15/(-6.5 * 10^(-3))) * ((R1_1/rho)^(1/4.252) - 1);
disp('FLIGHT ROOF: ');
disp(h);

%{
%Service flight roof
% Vz= 0.3
syms rho_hs;
Pr_ohs = rho_hs * S_r * (Vtip_r)^3 * (sigma_r * Cd0)/8;
Pr_ihs = Wo * (sqrt(Wo/(2 * rho_hs * S_r)) + 0.3);
Pd_hs = Pmc * rho_hs/rho;
E2 = ((1 + kt) * (Pr_ihs + Pr_ohs)) - Pd_hs;
R2 = double(solve(E2, rho_hs));
hs = (288.15/(-6.5 * 10^(-3))) * ((R2/rho_hs)^(1/4.252) - 1);
disp('Service flight roof is: ');
disp(hs);
%}

Answer 1

[John D'Errico]

Admittedly, it is a mess of stuff to paste in. :)

Wo = 3584.74 * 9.8;
Vm = 268.76/3.6; 
Rm = 680 * 10^3;
rho = 1.225;
l_tr = 6.71;
Pmc = 958.98 * 10^3;
Tmc = 751.5 * 10^3;
%NACA0012 Polars
Cl = 5.74;
Cd0 = 0.0085;
Cd2 = 0.263;
%Main Rotor
b_r = 4;
D_r = 11.5;
R_r = D_r/2; 
c_r = 0.33;
sigma_r = 0.073309651;
omega_r = 37.01;
S_r = pi * R_r^2;
Vtip_r = omega_r * R_r;
Tu_r = rho * S_r * (Vtip_r)^2;
Pu_r = rho * S_r * (Vtip_r)^3;
Vi0 = sqrt(Wo/(2*rho*S_r));
%Torque rotor
b_t = 2;
D_t = 2.21;
R_t = D_t/2;
c_t = 0.2214;
sigma_t = 0.155749906;
omega_t = 188.8;
S_t = pi * R_t^2;
Vtip_t = omega_t * R_t;
Tu_t = rho * S_t * (Vtip_t)^2;
Pu_t = rho * S_t * (Vtip_t)^3;
% This is the one that calculates the flight roof; 
%Power VPF sea level
%Main Rotor
Pr_i = Wo*Vi0;
Pr_o = (sigma_r*Cd0/8)*Pu_r;
Pr_VPF_sl = Pr_o + Pr_i;
%Torque Rotor
Tt_sl = Pr_VPF_sl/(l_tr*omega_r);
CTt_sl = Tt_sl/Tu_t;
lambda_t_sl = sqrt(CTt_sl/2);
Pt_i = CTt_sl*lambda_t_sl*Pu_t;
Pt_o = (sigma_t*Cd0/8)*Pu_t;
Pt_VPF_sl = Pt_i + Pt_o;
%Constant ratio between Main and torque rotor 
kt = Pt_VPF_sl/Pr_VPF_sl;
%FLIGHT ROOF
syms rho_h; 
Pr_oh = rho_h * S_r * (Vtip_r)^(3) * (sigma_r * Cd0)/8;
Pr_ih = Wo * (sqrt(Wo/(2 * rho_h *S_r)));
Pd_h = Pmc * rho_h/rho;
E1(rho_h) = ((1+kt) * (Pr_ih + Pr_oh)) - Pd_h;

Now, look at the result, in the first form.

vpa(E1,5)

ans(rho_h) = 

And that should indeed have a simple solution if we set it to zero.

syms rho_hs;
Pr_ohs = rho_hs * S_r * (Vtip_r)^3 * (sigma_r * Cd0)/8;
Pr_ihs = Wo * (sqrt(Wo/(2 * rho_hs * S_r)) + 0.3);
Pd_hs = Pmc * rho_hs/rho;
E2 = ((1 + kt) * (Pr_ihs + Pr_ohs)) - Pd_hs;

Now consider the last part.

vpa(E2,5)

ans = 

This is essentially equivalent to a cubic polynomial in the unknown. If we just try to throw it at solve, solve has some issues.

rho_sol = solve(E2,'returnconditions',true)

rho_sol = struct with fields:
        rho_hs: z1
    parameters: z1
    conditions: (-signIm(- (z1*pi^(1/2)*7914398259159293173485080988876800i)/381202755794727460222990550282303 + (z1*pi^(3/2)*3932748271564133464119341266544943i)/5566642320511124857678584929583104 + (pi^(1/2)*3i)/10) == 1 | (7914398259159293173485080…

As you can see, solve did not actually solve the problem. It just says the solution is in the form of a polynomial in z1, where z1 is pretty much the solution to your original problem.

First, I'll plot this function.

fplot(E2)

It appears there is a single real root.

R2 = vpasolve(E2)

R2 = 

0.80019627293622902959538461016129

as found by vpasolve. Do you really need to see an algebraic solution anyway? It would be a complete mess.

Answer 2

[Torsten]

%General Data
Wo = 3584.74 * 9.8;
Vm = 268.76/3.6; 
Rm = 680 * 10^3;
rho = 1.225;
l_tr = 6.71;
Pmc = 958.98 * 10^3;
Tmc = 751.5 * 10^3;
%NACA0012 Polars
Cl = 5.74;
Cd0 = 0.0085;
Cd2 = 0.263;
%Main Rotor
b_r = 4;
D_r = 11.5;
R_r = D_r/2; 
c_r = 0.33;
sigma_r = 0.073309651;
omega_r = 37.01;
S_r = pi * R_r^2;
Vtip_r = omega_r * R_r;
Tu_r = rho * S_r * (Vtip_r)^2;
Pu_r = rho * S_r * (Vtip_r)^3;
Vi0 = sqrt(Wo/(2*rho*S_r));
%Torque rotor
b_t = 2;
D_t = 2.21;
R_t = D_t/2;
c_t = 0.2214;
sigma_t = 0.155749906;
omega_t = 188.8;
S_t = pi * R_t^2;
Vtip_t = omega_t * R_t;
Tu_t = rho * S_t * (Vtip_t)^2;
Pu_t = rho * S_t * (Vtip_t)^3;
%Power VPF sea level
%Main Rotor
Pr_i = Wo*Vi0;
Pr_o = (sigma_r*Cd0/8)*Pu_r;
Pr_VPF_sl = Pr_o + Pr_i;
%Torque Rotor
Tt_sl = Pr_VPF_sl/(l_tr*omega_r);
CTt_sl = Tt_sl/Tu_t;
lambda_t_sl = sqrt(CTt_sl/2);
Pt_i = CTt_sl*lambda_t_sl*Pu_t;
Pt_o = (sigma_t*Cd0/8)*Pu_t;
Pt_VPF_sl = Pt_i + Pt_o;
%Constant ratio between Main and torque rotor 
kt = Pt_VPF_sl/Pr_VPF_sl;
%FLIGHT ROOF
syms rho_h real
Pr_oh = rho_h * S_r * (Vtip_r)^(3) * (sigma_r * Cd0)/8;
Pr_ih = Wo * (sqrt(Wo/(2 * rho_h *S_r)));
Pd_h = Pmc * rho_h/rho;
E1 = ((1+kt) * (Pr_ih + Pr_oh)) - Pd_h

E1 = 

E1 = simplify(expand(E1*sqrt(rho_h)))

E1 = 

format long
R1 = double(solve(E1,rho_h))

R1 = 
   0.789373365893613

h = (288.15/(-6.5 * 10^(-3))) * ((R1/rho)^(1/4.252) - 1);
disp('FLIGHT ROOF: ');

FLIGHT ROOF: 

disp(h);

     4.352898824088797e+03