# Solving a system of nonlinear equations with bounded variables

3 views (last 30 days)
Jonathan Butt on 3 Dec 2016
Answered: Alex Sha on 2 Feb 2020
Hello all,
I'm building a parametric model of a mass damper system. The majority of this work is simple enough in simulink however I have a complex damper which requires six nonlinear equations to be solved.
I have managed to use 'fsolve' and 'vpasolve' to do this however, when I try and apply bounds to the system to recreate the limits of the physical system the vpa system breaks down. This is problem number one.
One of the things that's really confusing me is that the my equation states that the sum of the flows is equal to the velocity multiplied by the area of by supporting rod. However when I sum these flows together after the solver has completed the value is a factor of 10 out from the value given in the equation.
I've moved to vpasolve because I could constrain/apply bound values in fsolve (not that i'm having much luck at doing this with vpasolve).
Code attached below, any help, mild abuse for lack of programming prowess or general feedback would be greatly appreciated. Cheers Jon
if true
%--------------------jonBoy's code------------------------
%constant system parameters
A_rod = 0.010;
A_b = 0.020;
A_0 = 0.022;
A_v = 0.020;
a = 0.003;
b = 0.005;
Cd = 0.60;
Cf = 0.30;
dens = 1000;
D_v = 0.10;
D_p = 0.110;
L = 0.010;
k = 5000;
Mu = 0.50;
F_sp = 5;
%system variables (FED IN FROM SIMULINK)
x_dot = 5;
Pc = 300000;
%Solving things
syms Q_b Q_v Q_lp y Pr Pv
[Q_b_sol, Q_v_sol, Q_lp_sol, y_sol, Pr_sol, Pv_sol] = vpasolve(...
[...
A_rod*x_dot == Q_b - Q_v - Q_lp ...
Q_b == (A_b * Cd * sqrt( (2*(Pc-Pr))/dens) ) , ...
Q_v == (A_0 * Cd * sqrt( (2*(Pc-Pv))/dens) ) , ...
Q_v == ( (a*pi*D_v*y)*Cd )*sqrt( (2*(Pv-Pr))/dens ) , ...
k*y == ( (Pv-Pr)*A_v ) - ( dens * Cf * ((Q_v/A_0)^2) ) + F_sp , ...
Q_lp ==( ( ( ((Pc-Pr)*b^3)/(12*Mu*L) ) + ( x_dot*(b/2) ) )*pi*D_p ), ...
],...
[Q_b Q_v Q_lp y Pr Pv])
%[-inf inf, -inf inf, -inf inf, 0 0.00002, -inf inf, -inf inf] )
%Display useful facts.....
disp(A_rod*x_dot);
flow = Q_b_sol+Q_v_sol+Q_lp_sol;
disp(flow);
disp(' all done ');
end

Torsten on 5 Dec 2016
Edited: Torsten on 5 Dec 2016
Didn't you forget to set a comma after the equation
A_rod*x_dot == Q_b - Q_v - Q_lp
?
Best wishes
Torsten.

Alex Sha on 2 Feb 2020
There are two solutions:
1:
q_b: 0.305291000916417
q_v: 0.0179815487538954
q_lp: 0.237309452162523
y: 1.25168611470519
pr: -23620.1223630127
pv: 299072.153077397
Fevl:
1.42941214420489E-15
5.55111512312578E-17
3.88578058618805E-16
4.75314232417645E-16
1.90993887372315E-11
-5.55111512312578E-17
2:
q_b: 0.0650963672810365
q_v: 0.000183595254128151
q_lp: 0.0149127720269083
y: 0.0598501155302159
pr: 285286.329745876
pv: 299999.903273596
Fevl:
-4.16333634234434E-17
-8.18789480661053E-16
2.17154311026035E-13
3.85393214737129E-15
-1.16529008664656E-11
-2.55004350968591E-16