# mathematica is not able to solve this system

I have this system of three equation that mathematica does not resolve.

x = 1005790/1000000;
y = 3012034/1000000;
varx = 9979050/1000000;
vary = 12985206/1000000;
covxy = 9986968/1000000;
z = 3614929/1000000;
varz = 25466937/1000000;
z3 = 323087800/1000000;
eq1 = z == av + at*x + a1*y
eq2 = varz == at^2 *varx + a1^2 *vary + 2*covxy *a1*at
eq3 = (av + x*at + a1*y)^3 +3*(av + x*at + a1*y)*(at^2 *varx + a1^2 * vary + 2* at*a1* covxy) == z3
FindRoot[{eq1 , eq2, eq3}, {{av, 0.4}, {at, 0.7}, {a1, 1.1}}]
Reduce[{eq1 , eq2, eq3}, {av, at, a1}]


Reduce return False and FindRoot is printing the error: FindRoot::jsing: Encountered a singular Jacobian at the point {av,at,a1} = {0.4,0.7,1.1}. Try perturbing the initial point(s). But even if I change the initial points by a Little, The program is unable to solve the equations. I have plotted the three equations minus the known term, and graphically it seems that there is a solution (I fixed av to the correct value, and I see that the three functions are near zero at the same point).

code to plot:

prim[av_, at_, a1_] := av + at*x + a1*y - z;
sec [av_, at_, a1_] := at^2 *varx + a1^2 *vary + 2*covxy *a1*at - varz;
terz [av_, at_, a1_] := (av + x*at + a1*y)^3 +
3*(av + x*at + a1*y)*(at^2 *varx + a1^2 * vary +
2* at*a1* covxy) - z3;
Plot3D[{prim[0.1, at, a1], sec[0.1, at, a1], terz[0.1, at, a1]}, {at,
0.4, 0.6}, {a1, 0.9, 1.1}]

• Can you add the code you used to plot the functions and the value of av you tried? The fact that Reduce returns False seems a strong indication that there is no solution. Aug 2, 2022 at 12:47
• Sure! I changed the question Aug 2, 2022 at 13:12
• There is no actual solution (the system is inconsistent). Aug 2, 2022 at 19:17

Try NMinimize:

eq = {-z + av + at*x + a1*y,-varz + at^2*varx + a1^2*vary + 2*covxy*a1*at, (av + x*at + a1*y)^3 +3*(av+ x*at + a1*y)*(at^2*varx + a1^2*vary +2*at*a1*covxy) - z3};

sol=NMinimize[eq . eq, {a1, at, av}]
(*{8.30607*10^-6, {a1 -> 0.966855, at ->0.539657, av -> 0.157078}}*)
eq /. sol[]
(* {-0.00286933, -0.000269104, 0.0000248339}*)


If we look examplary for positive solutions, two ContourPlot's show

sol1 = Solve[eq[] == 0, a1][[-1]]
Column[{
ContourPlot[eq[] == 0, {at, 0, 2}, {a1, 0, 2},
FrameLabel -> {at, a1}],
ContourPlot[
Evaluate[Thread[eq[[{1, 3}]] == 0 /. sol1]], {at, 0, 2}, {av, 0,
2}, ContourStyle -> {Blue, {Dashed, Red}},
FrameLabel -> {at, av}]}] that solutions a1[at],av [at]exist for ~.5<at<1.5

• Thank you! You have found exactly what I needed. The coefficients of the system are statistics, so the correct way is to minimize with the three equation. Now It's clear. Thanks Aug 2, 2022 at 13:32
• @AdrianoDelVincio You're welcome. See my addendum! Aug 2, 2022 at 13:42