Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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?

share|improve this question
add comment

2 Answers

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[[20]]//N

or scan through them:

Manipulate[sol[[i]] // N, {i, 1, Length[sol], 1}]

enter image description here

share|improve this answer
    
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 –  jorgen Sep 25 '13 at 14:30
add comment

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

share|improve this answer
    
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. –  jorgen 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. –  jorgen Sep 25 '13 at 15:43
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.