I have to solve a differential equation involving the error function Erf, but my code is not able to evaluate it in the range I need.

I think I can formulate my problem in the following simplified terms: let's consider the following functions

VV[x_]:=Exp[-2 x^2]

f[r_]:=cc[r^-2]/Sqrt[cc[r^-2]^2 + (1 + r) VV[r^-2]^2]

I have to evaluate the function f[r] for small values of r. However, for r<.2 I get the following message

Power::infy: "Infinite expression 1/0. encountered"

due to the presence of error functions. I tried to increase the precision but without success.

Do you have any suggestions to overcome this problem? Thank you very much for your help.

  • 3
    $\begingroup$ Try cc[x_] := -Erfc[x] instead. $\endgroup$ – user484 May 14 '15 at 16:36
  • $\begingroup$ thank you for the suggestion, I will try immediately $\endgroup$ – user9994 May 14 '15 at 16:36
  • $\begingroup$ Indeed, it works perfectly. Thank you very much! $\endgroup$ – user9994 May 14 '15 at 16:45

For large $x$, the value of $\operatorname{erf}(x)$ approaches $1$, so even if you were able to evaluate it you would encounter a catastrophic loss of precision when you subtracted $1$ from it.

For this reason, implementations typically also provide the complementary function $$\operatorname{erfc}(x)=1-\operatorname{erf}(x)$$ designed to provide better precision in these situations. Thus you can define

cc[x_] := -Erfc[x]

and everything works:

VV[x_] := Exp[-2 x^2]

(* bad *)
cc[x_] := Erf[x] - 1
f[r_] := cc[r^-2]/Sqrt[cc[r^-2]^2 + (1 + r) VV[r^-2]^2]

(* good *)
cc2[x_] := -Erfc[x]
f2[r_] := cc2[r^-2]/Sqrt[cc2[r^-2]^2 + (1 + r) VV[r^-2]^2]

Plot[{f2[r], f[r]}, {r, 0, 2}, PlotStyle -> {Thickness[0.01], Automatic}]

enter image description here

If you wanted to evaluate $\operatorname{erf}(x)+1$ for large negative $x$, you would have the same problem because $\operatorname{erf}(x)\to-1$ as $x\to-\infty$, and then you would want to use $$\operatorname{erf}(x)+1=-\operatorname{erf}(-x)+1=\operatorname{erfc}(-x)$$ instead.

| improve this answer | |
  • $\begingroup$ Many thanks, your explanation is extremely clear and the problem is now solved. $\endgroup$ – user9994 May 15 '15 at 7:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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