Error and Warning solving equations.

15 次查看（过去 30 天）
Tamara Laranga Barreiro
编辑： John D'Errico 2022-12-19
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
  3 个评论
显示 1更早的评论
Steven Lord
Steven Lord 2022-12-19
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
Tamara Laranga Barreiro 2022-12-19
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);
%}

请先登录，再进行评论。

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

采纳的回答

John D'Errico
John D'Errico 2022-12-19
编辑：John D'Errico 2022-12-19
Ran in:
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.

更多回答（1 个）

Torsten
Torsten 2022-12-19
Ran in:
%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

类别

在 Help Center 和 File Exchange 中查找有关 MuPAD 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by