More issues Integrate, NIntegrate, FindRoot

Question

I'm trying to solve a system of 3 non-linear equations using FindRoot and Integrate. If I start FindRoot[] close to the right answer, it works well but returns a bunch of error messages first.

f[x_?NumericQ, a_?NumericQ, b_?NumericQ, c_?NumericQ] := Exp[-(a*x^4 + b*x^2 + c)];
aa = FindRoot[Integrate[(x^{0, 2, 4})*f[x, a, b, c], {x, -Infinity, Infinity}] - {1, 2, 10} // N, {a, .01}, {b, .15}, {c, 1.33}]

(*  NIntegrate::inumr: "The integrand f[x,a,b,c] has evaluated to non-numerical values for all...   *)
(*  {a -> 0.0108063, b -> 0.141937, c -> 1.36499}  *)

But if I start a ways farther from the right answer, it give me pages of error messages and returns my starting values. I actually ended up having to write a random-walk algorithm to get close enough to the right answer for FindRoot[] to finish the job.

The command that failed completely was this:

aa = FindRoot[NIntegrate[(x^{0, 2, 4})*f[x, a, b, c], {x, -Infinity,Infinity}] - {1, 2, 10} // N, {a, .1}, {b, .5}, {c, .33}]

My crude random-walk algorithm handled this starting point just fine.

So even though I got an answer eventually, was there a better approach I could have taken, that wouldn't have required me to find a nearly-correct answer first? And why am I getting all those warning messages before I get the right answer?

Answer 1

The error messages come from Integrate/NIntegrate calls on functions with symbolic (non-numerical) parameters a, b, c. These calls come before FindRoot substitutes numerical values for these parameters. (Taking N of Integrate here is effectively the same as calling NIntegrate.)

In this case the best approach it seems to me is to do the integral first, since it can be done.

obj = Integrate[
  x^{0, 2, 4} Exp[-(a*x^4 + b*x^2 + c)], {x, 0, Infinity}, 
  Assumptions -> a > 0 && b > 0 && c > 0]
(*
{1/4 Sqrt[b/a] E^(b^2/(8 a) - c) BesselK[1/4, b^2/(8 a)],
 (1/(16 Sqrt[2] Sqrt[a^3 b])) *
  E^(b^2/(8 a) - c) π (-b^2 BesselI[-(1/4), b^2/(8 a)] +
   (4 a + b^2) BesselI[1/4, b^2/(8 a)] + 
    b^2 (-BesselI[3/4, b^2/(8 a)] + BesselI[5/4, b^2/(8 a)])),
 (Sqrt[b] E^(b^2/(8 a) - c) ((2 a + b^2) BesselK[1/4, b^2/(8 a)] - 
    b^2 BesselK[3/4, b^2/(8 a)]))/(32 a^(5/2))}
*)

aa = FindRoot[obj - {1, 2, 10}, {a, .01}, {b, .15}, {c, 1.33}]
(*  {a -> 0.0108063, b -> 0.141937, c -> 0.671843}  *)

Checks:

obj /. aa
(*  {1., 2., 10.}  *)

NIntegrate[x^{0, 2, 4} Exp[-(a*x^4 + b*x^2 + c)] /. aa, {x, 0, Infinity}]
(*  {1., 2., 10.}  *)

---

Further explanation:

Integrating vector expressions can be tricky. See NIntegrate over a list of functions and linked questions.

The integration is done on each component separately. NIntegrate must see the components in the argument expression. The following, which is the normal way to use NumericQ in numeric solvers like NIntegrate, does not work because NIntegrate decides the integrand is not a List and gets confused when the values are not numbers:

i1[x_?NumericQ, a_?NumericQ, b_?NumericQ, c_?NumericQ] :=
  (x^{0, 2, 4})*f[x, a, b, c]
FindRoot[NIntegrate[i1[x, a, b, c], {x, -Infinity, Infinity}] - {1, 2, 10},
 {a, .01}, {b, .15}, {c, 1.33}]

NIntegrate::inumr: The integrand i1[x,a,b,c] has evaluated to non-numerical values.... >>

The proper way is to wrap the function up for FindRoot instead:

obj2[a_?NumericQ, b_?NumericQ, c_?NumericQ] := 
 NIntegrate[(x^{0, 2, 4})*f[x, a, b, c], {x, -Infinity, Infinity}]

But there's still one pitfall. If we substitute this for the integration there's a problem:

FindRoot[obj2[a, b, c] - {1, 2, 10}, {a, .01}, {b, .15}, {c, 1.33}]

FindRoot::nveq: The number of equations does not match the number of variables in FindRoot[obj2[a,b,c]-{1,2,10},{a,0.01},{b,0.15},{c,1.33}]. >>

The problem is that the argument to FindRoot is evaluated symbolically before obj2 evaluates. You get a vector of obj2 calls:

obj2[a, b, c] - {1, 2, 10}
(*  {-1 + obj2[a, b, c], -2 + obj2[a, b, c], -10 + obj2[a, b, c]}  *)

When FindRoot substitutes values for a, b, and c, this evaluates and you get a matrix of values.

obj2[a, b, c] - {1, 2, 10} /. {a -> 0.01, b -> 0.15, c -> 1.33}
(*
  {{0.0273318, 1.04919, 9.31436},
   {-0.972668, 0.0491918, 8.31436},
   {-8.97267, -7.95081, 0.314357}}
*)

The proper way to use obj2 is to use == instead of -:

FindRoot[obj2[a, b, c] == {1, 2, 10}, {a, .01}, {b, .15}, {c, 1.33}]
(*  {a -> 0.0108063, b -> 0.141937, c -> 1.36499}  *)

Answer 2

The documentation of FindRoot says

FindRoot first localizes the values of all variables, then evaluates f with the variables being symbolic, and then repeatedly evaluates the result numerically.

So the system first evaluates

Integrate[(x^{0, 2, 4})*f[x, a, b, c], {x, -Infinity, Infinity}] - {1, 2, 10} // N

which gives the errors you see (NIntegrate does the same btw). First we make your problem compatible with this. This mostly means that arguments within these functions should not evaluate to give error messages for non-numerical arguments, see https://mathematica.stackexchange.com/a/26037/6804 and http://support.wolfram.com/kb/12502 .

This works:

ClearAll[g, f, a, b, c];

f[x_?NumericQ, a_?NumericQ, b_?NumericQ, c_?NumericQ] := 
  Exp[-(a*x^4 + b*x^2 + c)];
g[a_?NumericQ, b_?NumericQ, c_?NumericQ] := 
  NIntegrate[(x^{0, 2, 4})*f[x, a, b, c], {x, -Infinity, 
     Infinity}] - {1, 2, 10};

FindRoot[g[a, b, c], {a, .01}, {b, .15}, {c, 1.33}]

as you observed, FindRoot[g[a, b, c], {a, .1}, {b, .5}, {c, .33}] doesn't. The messages states it clearly: NIntegrate::inumri: The integrand f[x,-3.12487,7.94363,-0.799115] has evaluated to Overflow, Indeterminate, or Infinity for all sampling points in the region with boundaries {{0.,17545.1}}. >>, since f[17545.130534855336`, -3.124874458858545`, 7.943629253761611`, -0.7991150536849108`] === Overflow[].

Your exponential becomes gigantic for negative a and large x which can very well occur in the algorithm FindRoot uses which can involve evaluating the function anywhere it sees fit. Sometimes it will not be able to find the solution...

Instead of FindRoot, let's use NMinimize which supports constraints - I'm assuming you only want positive a and b, so I use a > 0 && b > 0. If there is a root, NMinimize should find it too. We have to change the problem a bit because NMinimize takes a single valued function: apply Norm. The NelderMead method seems to find the desired solution.

NMinimize[{Norm@g[a, b, c], a > 0 && b > 0}, {a, b, c},
 MaxIterations -> 100,
 StepMonitor :> Print[Norm@g[a, b, c], {a, b, c}],
 Method -> "NelderMead"]
(*=>*)
{0.000163492, {a -> 0.0107889, b -> 0.142143, c -> 1.36467}}

Hoever note that this kind of global rootfinding/optimization is in general very hard.