Trouble using Solve and NSolve with functions involving Erf

Question

I have the following functions:

R[k_, x_, t_] := -.5*(k - x)*(1 + Erf[-(k - x)/t])
L[k_, c_, x_, t_] := .5*c*(k - x)*(1 + Erf[(k - x)/t])

I'm interested in finding the k that solves some equation using these two functions. For example:

Sum[(R[k, i, .000001]^1 + L[k, 1, i, .000001]^2)*
  Binomial[40, i]*.25^i*.75^(40 - i), {i, 0, 40}] == 10

The lefthand side plots well, but I am having no luck with

Solve[
  Sum[(R[k, i, .00001]^1 + L[k, 1, i, .00001]^2)*
    Binomial[40, i]*.25^i*.75^(40 - i), {i, 0, 40}] == 10, k]

or with using NSolve or N[Erf[...]] in the function definitions. The error message is

Solve::inex: Solve was unable to solve the system with inexact coefficients or the system obtained by direct rationalization of inexact numbers present in the system. Since many of the methods used by Solve require exact input, providing Solve with an exact version of the system may help.

Any ideas as to why this won't work or how I can get this to work?

Answer 1

The equations have some pretty large and small numbers at work in them. Sometimes that can be a problem because of the limits of a finite machine. Some of the limits are

$MaxNumber                                                                                  (*$*)
$MinNumber                                                                                  (*$*)
$MachineEpsilon
(*
  1.605216761933662*10^1355718576299609
  6.229688249675322*10^-1355718576299610
  2.22045*10^-16
*)

($MachineEpsilon is the limit on the ratio of terms x > y such that x + y is different than x. Even when the ratio close, only a few higher-order bits of the smaller number are significant in the sum. For sensitive system with lots of such sums, this might be an issue, although it does not seem to be so in this case.)

So first, let's convert the formulas to exact expressions. I also substitute Erfc for 1 - Erf, to avoid round-off error, since the arguments to Erf range in the positive and negative millions and the values of Erf are identical to ± 1 at machine precision.

R[k_, x_, t_] := -(1/2)*(k - x)*Erfc[(k - x)/t];
L[k_, c_, x_, t_] := (1/2)*c*(k - x)*Erfc[-(k - x)/t];

Clear[f]
f[k_] = Sum[(R[k, i, [email protected]]^1 + L[k, 1, i, [email protected]]^2)*
    Binomial[40, i]*(1/4)^i*(3/4)^(40 - i), {i, 0, 40}];

We can get an idea of what's going on by plotting (ignoring the underflow messages that come from large arguments to Erfc):

Plot[f[k], {k, -10, 15}]

There appear to be just two roots, one at near 0 and one around 11 or 12. We can search for them with NDSolve:

Last@ Reap@
  NDSolve[{k'[t] == 1, k[-10] == -10, 
    WhenEvent[f[t] == 10, 
     Sow[x /. FindRoot[f[x] == 10, {x, t}]]]}, {}, {t, -10, 20}]

General::unfl: Underflow occurred in computation. >>
...

(*
  {{-8.91463*10^-16, 11.8286}}
*)

I couldn't raise the working precision:

Last@Reap@
  NDSolve[{k'[t] == 1, k[-10] == -10, 
    WhenEvent[f[t] == 10, 
     Sow[x /. FindRoot[f[x] == 10, {x, t}, WorkingPrecision -> 20]]]}, {},
   {t, -10, 20}]

General::unfl: Underflow occurred in computation. >> ...
General::nomem: The current computation was aborted because there was insufficient memory available to complete the computation.
Throw::sysexc: Uncaught SystemException returned to top level. Can be caught with Catch[..., _SystemException]. >>

SystemException["MemoryAllocationFailure"]

Evaluating an Erfc with a very large argument at arbitrary precision seems to takes gigabytes.

Now let's look at the underflow warnings. If we look at the last term in f[0], we get

Last@f[0]
(*  (20 Erfc[-40000000] + 400 Erfc[40000000]^2)/1208925819614629174706176  *)

Now, Erfc[40000000] is a little less than Sqrt[$MinNumber]:

Erfc[40000000.]
(*  1.686668849131491*10^-694871171045211  *)

So squaring it will result in underflow. Since 400 times the square is being added to 20 Erfc[-40000000], it is effectively zero anyway.

If we look at all the terms of f at the zeros, we get

f[0] /. Plus -> Inactive[Plus] /. x_Integer | x_Rational :> N[x] // 
 Short[List @@ #, 20] &

General::unfl: Underflow occurred in computation. >>

(*
  {0.000134088 (1. + 2.495050185825339*10^-868588963820),
  0.000871571 (2. + 7.69*10^-3474355855240),
  0.00367997 (3. + 2.333625160219760*10^-7817300674272),
  <<34>>,
  5.80681*10^-21 (38. + 4.412141103756946*10^-1254242463736605),
  9.92617*10^-23 (39. + 5.915985858762761*10^-1321123813949706),
  8.27181*10^-25 (40. + Underflow[])}
*)

f[Rationalize[11.828562582665537`, 0]] /. Plus -> Inactive[Plus] /. 
  x_Integer | x_Rational :> N[x] // Short[List @@ #, 20] &
(*
  {0.105721 (0.171437 + 2.523363170467293*10^-25528514162),
   0.13124 (0.686516 + -1.420977915095513*10^-298150090304),
   0.0759025 (1.17144 + 3.100452009927070*10^-1191934775329),
   << 36 >>,
   9.92617*10^-23 (27.1714 + 1.248075143733985*10^-641267950157834),
   8.27181*10^-25 (28.1714 + 1.153984027103677*10^-689338160464552)}
*)

The upshot is that the two roots we found seem reliable. In f[0.] the error function is pegged at the extremes of its range, but in f[11.8] only half are effectively zero. Although f[0.] evaluates to 10. (exactly), it's difficult to see whether the root is exactly k == 0. Considering this,

f[0] /. Erfc[x_?Negative] :> 2 - Erfc[-x] /. _Erfc -> 0
(* 10 *)

all the error function terms would have to cancel out for f[0] to be exactly 10. But looking at the linear and squared Erfc terms, they do not appear to cancel, although they are too small to be considered practically different from 0.

foo - 10 /. {Erfc[_]^2 -> 0} // Simplify // N                  (* linear Erfc terms *)
foo - 10 /. {x : Erfc[_]^2 :> x, _Erfc -> 0} // Simplify // N  (* squared Erfc terms *)
(*
  -2.118062622431921*10^-434294481914
  3.345710199166185*10^-868588963824
*)

(Note that the first rule in {x : Erfc[_]^2 :> x, _Erfc -> 0} prevents the Erfc in a squared term from being replaced by 0.)

Answer 2

Solve's core functionality is symbolic algebraic (polynomial) equations. Your equation is far from algebraic, and Solve can't find a transformation that would make it algebraic. Since you seem to want a numerical answer, a numerical method is probably more appropriate:

FindRoot[Sum[(R[k, i, .00001]^1 + L[k, 1, i, .00001]^2)* Binomial[40, i]*.25^i*.75^(40 - i), {i, 0, 40}] == 10, {k, 1}] (* {k -> 6.36196*10^-16} *)

Guessing {k,0} would be cheating ;-)