# Trying to solve a system of equations with 14 equations and 14 unknowns

2 views (last 30 days)
Noah Adams on 19 Nov 2019
Commented: Walter Roberson on 1 Dec 2019
Below is my code for my system. I've tried a couple different ways and this is the way im on now. I get the error below when I try to run it. Ive made adjustments using ".^" and get the same error. I feel like my error is syntax at this point but if there is an easier way to solve this I'm open to ideas. In the end I am looking to solve for D. Thank you for the help!
Undefined operator '.^' for input arguments of type
'function_handle'.
Error in ^ (line 21)
Z = X.^Y;
Error in Group_Project_Test2 (line 30)
e5 = P_pmp==P_dim*rho*(pmp_spd^3)*D^5;
z1=30.48;
z2=45.72;
OL=1287.48;
Lx=OL/10;
L45=Lx*cos(pi/4);
Ly=(z2-z1);
Vfr_dim=.0364;
P_dim=.2758;
Vfr=(1.5*10^6)*(1/264.17)/(8*3600);
Kl_45=.4;
Kl_vlv=.4;
Kl_90=.2;
Kl_ent=.5;
Kl_ext=2;
g=9.81;
u=.001307;
rho=999.7;
Mfr=Vfr*rho;
x0=0;
syms D A V Re f P_pmp w_pmp Hl_min_45 Hl_min_vlv Hl_min_90 Hl_maj_45 ...
Hl_maj_x Hl_maj_y Hl_tot;
pmp_spd=@(D) Vfr*(Vfr_dim * D^3);
e1 = A==(pi/4)*D^2;
e2 = V==Vfr*(A)^-1;
e3 = Re==(rho*D*V)/u;
e4 = f==64*(Re^-1);
e5 = P_pmp==P_dim*rho*(pmp_spd^3)*D^5;
e6 = w_pmp==P_pmp*(Mfr^-1);
e7 = Hl_min_45==(Kl_45*(V^2))/((2*g));
e8 = Hl_min_vlv==(Kl_vlv*(V^2))/((2*g));
e9 = Hl_min_90==(Kl_90*(V^2))/((2*g));
e10 = Hl_maj_45==f*(L45/D)*((V^2)/(2*g));
e11 = Hl_maj_x==f*(Lx/D)*((V^2)/(2*g));
e12 = Hl_maj_y==f*(Ly/D)*((V^2)/(2*g));
e13 = Hl_tot==(8*Hl_min_45)+(4*Hl_min_vlv)+(2*Hl_min_90)+(4*Hl_maj_45)+(6*Hl_maj_x)+Hl_maj_y;
e14= z1 - z2 + w_pmp - Hl_tot ==0;
S=solve([e1,e2,e3,e4,e5,e6,e7,e8,e9,e10,e11,e12,e13,e14]...
,[D,A,V,Re,f,P_pmp,w_pmp,Hl_min_45,Hl_min_vlv,Hl_min_90,Hl_maj_45, ...
Hl_maj_x,Hl_maj_y,Hl_tot]);

Walter Roberson on 20 Nov 2019
pmp_spd=@(D) Vfr*(Vfr_dim * D^3);
That creates a function handle. The only things you can do with a function handle are:
• display it
• create a character vector version of it
• copy it (such as passing it to a function)
• invoke it (typically passing in parameters)
• use functions() to inquire about it
• ask how many inputs or outputs it has with nargin() and nargout()
e5 = P_pmp==P_dim*rho*(pmp_spd^3)*D^5;
but there you try to raise the function handle to the third power. The ^3 operation is not any of the above permitted operations.
Perhaps you want
e5 = P_pmp==P_dim*rho*(pmp_spd(D)^3)*D^5;
However! solve() is for requesting the complete set of indefinitely precise exact closed-form solutions whenever those can be computed. In any case in which you have floating point numbers in your calculation, using solve() is almost certainly the wrong thing to do. For example look at your line
g=9.81;
Why are you asking for an infinitely precise solution to your equations when you know that 9.81 is just an approximation, with g being 9.789 at the equator and 9.832 at the poles? What is the point of doing computations to a minimum of 117 digits when your inputs are only valid to 2 digits? Especially when the computation involves the roots of a polynomial of degree 18, 10191242662004313*Pi^2*z^18 - 300388392382226505177600*Pi^2*z^4 - 17069893026094080000*Pi*sqrt(2) - 52219967275837440000*Pi - 3249025000000000000000 (this can be reduced to a polynomial of degree 9)
You are doing a scientific calculation. You should be doing appropriate error analysis. How do you justify OL=1287.48 being 6 digits when your 3 digit g is not even correct?
vpasolve(). And not to very many digits either. And when vpasolve() complains about numeric instability, reflect upon the fact that your inputs are potentially not precise enough to justify any output.
##### 2 CommentsShowHide 1 older comment
Walter Roberson on 20 Nov 2019
Make the change I suggested for the function handle ; also use vpasolve() instead of solve()
You will need to weed through the 18 solutions.
Perhaps you can constrain the solutions? Perhaps some of the variables are known to be real valued or even nonnegative?

Alex Sha on 1 Dec 2019
One result I got:
a: 9.16543127083558
d: -3.41610736637613
v: 0.0215111341146459
re: -56206.8082673153
f: -0.00113861533204565
p_pmp: 3003.68991741061
w_pmp: 15.239484009305
hl_min_45: -2.79268550704629E-5
hl_min_vlv: -9.15501741224959E-6
hl_min_90: -4.41086849672573E-6
hl_maj_45: -1.78570364039789E-5
hl_maj_x: -2.69928669675165E-5
hl_maj_y: -4.40779091820182E-6
hl_tot: -0.000511322453475079
Walter Roberson on 1 Dec 2019
Slightly more precisely:
A = 9.165487319777216
D = 3.416117811538274
Hl_tot = 0.0001316926473876508
P_pmp = 3003.818495417170
Re = 56206.51501771111
V = 0.02151095611165516
f = 0.001138658035991612
w_pmp = 15.24013169264739
Hl_maj_45 = 7.156594239983986E-7
Hl_maj_x = 1.012095263458652E-6
Hl_maj_y = 1.198024964668178E-7
Hl_min_45 = 9.433664278033732E-6
Hl_min_90 = 4.716832139016866E-6
Hl_min_vlv = 9.433664278033732E-6
These differ from Alex's solutions starting in the 10th decimal place. These were produced by a symbolic mathematics package operating at a higher number of decimal places and then reduced down to 16 decimals.
This is the only real-valued solution with positive D. Mathematically there are another 17 solutions to the equations, but if we take the fact that D represents a pipe diameter to impose a constraint that D must be real valued and positive then this is the only solution.