Message Boards Message Boards

How to solve this difficult nonlinear equation involving InverseErf?

Posted 7 months ago

Hello,
I have a difficulty with solving the following nonlinear equation involving InverseErf[] using FindRoot:

FindRoot[Sqrt[p*r]/InverseErf[(2/Sqrt[π])*Sqrt[p]]==1,{r,1},WorkingPrecision->50,PrecisionGoal->16]

The problem is that when the argument of InverseErf[] approaches unity, the value of the function tends to infinity, but I need to obtain accurate values of "r" just close to this limit of the argument of InverseErf[]. Unfortunately, it seems it is not possible to have the parameter "p" greater than about p=78/100, which is still too small for my needs. The function InverseErf[] is supposed to return correct values for any argument between -1 and 1, but apparently it fails already when p=79/100. Is there any way to overcome this difficulty?
Leslaw.

POSTED BY: Leslaw Bieniasz
4 Replies
Posted 7 months ago

Sorry for the nonsensical formulation of my post; indeed the equation can be solved analytically when the InverseErf[] function is used, and I also incorrectly evaluated a limit on p. I was mislead by my initial attempts to use Erf[] instead of InverseErf[], and I would be happy to understand why such an approach causes warnings in FindRoot. My code is:

F[p_]:=If[p==0,1,FindRoot[Erf[Sqrt[p*r]]==((2/Sqrt[π])*Sqrt[p]),{r,2},WorkingPrecision->100,PrecisionGoal->16][[1]][[2]]];
Plot[F[p],{p,0,π/4},PlotRange->{0,30},PlotPoints->5000]

This produces a warning:

FindRoot::precw: The precision of the argument function (Erf[0.000396571 Sqrt[r]]==0.000447482) is less than WorkingPrecision (100.`).

Is there any way to get rid of this warning? I cannot figure out why it occurs and how serious it is supposed to be. The online help suggests it can be ignored, but is this true? I need to be sure that the results have a desired precision. The equation is a special case of a more complicated one, for which the approach with InverseErf[] cannot be used, so that the above code appears to be the only method.

Leslaw

POSTED BY: Leslaw Bieniasz
Posted 7 months ago

This

F[p_]:=(Print[{p,Precision[p]}];If[p==0,1,FindRoot[Erf[Sqrt[p*r]]==((2/Sqrt[\[Pi]])*Sqrt[p]),{r,2},WorkingPrecision->100,PrecisionGoal->16][[1]][[2]]]);
Plot[F[p],{p,0,\[Pi]/4},PlotRange->{0,30},PlotPoints->5000]

shows that Plot is using MachinePrecision for each step of the calculation and passing those values to F[p], after recognizing that 0 and Pi are exact. That does not seem surprising. You are asking FindRoot to work with 100 digits of precision in the calculations it does and it can't produce that precision given the less precise MachinePrecision value for p. I believe that is the reason for the warning.

So I will try to sidestep that warning with

F[p_]:=(Print[{p,Precision[p]}];If[p==0,1,FindRoot[Erf[Sqrt[p*r]]==((2/Sqrt[\[Pi]])*Sqrt[p]),{r,2},WorkingPrecision->100,PrecisionGoal->16][[1]][[2]]]);
Plot[F[p],{p,0,\[Pi]/4},PlotRange->{0,30},PlotPoints->5000,WorkingPrecision->100]

which instructs Plot to use 100 digit math. That shows that Plot has made all the subsequent calculations done with 100 digits and thus pass 100 digit precision arguments to F[p], but that first step remains MachinePrecision which still triggers that warning. I don't have an explanation for that.

The warning remains even if I change your PrecisionGoal to 32 or 100 or any other value.

I suppose you could do this

F[p_]:=(q=SetPrecision[p,100];
  Print[{q,Precision[q]}];
  If[q==0,1,FindRoot[Erf[Sqrt[q*r]]==((2/Sqrt[\[Pi]])*Sqrt[q]),{r,2},WorkingPrecision->100,PrecisionGoal->16][[1]][[2]]]);
  Plot[F[p],{p,0,\[Pi]/4},PlotRange->{0,30},PlotPoints->5000,WorkingPrecision->100]

but that worries me. Maybe there is a better way if the precision is essential for the FindRoot

POSTED BY: Bill Nelson

Your equation is trivial to solve symbolically:

sol = Solve[Sqrt[p*r]/InverseErf[(2/Sqrt[\[Pi]])*Sqrt[p]] == 1, r]
Plot[Evaluate[r /. sol[[1]]], {p, 0, Pi/4}, PlotRange -> All]
POSTED BY: Gianluca Gorni
Posted 7 months ago

As you note, InverseErf[] is defined between -1 and 1

So what is the value of p at that upper bound in your expression?

Solve[2/Sqrt[Pi]*Sqrt[p]==1,p]
(*p->Pi/4 ~= 0.78539816339744830962*)

So you are not going to get anything bigger than p==Pi/4.

Well what happens if you get really really close to but still slightly less than Pi/4?

p=Pi/4-1/1000000;
FindRoot[Sqrt[p*r]/InverseErf[(2/Sqrt[\[Pi]])*Sqrt[p]]==1,{r,1}]
(*{r->15.7869}*)
p=Pi/4-1/1000000000000;
FindRoot[Sqrt[p*r]/InverseErf[(2/Sqrt[\[Pi]])*Sqrt[p]]==1,{r,1}]
(*{r->32.9328}*)
p=Pi/4-1/1000000000000000000;
FindRoot[Sqrt[p*r]/InverseErf[(2/Sqrt[\[Pi]])*Sqrt[p]]==1,{r,1},WorkingPrecision->64]
(*{r->50.26192019294443484109321375586527076626636935870604968790240329}*)
p=Pi/4-1/1000000000000000000000000;
FindRoot[Sqrt[p*r]/InverseErf[(2/Sqrt[\[Pi]])*Sqrt[p]]==1,{r,1},WorkingPrecision->64]
(*{r->67.66701125584564001396349443097703457649792304560596948761452092}*)
p=Pi/4-1/1000000000000000000000000000000000000000000000000;
FindRoot[Sqrt[p*r]/InverseErf[(2/Sqrt[\[Pi]])*Sqrt[p]]==1,{r,1},WorkingPrecision->64]
(*r->137.5829498581858273409042572545363711576114553405980763677426141*)

and those aren't even REALLY pushing the boundaries of precision yet, where you need to start pulling out all the tricks needed to get accurate answers to much more challenging problems.

POSTED BY: Bill Nelson
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract