How to Solve a System of Equations with Multiple Variables

I'm working through a propulsion homework problem. I've gotten to the end where I've ended up with four equations and four unknowns in a fairly complex system of equations:

2=x+y;
2=2x+y+2z+w;
0.046328=(w^2/z)*(2/(x+y+z+w));
.65116=((y*z^.5)/x)*(2/(x+y+z+w))^.5

Obviously this is a difficult system to solve by hand. How do I write a code in Matlab to give me the values of x, y, z, and w?

[jgg]

To make it easier, you should first eliminate the first equation from all the other ones, getting rid of (probably) y.

Answer 1

[jgg]

You can do this using numerical optimization by setting up an equation in which the objective function is the norm of the equations. For instance, if you want to solve:

x + y = 2
2x + y = 3

You can consider the function F which evaluates:

F(1) = abs(x + y - 2)
F(2) = abs(2x + y - 3)

A solution to the original system of equations would also be a solution such that F = 0. You can implement this using any solver you'd like in Matlab. For instance, I eliminated the first equation from your problem, then set up this objective function:

 function [ out ] = func(X)

 x = X(1);
 z = X(2);
 w = X(3);

 C1 = 0.046328;
 C2 = 0.65116;

 m = (2/(2+z+1));

 f1 = x + 2*z + w;
 f2 = ((w^2)/z)*m - C1;
 f3 = ((2-x)/(x))*sqrt(z)*sqrt(m) - C2;

 out = abs(f1) + abs(f2) + abs(f3);
 end

You can then call it from the command line, or in a script with:

 options = optimset('MaxFunEvals', 100000);
 fminsearch(@func,[1,1,1],options)

You'll need to check a few things though:

• I might have made a mistake simplifying or setting up the problem. Check it!
• You didn't mention any constraints, so I used a gradient-free unconstrained solver. If you have restrictions on your variables, you might want to take those into account by using a solver like fmincon.
• All of the built-in solvers are local solvers. I chose an arbitrary starting point (1,1,1) to estimate from, but that might not the global optimum. You'll want to either prove your function is globally concave, or try lots of starting points.
• You'll also probably want to play with the tolerances on the objective. It doesn't look that close to a solution from my guess here.
• Alternatively, you make use a global optimizer if you have the toolbox for it, such as ga.

[John D'Errico]

Optimization won't help when there are no real solutions.

Answer 2

[John D'Errico]

There appear to be no real solutions.

syms x y z w
E1 = 2==x+y;
E2 = 2==2*x + y + 2*z + w;
E3 = 0.046328 == (w^2/z)*(2/(x+y+z+w));
E4 = .65116 == ((y*z^.5)/x)*(2/(x+y+z+w))^.5;
res = solve(E1,E2,E3,E4);

vpa(res.x)
ans =
2.0 + 1.0868394327073291299733800977949i
2.0 - 1.0868394327073291299733800977949i
2.0 + 0.9941240974066114685606129952539i
2.0 - 0.9941240974066114685606129952539i
2.0 + 1.0868394327073291299733800977949i
2.0 - 1.0868394327073291299733800977949i
2.0 + 0.9941240974066114685606129952539i
2.0 - 0.9941240974066114685606129952539i

Answer 3

[Walter Roberson]

That system has 8 solutions. None of the solutions are real-valued in even one of the components (but y might possibly have a purely-imaginary solution). Are you expecting real-valued solutions?

Your system of equations is very sensitive to the precision of your 0.046328 and your .65116 . Are we to take it that .65116 is 65116/100000 exactly, or are we to understand it as 65116/100000 +/- 5/1000000 ? If we say that the "real" value is 65116/100000 + delta2 for some unknown delta2 in the range +/- 5/1000000, then the solution to the equations involves the root of an expression involving delta2^20 multiplied by large numbers.

I have done a bunch more investigating. If we take the values as being approximate, then it is possible to determine some of the conditions under which real-valued solutions can be found. If we call the error on the 0.046328 as "delta1" and the error on the .65116 as "delta2", then:

If we assume that delta2 is close to 0, then delta1 needs to be at least +5.8114 -- very large in comparison, making it unlikely that the .65116 is fundamentally correct.

If we assume that delta1 is close to 0, then delta2 needs to be at least +0.6595228062 or -0.7587797008 (actually that's a double root, so safer would be -2.177082208). Adding +0.6595228062 to .65116 is slightly more than a factor of 2, suggesting perhaps a relatively minor miscalculation in determining it.

[Justin Bruh]

The 0.046328 and 0.65116 values are pulled from a table in the book's appendix based on the parameters given in the problem. I've also been trying with an example from the book and I'm getting similar non-real solutions from Matlab. The example uses:

2=2x+2y+w
5=x+2Z+w
0.5816=((y*(z^.5))/x)*(1/(x+y+z+w))^.5
0.9570=((w*(y^.5))/x)*(1/(x+y+z+w))^.5

And it lists results:

x=0.271
y=0.213
z=1.849
w=1.032

Also, the answers are "number of moles" so they must be real and non-negative.

[Justin Bruh]

Sorry, that posted poorly. Equations:

 2=2x+2y+w; 
5=x+2Z+w; 
0.5816=((y*(z^.5))/x)*(1/(x+y+z+w))^.5;
0.9570=((w*(y^.5))/x)*(1/(x+y+z+w))^.5;

results:

 x=0.271 
y=0.213 
z=1.849 
w=1.032

[Walter Roberson]

In those equations you just posted, two of the equations include sub-expressions of the form (1/(x+y+z+w))^.5 . In the equations you first posted, only one of the equations included a sub-expression of that form. So besides some differences in numeric values and a couple of changes of variables to find corresponding terms, the pattern is just completely different for one of the equations. It looks to me as if perhaps a square root was not taken on

0.046328=(w^2/z)*(2/(x+y+z+w));