# Approximate solution to system of polynomial equations

I'd like to solve approximately the following system of equations for $a_1$, $a_2$, $\alpha_1$ and $\alpha_2$ all real and nonzero:

$a_1(1+4\alpha_2)^3+8a_2(1+4\alpha_1)^2+\frac{\sqrt{\pi}}{64}(1+4\alpha_1)^2(1+4\alpha_2)^3=0$

$a_1(1+4\alpha_2)^5+16a_2(1+4\alpha_1)^4-\frac{11\sqrt{\pi}}{512}(1+4\alpha_1)^4(1+4\alpha_2)^5=0$

$a_1(1+4\alpha_2)^7+24a_2(1+4\alpha_1)^6-\frac{65\sqrt{\pi}}{8192}(1+4\alpha_1)^6(1+4\alpha_2)^7=0$

$a_1(1+4\alpha_2)^9+32a_2(1+4\alpha_1)^8-\frac{291\sqrt{\pi}}{65536}(1+4\alpha_1)^8(1+4\alpha_2)^9=0$

where I know that some approximate solutions are

$a_1 = 117.429,\ \ a_2 = -755.468,\ \ \alpha_1 = 1.3177,\ \ \alpha_2 = 2.9026$

I'd like to improve on these or at least confirm them. Plotting the equations in one parameter while keeping the other three as above I see that there are indeed roots there, e.g:

a2 = -755.468; alpha1 = 1.3177; alpha2 = 2.9026;


and:

Plot[a1/(1 + 4 alpha1)^2 + (8 a2)/(1 + 4 alpha2)^3 +  Sqrt[Pi]/64, {a1, 110, 120},
PlotRange -> {{110, 120}, {-0.1, 0.1}}]


I've tried using FindRoot with different parameters (accuracies, iterations, methods), but I end up with Failed to converge to the requested accuracy or precision within # iterations, and results which are completely different. One example:

FindRoot[{a1 (1 + 4alpha2)^3 + 8 a2 (1 + 4alpha1)^2 +
Sqrt[Pi] (1 +4 alpha1)^2 (1 +4 alpha2)^3/64 ,
a1 (1 +4 alpha2)^5 + 16 a2 (1 +4 alpha1)^4 -
11 Sqrt[Pi] (1 + 4alpha1)^4 (1 + 4alpha2)^5/512 ,
a1 (1 +4 alpha2)^7 + 24 a2 (1 + 4alpha1)^6 -
65 Sqrt[Pi] (1 +4 alpha1)^6 (1 + 4alpha2)^7/8192 ,
a1 (1 +4 alpha2)^9 + 32 a2 (1 + 4alpha1)^8 -
291 Sqrt[Pi] (1 + 4alpha1)^8 (1 + 4alpha2)^9/65536}, {{a1, 110,
120}, {a2, -760, -750}, {alpha1, 1.3, 1.4}, {alpha2, 2.9, 3}},
MaxIterations -> 100000, WorkingPrecision -> 50]


Gives the error above and

$a_1 \rightarrow 0.631165832,$

$a_2 \rightarrow -70.5406787,$

$\alpha_1 \rightarrow 0.30522891,$

$\alpha_2 \rightarrow 3.7582495135$.

Any ideas how I can improve on this?

You can actually solve the whole system:

sol = Solve[{a1 (1 + alpha2)^3 + 8 a2 (1 + alpha1)^2 +
Sqrt[Pi] (1 + alpha1)^2 (1 + alpha2)^3/64 == 0,
a1 (1 + alpha2)^5 + 16 a2 (1 + alpha1)^4 -
11 Sqrt[Pi] (1 + alpha1)^4 (1 + alpha2)^5/512 == 0,
a1 (1 + alpha2)^7 + 24 a2 (1 + alpha1)^6 -
65 Sqrt[Pi] (1 + alpha1)^6 (1 + alpha2)^7/8192 == 0,
a1 (1 + alpha2)^9 + 32 a2 (1 + alpha1)^8 -
291 Sqrt[Pi] (1 + alpha1)^8 (1 + alpha2)^9/65536 == 0}, {a1, a2,
alpha1, alpha2}];


This returns a very large equation with 39 different solutions. You can look at the individual solutions:

sol[]//N


or scan through them:

Manipulate[sol[[i]] // N, {i, 1, Length[sol], 1}] • Thanks for the answer! It highlighted what I forgot to mention in the OP: the coefficients must be real and nonzero. With that constraint Solve doesn't give any answers. Sorry about that Sep 25 '13 at 14:30

For confirmation of your approximate root you can simply replace it into the expression and see if everything is close to 0:

expr = {a1 (1 + 4 alpha2)^3 + 8 a2 (1 + 4 alpha1)^2 +
Sqrt[Pi] (1 + 4 alpha1)^2 (1 + 4 alpha2)^3/64,
a1 (1 + 4 alpha2)^5 + 16 a2 (1 + 4 alpha1)^4 -
11 Sqrt[Pi] (1 + 4 alpha1)^4 (1 + 4 alpha2)^5/512,
a1 (1 + 4 alpha2)^7 + 24 a2 (1 + 4 alpha1)^6 -
65 Sqrt[Pi] (1 + 4 alpha1)^6 (1 + 4 alpha2)^7/8192,
a1 (1 + 4 alpha2)^9 + 32 a2 (1 + 4 alpha1)^8 -
291 Sqrt[Pi] (1 + 4 alpha1)^8 (1 + 4 alpha2)^9/65536};

knownApprox = {a1 -> 117.429, a2 -> -755.468, alpha1 -> 1.3177, alpha2 -> 2.9026};
expr/.knownApprox//Abs//Max
(* 1.50862*10^14    oops *)


At least NMinimize can find a solution by minimizing the sum of the absolute values:

NMinimize[Total@Abs@expr, {a1, a2, alpha1, alpha2}]
(* {3.45995*10^-6,
{a1 -> -0.50239, a2 -> 1.15054, alpha1 -> -0.251652, alpha2 -> -0.226791}} *)

(* You can put in constraints if you want a solution within some region *)
NMinimize[{Total@Abs@expr, a1 > 110 && a2 < -700}, {a1, a2, alpha1, alpha2}]


There is lots of good documentation on NMinimize in the Optimization tutorails

• Thanks! This shows a mistake I made: the system of equations in the OP have been manipulated as if the $=0$'s are exact; namely I multiplied the first equation with $(1+4\alpha_1)^2(1+4\alpha_2)^3$ and correspondingly for the rest (the original equations are like the one in my Plot in the OP). Using the original equations with the stated coefficients the sum of the absolute values is $0.0204887$, slightly better than $0.060103232$, which NMinimize gives us with WorkingPrecision $50$ and $10^6$ iterations. Still not sure how the stated coefficients came about. Sep 25 '13 at 15:30
• Actually, using FindRoot with the original equations gives answers pretty close to the stated coefficients, presumably I can get closer with higher number of iterations. Sep 25 '13 at 15:43