# Trouble using Solve and NSolve with functions involving Erf

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?

• Welcome to Mathematica.SE! I suggest that: 1) You take the introductory Tour now! 2) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! 3) As you receive help, try to give it too, by answering questions in your area of expertise. – bbgodfrey Mar 4 '15 at 16:40
• The message seems clear. You are using inexact numbers when you shouldn't. Try specifying .00001 as 1/100000, .25 as 1/4, etc. – m_goldberg Mar 4 '15 at 16:44
• @m_goldberg Could you please link a question (if you have one in mind) which discusses how Mathematica views 0.25 as different from 1/4? I'm halfway through the Shifrin text and that nuance hasn't come up yet (I don't think). – Shane Mar 4 '15 at 16:54
• You have got to find them all and get rid of them. Even the ones in the definitions of R and L. I did that and I do not get your message. I get get "Solve::nsmet: This system cannot be solved with the methods available to Solve. >>" instead, but that is a whole different matter. – m_goldberg Mar 4 '15 at 17:01
• I suppose mathematica.stackexchange.com/questions/24444/… is relevant. But in what sense are 1/4 and 0.25 distinct? 1/4 == 0.25 evaluates to true, after all. – Shane Mar 4 '15 at 17:02

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, Rationalize@.000001]^1 + L[k, 1, i, Rationalize@.000001]^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.)

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 ;-)