# Complicated system of two non linear equations in three variables

I have a system of two equations in three variables. The two expressions involved in the system are obtained diagonalizing the following matrix with Eigensystem[]:

m = {{M1, 0, -4.28345, 42.8361}, {0, M2, 7.9998, -80.001}, {-4.28345, 7.9998, 0, -mu}, {42.8361, -80.001, -mu, 0}},


where all the involved parameters are fixed except for mu, M1 and M2. If you proceed in the diagonalization you'll find some pretty heavy expressions for the eigenvalues e1, e2, e3 and e4.

Now, I need to solve the system:

e2 == 50 && e3 == 200


and obtain M1(mu) and, separately, M2(mu).

NSolve[] or FindIstance[] juststarted the evaluation and never end it. I also tried to approximate the eigenvalue expressions for mu -> inf or mu -> 0 in order to cover at least some of the range of mu, but the resulting expression was still a too complicated function of all three parameters and numerical algorithms kept not working.

I find it hard to believe that the system is unsolvable, so if anyone has any suggestion or, also, a computer powerful enough to make NSolve[] work, is very welcome.

• It looks like a good application for the Schur complement (also known in physics as the Feshbach map). Aug 12, 2023 at 18:01

If two of the eigenvalues are 50 and 200, and the characteristic polynomial of the matrix is computed in the new variable t, then it must be divisible by (t-50)*(t-200). That can be used to obtain a remainder polynomial in t, the coefficients of which must equate to zero. Solving that for the variables of interest is straightforward.

mat = {{m1, 0, -4.28345, 42.8361}, {0, m2,
7.9998, -80.001}, {-4.28345, 7.9998,
0, -mu}, {42.8361, -80.001, -mu, 0}};

cp = CharacteristicPolynomial[mat, t];
rootpoly = (t - 50)*(t - 200);
{quo, rem} = PolynomialQuotientRemainder[cp, rootpoly, t]

(* Out[310]= {44182.6 - 250. m1 - 250. m2 + 1. m1 m2 -
1. mu^2 + (250. - 1. m1 - 1. m2) t + 1. t^2, -4.41826*10^8 +
2.5*10^6 m1 + 2.5*10^6 m2 - 10000. m1 m2 + 1279.98 m1 mu +
366.973 m2 mu + 10000. mu^2 -
m1 m2 mu^2 + (8.54564*10^6 - 46035.8 m1 - 50646.7 m2 + 250. m1 m2 -
1646.96 mu - 250. mu^2 + 1. m1 mu^2 + 1. m2 mu^2) t} *)


The solution is not pretty though.

Quiet[Solve[CoefficientList[rem, t] == 0, {m1, m2}]] // ExpandAll

(* Out[312]= {{m1 ->
1.13358*10^22/(1.02189*10^20 + 1.98614*10^17 mu -
2.23666*10^16 mu^2 +
6.20675*10^11 mu^4) + (1.99866*10^19 mu)/(1.02189*10^20 +
1.98614*10^17 mu - 2.23666*10^16 mu^2 +
6.20675*10^11 mu^4) - (2.65203*10^18 mu^2)/(1.02189*10^20 +
1.98614*10^17 mu - 2.23666*10^16 mu^2 +
6.20675*10^11 mu^4) + (2.27771*10^14 mu^3)/(1.02189*10^20 +
1.98614*10^17 mu - 2.23666*10^16 mu^2 +
6.20675*10^11 mu^4) + (7.75844*10^13 mu^4)/(1.02189*10^20 +
1.98614*10^17 mu - 2.23666*10^16 mu^2 +
6.20675*10^11 mu^4) - (3.44366*10^-23 \[Sqrt](5.41827*10^88 +
2.24318*10^86 mu - 2.13501*10^85 mu^2 -
4.49136*10^82 mu^3 + 2.77084*10^81 mu^4 +
1.31209*10^78 mu^5 - 1.25077*10^77 mu^6 -
1.22042*10^58 mu^7 + 1.82731*10^72 mu^8))/(1.02189*10^20 +
1.98614*10^17 mu - 2.23666*10^16 mu^2 + 6.20675*10^11 mu^4),
m2 -> 4.50593*10^44/(2.1364*10^42 + 1.65355*10^39 mu -
7.32605*10^38 mu^2 +
1.80237*10^34 mu^4) + (2.96843*10^41 mu)/(2.1364*10^42 +
1.65355*10^39 mu - 7.32605*10^38 mu^2 +
1.80237*10^34 mu^4) - (1.54024*10^41 mu^2)/(2.1364*10^42 +
1.65355*10^39 mu - 7.32605*10^38 mu^2 +
1.80237*10^34 mu^4) + (2.96843*10^37 mu^3)/(2.1364*10^42 +
1.65355*10^39 mu - 7.32605*10^38 mu^2 +
1.80237*10^34 mu^4) + (4.50593*10^36 mu^4)/(2.1364*10^42 +
1.65355*10^39 mu - 7.32605*10^38 mu^2 + 1.80237*10^34 mu^4) -
3.36384*10^64/((1.02189*10^20 + 1.98614*10^17 mu -
2.23666*10^16 mu^2 + 6.20675*10^11 mu^4) (2.1364*10^42 +
1.65355*10^39 mu - 7.32605*10^38 mu^2 +
1.80237*10^34 mu^4)) - (1.24688*10^62 mu)/((1.02189*10^20 +
1.98614*10^17 mu - 2.23666*10^16 mu^2 +
6.20675*10^11 mu^4) (2.1364*10^42 + 1.65355*10^39 mu -
7.32605*10^38 mu^2 +
1.80237*10^34 mu^4)) + (1.51171*10^61 mu^2)/((1.02189*10^20 \
+ 1.98614*10^17 mu - 2.23666*10^16 mu^2 +
6.20675*10^11 mu^4) (2.1364*10^42 + 1.65355*10^39 mu -
7.32605*10^38 mu^2 +
1.80237*10^34 mu^4)) + (2.7601*10^58 mu^3)/((1.02189*10^20 +
1.98614*10^17 mu - 2.23666*10^16 mu^2 +
6.20675*10^11 mu^4) (2.1364*10^42 + 1.65355*10^39 mu -
7.32605*10^38 mu^2 +
1.80237*10^34 mu^4)) - (2.15835*10^57 mu^4)/((1.02189*10^20 \
+ 1.98614*10^17 mu - 2.23666*10^16 mu^2 +
6.20675*10^11 mu^4) (2.1364*10^42 + 1.65355*10^39 mu -
7.32605*10^38 mu^2 +
1.80237*10^34 mu^4)) - (6.59765*10^53 mu^5)/((1.02189*10^20 \
+ 1.98614*10^17 mu - 2.23666*10^16 mu^2 +
6.20675*10^11 mu^4) (2.1364*10^42 + 1.65355*10^39 mu -
7.32605*10^38 mu^2 +
1.80237*10^34 mu^4)) + (9.81906*10^52 mu^6)/((1.02189*10^20 \
+ 1.98614*10^17 mu - 2.23666*10^16 mu^2 +
6.20675*10^11 mu^4) (2.1364*10^42 + 1.65355*10^39 mu -
7.32605*10^38 mu^2 +
1.80237*10^34 mu^4)) - (4.10528*10^48 mu^7)/((1.02189*10^20 \
+ 1.98614*10^17 mu - 2.23666*10^16 mu^2 +
6.20675*10^11 mu^4) (2.1364*10^42 + 1.65355*10^39 mu -
7.32605*10^38 mu^2 +
1.80237*10^34 mu^4)) - (1.39836*10^48 mu^8)/((1.02189*10^20 \
+ 1.98614*10^17 mu - 2.23666*10^16 mu^2 +
6.20675*10^11 mu^4) (2.1364*10^42 + 1.65355*10^39 mu -
7.32605*10^38 mu^2 +
1.80237*10^34 mu^4)) + (1.02189*10^20 \[Sqrt](5.41827*10^88 \
+ 2.24318*10^86 mu - 2.13501*10^85 mu^2 - 4.49136*10^82 mu^3 +
2.77084*10^81 mu^4 + 1.31209*10^78 mu^5 -
1.25077*10^77 mu^6 - 1.22042*10^58 mu^7 +
1.82731*10^72 mu^8))/((1.02189*10^20 + 1.98614*10^17 mu -
2.23666*10^16 mu^2 + 6.20675*10^11 mu^4) (2.1364*10^42 +
1.65355*10^39 mu - 7.32605*10^38 mu^2 +
1.80237*10^34 mu^4)) + (1.98614*10^17 mu \
\[Sqrt](5.41827*10^88 + 2.24318*10^86 mu - 2.13501*10^85 mu^2 -
4.49136*10^82 mu^3 + 2.77084*10^81 mu^4 +
1.31209*10^78 mu^5 - 1.25077*10^77 mu^6 -
1.22042*10^58 mu^7 + 1.82731*10^72 mu^8))/((1.02189*10^20 +
1.98614*10^17 mu - 2.23666*10^16 mu^2 +
6.20675*10^11 mu^4) (2.1364*10^42 + 1.65355*10^39 mu -
7.32605*10^38 mu^2 +
1.80237*10^34 mu^4)) - (2.23666*10^16 mu^2 \
\[Sqrt](5.41827*10^88 + 2.24318*10^86 mu - 2.13501*10^85 mu^2 -
4.49136*10^82 mu^3 + 2.77084*10^81 mu^4 +
1.31209*10^78 mu^5 - 1.25077*10^77 mu^6 -
1.22042*10^58 mu^7 + 1.82731*10^72 mu^8))/((1.02189*10^20 +
1.98614*10^17 mu - 2.23666*10^16 mu^2 +
6.20675*10^11 mu^4) (2.1364*10^42 + 1.65355*10^39 mu -
7.32605*10^38 mu^2 +
1.80237*10^34 mu^4)) + (6.20675*10^11 mu^4 \
\[Sqrt](5.41827*10^88 + 2.24318*10^86 mu - 2.13501*10^85 mu^2 -
4.49136*10^82 mu^3 + 2.77084*10^81 mu^4 +
1.31209*10^78 mu^5 - 1.25077*10^77 mu^6 -
1.22042*10^58 mu^7 + 1.82731*10^72 mu^8))/((1.02189*10^20 +
1.98614*10^17 mu - 2.23666*10^16 mu^2 +
6.20675*10^11 mu^4) (2.1364*10^42 + 1.65355*10^39 mu -
7.32605*10^38 mu^2 + 1.80237*10^34 mu^4))}, {m1 ->
1.13358*10^22/(1.02189*10^20 + 1.98614*10^17 mu -
2.23666*10^16 mu^2 +
6.20675*10^11 mu^4) + (1.99866*10^19 mu)/(1.02189*10^20 +
1.98614*10^17 mu - 2.23666*10^16 mu^2 +
6.20675*10^11 mu^4) - (2.65203*10^18 mu^2)/(1.02189*10^20 +
1.98614*10^17 mu - 2.23666*10^16 mu^2 +
6.20675*10^11 mu^4) + (2.27771*10^14 mu^3)/(1.02189*10^20 +
1.98614*10^17 mu - 2.23666*10^16 mu^2 +
6.20675*10^11 mu^4) + (7.75844*10^13 mu^4)/(1.02189*10^20 +
1.98614*10^17 mu - 2.23666*10^16 mu^2 +
6.20675*10^11 mu^4) + (3.44366*10^-23 \[Sqrt](5.41827*10^88 +
2.24318*10^86 mu - 2.13501*10^85 mu^2 -
4.49136*10^82 mu^3 + 2.77084*10^81 mu^4 +
1.31209*10^78 mu^5 - 1.25077*10^77 mu^6 -
1.22042*10^58 mu^7 + 1.82731*10^72 mu^8))/(1.02189*10^20 +
1.98614*10^17 mu - 2.23666*10^16 mu^2 + 6.20675*10^11 mu^4),
m2 -> 4.50593*10^44/(2.1364*10^42 + 1.65355*10^39 mu -
7.32605*10^38 mu^2 +
1.80237*10^34 mu^4) + (2.96843*10^41 mu)/(2.1364*10^42 +
1.65355*10^39 mu - 7.32605*10^38 mu^2 +
1.80237*10^34 mu^4) - (1.54024*10^41 mu^2)/(2.1364*10^42 +
1.65355*10^39 mu - 7.32605*10^38 mu^2 +
1.80237*10^34 mu^4) + (2.96843*10^37 mu^3)/(2.1364*10^42 +
1.65355*10^39 mu - 7.32605*10^38 mu^2 +
1.80237*10^34 mu^4) + (4.50593*10^36 mu^4)/(2.1364*10^42 +
1.65355*10^39 mu - 7.32605*10^38 mu^2 + 1.80237*10^34 mu^4) -
3.36384*10^64/((1.02189*10^20 + 1.98614*10^17 mu -
2.23666*10^16 mu^2 + 6.20675*10^11 mu^4) (2.1364*10^42 +
1.65355*10^39 mu - 7.32605*10^38 mu^2 +
1.80237*10^34 mu^4)) - (1.24688*10^62 mu)/((1.02189*10^20 +
1.98614*10^17 mu - 2.23666*10^16 mu^2 +
6.20675*10^11 mu^4) (2.1364*10^42 + 1.65355*10^39 mu -
7.32605*10^38 mu^2 +
1.80237*10^34 mu^4)) + (1.51171*10^61 mu^2)/((1.02189*10^20 \
+ 1.98614*10^17 mu - 2.23666*10^16 mu^2 +
6.20675*10^11 mu^4) (2.1364*10^42 + 1.65355*10^39 mu -
7.32605*10^38 mu^2 +
1.80237*10^34 mu^4)) + (2.7601*10^58 mu^3)/((1.02189*10^20 +
1.98614*10^17 mu - 2.23666*10^16 mu^2 +
6.20675*10^11 mu^4) (2.1364*10^42 + 1.65355*10^39 mu -
7.32605*10^38 mu^2 +
1.80237*10^34 mu^4)) - (2.15835*10^57 mu^4)/((1.02189*10^20 \
+ 1.98614*10^17 mu - 2.23666*10^16 mu^2 +
6.20675*10^11 mu^4) (2.1364*10^42 + 1.65355*10^39 mu -
7.32605*10^38 mu^2 +
1.80237*10^34 mu^4)) - (6.59765*10^53 mu^5)/((1.02189*10^20 \
+ 1.98614*10^17 mu - 2.23666*10^16 mu^2 +
6.20675*10^11 mu^4) (2.1364*10^42 + 1.65355*10^39 mu -
7.32605*10^38 mu^2 +
1.80237*10^34 mu^4)) + (9.81906*10^52 mu^6)/((1.02189*10^20 \
+ 1.98614*10^17 mu - 2.23666*10^16 mu^2 +
6.20675*10^11 mu^4) (2.1364*10^42 + 1.65355*10^39 mu -
7.32605*10^38 mu^2 +
1.80237*10^34 mu^4)) - (4.10528*10^48 mu^7)/((1.02189*10^20 \
+ 1.98614*10^17 mu - 2.23666*10^16 mu^2 +
6.20675*10^11 mu^4) (2.1364*10^42 + 1.65355*10^39 mu -
7.32605*10^38 mu^2 +
1.80237*10^34 mu^4)) - (1.39836*10^48 mu^8)/((1.02189*10^20 \
+ 1.98614*10^17 mu - 2.23666*10^16 mu^2 +
6.20675*10^11 mu^4) (2.1364*10^42 + 1.65355*10^39 mu -
7.32605*10^38 mu^2 +
1.80237*10^34 mu^4)) - (1.02189*10^20 \[Sqrt](5.41827*10^88 \
+ 2.24318*10^86 mu - 2.13501*10^85 mu^2 - 4.49136*10^82 mu^3 +
2.77084*10^81 mu^4 + 1.31209*10^78 mu^5 -
1.25077*10^77 mu^6 - 1.22042*10^58 mu^7 +
1.82731*10^72 mu^8))/((1.02189*10^20 + 1.98614*10^17 mu -
2.23666*10^16 mu^2 + 6.20675*10^11 mu^4) (2.1364*10^42 +
1.65355*10^39 mu - 7.32605*10^38 mu^2 +
1.80237*10^34 mu^4)) - (1.98614*10^17 mu \
\[Sqrt](5.41827*10^88 + 2.24318*10^86 mu - 2.13501*10^85 mu^2 -
4.49136*10^82 mu^3 + 2.77084*10^81 mu^4 +
1.31209*10^78 mu^5 - 1.25077*10^77 mu^6 -
1.22042*10^58 mu^7 + 1.82731*10^72 mu^8))/((1.02189*10^20 +
1.98614*10^17 mu - 2.23666*10^16 mu^2 +
6.20675*10^11 mu^4) (2.1364*10^42 + 1.65355*10^39 mu -
7.32605*10^38 mu^2 +
1.80237*10^34 mu^4)) + (2.23666*10^16 mu^2 \
\[Sqrt](5.41827*10^88 + 2.24318*10^86 mu - 2.13501*10^85 mu^2 -
4.49136*10^82 mu^3 + 2.77084*10^81 mu^4 +
1.31209*10^78 mu^5 - 1.25077*10^77 mu^6 -
1.22042*10^58 mu^7 + 1.82731*10^72 mu^8))/((1.02189*10^20 +
1.98614*10^17 mu - 2.23666*10^16 mu^2 +
6.20675*10^11 mu^4) (2.1364*10^42 + 1.65355*10^39 mu -
7.32605*10^38 mu^2 +
1.80237*10^34 mu^4)) - (6.20675*10^11 mu^4 \
\[Sqrt](5.41827*10^88 + 2.24318*10^86 mu - 2.13501*10^85 mu^2 -
4.49136*10^82 mu^3 + 2.77084*10^81 mu^4 +
1.31209*10^78 mu^5 - 1.25077*10^77 mu^6 -
1.22042*10^58 mu^7 + 1.82731*10^72 mu^8))/((1.02189*10^20 +
1.98614*10^17 mu - 2.23666*10^16 mu^2 +
6.20675*10^11 mu^4) (2.1364*10^42 + 1.65355*10^39 mu -
7.32605*10^38 mu^2 + 1.80237*10^34 mu^4))}} *)


The scale differences in coefficients lead me to suspect the root expresssions might be not well behaved numerically. If that is in fact the case, it might be more useful to instead use a function that computes m1 and m2 numerically when given numerical values for mu.

• That's a very good idea, it worked perfectly! I checked using some expected values and the results are correct.
– Anna
Aug 13, 2023 at 7:44
• I just noticed that there is a problem: in this way you don't find the solutions for the specific eigenvalues 3 and 2 but for two random eigenvalues. Thing is, for the purpose of my research, the ordering of the four eigenvalues is established and relevant. Setting a low value of mu and evaluating the four eigenvalues using the obtained formulas for M1 and M2, gives that e2 is indeed 50, while it is e4 the one with the value 200; yet, this it's not always true, in fact setting a high mu the two eigenvalues involved are the 3rd and the 2nd as they should but it's by a welcome chance.
– Anna
Aug 16, 2023 at 10:26
• Matrix eigenvalues are a set. The result from the Eigenvalues function imposes an ordering on them, but it is not canonical in any mathematical sense. So more work would be needed to impose a particular ordering relevant to the work at hand. Aug 16, 2023 at 15:39
Clear["Global*"]


Use exact numbers

m = {{M1, 0, -4.28345, 42.8361}, {0, M2, 7.9998, -80.001}, {-4.28345, 7.9998,
0, -mu}, {42.8361, -80.001, -mu, 0}} // Rationalize;

eig = Eigenvalues[m];

eqns23 = {eig[[2]] == 50, eig[[3]] == 200};

sol23 = Solve[eqns23, {M1, M2}] // FullSimplify;

(* Solve::nongen: There may be values of the parameters for which some or all solutions are not valid. *)

eig /. sol23 /. ({mu -> #} & /@ Range[3]) // FullSimplify


The target values equate to eig[[3]] and eig[[4]]

eqns34 = {eig[[3]] == 50, eig[[4]] == 200};

sol34 = Solve[eqns34, {M1, M2}] // FullSimplify;

(* Solve::nongen: There may be values of the parameters for which some or all solutions are not valid. *)

((And @@ eqns34) /. sol34 /. ({mu -> #} & /@ Range[10]) // FullSimplify)

(* {{True, True}, {True, True}, {True, True}, {True, True}, {True, True}, {True,
True}, {True, True}, {True, True}, {True, True}, {True, True}} *)

SeedRandom[1234];

((And @@ eqns34) /.
sol34 /. ({mu -> #} & /@ Rationalize[RandomReal[1, 10], 0]) //
FullSimplify)

(* {{True, True}, {True, True}, {True, True}, {True, True}, {True, True}, {True,
True}, {True, True}, {True, True}, {True, True}, {True, True}} *)
`
• Thanks for your answer. This procedure seems to highlight my problem: how is it possible that solving eqns32, actually gives the solutions for 3,4...
– Anna
Aug 16, 2023 at 10:51