Getting a wrong answer using solve()

Matlab gives wrong solutions for a system of two equations. The code is given below,

syms pa p1a;
c=1.056;
Eq1= -8 * p1a ^ 2 + (12 * c - 10 * pa - 6) * p1a - 2 * c ^ 2 + 3 * pa ^ 2 + 4 * c - 4 * pa - 1;
Eq2= -2 * pa ^ 2 + (6 * c - 8 * p1a - 3) * pa - 2 * c ^ 2 + 8 * p1a ^ 2 + 4 * c - 6 * p1a - 1;
[solutions_p1a, solutions_pa] =...
solve(Eq1==0,Eq2==0,p1a,pa);
solutions = [solutions_p1a, solutions_pa];

Matlab gives the following answer:

solutions =
 [                                           7.9663553729228043315977932331051,                                     0.050975475760590065640643393645407]
[                                         -0.19156088730325461282546227710482,                                      -1.1108531940976453676647830585577]
[ 0.086183245911339552040734591024871 + 0.012548420079599710197932923828486*i, 1.2388776346787317326447228936806 + 1.0340545303434847217578152453505*i]
[ 0.086183245911339552040734591024871 - 0.012548420079599710197932923828486*i, 1.2388776346787317326447228936806 - 1.0340545303434847217578152453505*i]
 

Substituting pa= 0.050975475760590065640643393645407 and p1a=7.9663553729228043315977932331051 in the equations gives: Eq1=-457.81428633120934568565749111782 and Eq2=457.81428633120934568565749111782.

Clearly, the answer is wrong. Is there anyone who could tell me why? Really appreciated!

Answer 1

[Alan Stevens]

Running your code as is (I copied and pasted, with no changes), I get

>> syms pa p1a;
c=1.056;
Eq1= -8 * p1a ^ 2 + (12 * c - 10 * pa - 6) * p1a - 2 * c ^ 2 + 3 * pa ^ 2 + 4 * c - 4 * pa - 1;
Eq2= -2 * pa ^ 2 + (6 * c - 8 * p1a - 3) * pa - 2 * c ^ 2 + 8 * p1a ^ 2 + 4 * c - 6 * p1a - 1;
[solutions_p1a, solutions_pa] =...
solve(Eq1==0,Eq2==0,p1a,pa);
solutions = [solutions_p1a, solutions_pa]
 
solutions =
 
[ (765625*root(z^4 - (17369*z^3)/12250 - (60222*z^2)/765625 + (6110509*z)/19140625 - 122259598/2392578125, z, 1)^3)/276196 - (6341875*root(z^4 - (17369*z^3)/12250 - (60222*z^2)/765625 + (6110509*z)/19140625 - 122259598/2392578125, z, 1)^2)/1657176 - (253015*root(z^4 - (17369*z^3)/12250 - (60222*z^2)/765625 + (6110509*z)/19140625 - 122259598/2392578125, z, 1))/828588 + 86656187/103573500, root(z^4 - (17369*z^3)/12250 - (60222*z^2)/765625 + (6110509*z)/19140625 - 122259598/2392578125, z, 1)]
[ (765625*root(z^4 - (17369*z^3)/12250 - (60222*z^2)/765625 + (6110509*z)/19140625 - 122259598/2392578125, z, 2)^3)/276196 - (6341875*root(z^4 - (17369*z^3)/12250 - (60222*z^2)/765625 + (6110509*z)/19140625 - 122259598/2392578125, z, 2)^2)/1657176 - (253015*root(z^4 - (17369*z^3)/12250 - (60222*z^2)/765625 + (6110509*z)/19140625 - 122259598/2392578125, z, 2))/828588 + 86656187/103573500, root(z^4 - (17369*z^3)/12250 - (60222*z^2)/765625 + (6110509*z)/19140625 - 122259598/2392578125, z, 2)]
[ (765625*root(z^4 - (17369*z^3)/12250 - (60222*z^2)/765625 + (6110509*z)/19140625 - 122259598/2392578125, z, 3)^3)/276196 - (6341875*root(z^4 - (17369*z^3)/12250 - (60222*z^2)/765625 + (6110509*z)/19140625 - 122259598/2392578125, z, 3)^2)/1657176 - (253015*root(z^4 - (17369*z^3)/12250 - (60222*z^2)/765625 + (6110509*z)/19140625 - 122259598/2392578125, z, 3))/828588 + 86656187/103573500, root(z^4 - (17369*z^3)/12250 - (60222*z^2)/765625 + (6110509*z)/19140625 - 122259598/2392578125, z, 3)]
[ (765625*root(z^4 - (17369*z^3)/12250 - (60222*z^2)/765625 + (6110509*z)/19140625 - 122259598/2392578125, z, 4)^3)/276196 - (6341875*root(z^4 - (17369*z^3)/12250 - (60222*z^2)/765625 + (6110509*z)/19140625 - 122259598/2392578125, z, 4)^2)/1657176 - (253015*root(z^4 - (17369*z^3)/12250 - (60222*z^2)/765625 + (6110509*z)/19140625 - 122259598/2392578125, z, 4))/828588 + 86656187/103573500, root(z^4 - (17369*z^3)/12250 - (60222*z^2)/765625 + (6110509*z)/19140625 - 122259598/2392578125, z, 4)]
 

Using vpasolve, I get

>> syms pa p1a;
c=1.056;
Eq1= -8 * p1a ^ 2 + (12 * c - 10 * pa - 6) * p1a - 2 * c ^ 2 + 3 * pa ^ 2 + 4 * c - 4 * pa - 1;
Eq2= -2 * pa ^ 2 + (6 * c - 8 * p1a - 3) * pa - 2 * c ^ 2 + 8 * p1a ^ 2 + 4 * c - 6 * p1a - 1;
[solutions_p1a, solutions_pa] =...
vpasolve(Eq1==0,Eq2==0,p1a,pa);
solutions = [solutions_p1a, solutions_pa]
 
solutions =
 
[  0.29377663793901891172446517857323,  0.39303774024642807767205311835288]
[  0.64327467918062628134033133387723,  0.20094137279799901439637447288045]
[ -0.26835932140677888831670790351722, -0.49177389771394883707711632108608]
[  0.12363453489937859321109506453615,   1.3156723356899299082739948523017]

This was with MATLAB R2018b.

[Walter Roberson]

I get the same as Alan in R2020a.

[sara s]

I upgraded my matlab. I really appreciate your help and time.

[John D'Errico]

It probably was not the upgrade that fixed it, but just making sure that you had the correct system of equations and constants, etc. Very often, when we see someone thinking they had the wrong solution, they had just changed some constant, or made some subtle fix to their equations. So it is the correct solution, but a solution to a completely different problem than they thought they solved.