Approximating the erf function

Question

I was trying to find an approximate solution to the following:

$\DeclareMathOperator\erf{erf}$

$$\frac12 \sqrt{\pi} \erf\left (\frac{x-2}{\sqrt{10}}\right) + \frac12 \sqrt{\pi} \erf \left(\frac{x+2}{\sqrt{10}}\right) = \frac25 \sqrt{\pi}$$

This naturally equals

$$\int_0^{\frac{x-2}{\sqrt{10}}} e^{-t^2} dt+ \int_0^{\frac{x+2}{\sqrt{10}}} e^{-t^2} dt = \frac25 \sqrt\pi$$

What I tried then was to use the taylor approximation and then solve the polynomial equations.. but the results I obtained were not consistently getting close to the correct solution ($x \sim 1.71$), obtained with wolfram alpha

In fact, using $\displaystyle \int_0^x e^{-t^2}dt = x - \frac{x^3}{3} + \frac{x^5}{10} - \frac{x^7}{42} + \dots$

I obtained the following approximations ($x_i$ indicates the result obtained considering the terms until $x^i$) $$x_1 \sim 1.12\ \ \ \ x_3 = -4.975 \ \ \ x_5 = 1.59 \ \ \ x_7 = -4.718 | 3.729 | 1.755 \ \ \ \ x_9 = 1.70$$

Questions

1) If I take into account the whole series, what assures me that only one solution will be found (with $x_7$, for example, I find $3$ real solutions)

2) I understand that as long as I am near $0$, the solution will be good approximation. But a priori I don't know what values $x$ is going to take, so the error can be large as large as $x$ and the Taylor approximation is useless.

3) Why is (eventually) the error going to $0$ if one takes the whole series?

I don't think I have a clear understanding of the passage between Taylor polynomial (valid only near a point $x_0$) and series (why exactly they are convergent everywhere (well, at least in the case of $e^x$) if we consider an an infinite amount of terms)

Answer 1

As you suspect, Taylor series expansions are only accurate near one point. To obtain a polynomial approximation that's good over an entire interval, the standard technique is to use Chebyshev approximation. The easiest approach is to just do interpolation at some carefully chosen nodes (the Chebyshev nodes). If you have access to Matlab, there is an add-on called Chebfun that does a very good job of constructing these sorts of approximations.

Answer 2

$\DeclareMathOperator\erf{erf}$Approximations of $\erf$ according to the Handbook of Mathematical Functions:

http://people.math.sfu.ca/~cbm/aands/page_299.htm

Answer 3

use this $$\text{Solve}\left[\frac{x \left(-2 \left(\sqrt{2 \left(\sqrt{2}+2\right) \pi } \left(\text{erf}\left(\frac{\sqrt{2-\sqrt{2}}-4}{2 \sqrt{10}}\right)+\text{erf}\left(\frac{\sqrt{2-\sqrt{2}}+4}{2 \sqrt{10}}\right)\right)-\sqrt{2 \left(2-\sqrt{2}\right) \pi } \left(\text{erf}\left(\frac{\sqrt{\sqrt{2}+2}-4}{2 \sqrt{10}}\right)+\text{erf}\left(\frac{\sqrt{\sqrt{2}+2}+4}{2 \sqrt{10}}\right)\right)\right) x^2-\left(\sqrt{2}-1\right) \sqrt{\left(2-\sqrt{2}\right) \pi } \left(\text{erf}\left(\frac{\sqrt{\sqrt{2}+2}-4}{2 \sqrt{10}}\right)+\text{erf}\left(\frac{\sqrt{\sqrt{2}+2}+4}{2 \sqrt{10}}\right)\right)+\left(\sqrt{2}+1\right) \sqrt{\left(\sqrt{2}+2\right) \pi } \left(\text{erf}\left(\frac{\sqrt{2-\sqrt{2}}-4}{2 \sqrt{10}}\right)+\text{erf}\left(\frac{\sqrt{2-\sqrt{2}}+4}{2 \sqrt{10}}\right)\right)\right)}{2 \sqrt{2}}=\frac{2 \sqrt{\pi }}{5},x\right]$$ gives x aproximately 1.7106147726585592

Answer 4

$\DeclareMathOperator\erf{erf}$Since $\erf'(x) =\dfrac{2}{\sqrt{\pi}}e^{-x^2/2} $, once you have a routine for $\erf(x)$, you can use Newton's iteration $\left(x_{n+1} =x_n-\dfrac{f(x_n)}{f'(x_n)} \right)$ to find the root of $f(x) = 0$.

Answer 5

There actually is a closed form if you allow functions implemented in a Python Library named SciPy. Notice that this function is just a special case of the Marcum Q function:

$$Q_m(a,b)=1-e^{-\frac{a^2}2}\sum_{k=0}\frac{\text P\left(m+k,\frac{b^2}2\right)}{k!}\left(\frac{a^2}2\right)^k=1-a^{1-m}\int _0^b x^v e^{-\frac{x^2+a^2}2}\text I_{v-1}(ax)dx$$

where appears the Bessel I and the Lower Normalized Incomplete Gamma “P” functions. Also, notice that it is the Non-Central Chi-Squared Distribution which is just the standard scipy.special.chndtr

$$\text{chndtr}(x,d,n)=1-Q_\frac n2(\sqrt d,\sqrt x)$$

and the special case follows after integrating using the function’s definition. Note that Wolfram Alpha’s “NonCentralChiSquareDistribution” command gives a positive sign in the integral whereas simply using the Marcum Q function definition gives the negative sign due to Wolfram Alpha’s not being the compliment unlike Wikipedia’ definition:

$$\text{chndtr}(x,p,1)=\int_0^x \frac{\cosh(\sqrt d\sqrt x)}{\sqrt{2\pi x}e^\frac{d+x}2}=\frac12\left(\text{erf}\left(\frac{\sqrt x-\sqrt d}{\sqrt 2}\right) +\text{erf}\left(\frac{\sqrt x+\sqrt d}{\sqrt 2}\right)\right)\mathop=^\text{set}p$$

Now use SciPy’s chndtrix function and a small change of variable to get a “closed form” since these functions are only in a some code databases:

$$\text{erf}\left(\frac{\sqrt x+\sqrt d}{\sqrt 2}\right)+\text{erf}\left(\frac{\sqrt x+\sqrt d}{\sqrt 2}\right)\mathop=2p\implies x=\text{chndtrix}(p,d,1) $$

Now to test it on SciPy software, which I do not know how to access. In short, the solution to your equation is a quantile of a Noncentral Chi Squared Distribution with a Noncentral Parameter of $1$. Let me finish with a few variable changes. Please correct me and give me feedback!

Here is the solution:

$$\text{erf}\left(\sqrt \frac x2-\sqrt \frac d2\right) +\text{erf}\left(\sqrt \frac x2+\sqrt \frac d2\right)= \frac2{\sqrt\pi}\frac {2\sqrt\pi}5=2p\implies x=\text{chndtrix}(p,d,1)$$

$$\mathop\iff ^{x\to a x}$$

$$ \text{erf}\left(\sqrt \frac {ax}2-\sqrt \frac {d}2\right) +\text{erf}\left(\sqrt \frac {ax}2+\sqrt \frac {d}2\right)=\frac 45=2p\implies ax=\text{chndtrix}(p,d,1) $$

set:

$$\sqrt{\frac{a x}2}=\frac x{\sqrt {10}},\sqrt{\frac d2}=\frac 2{\sqrt{10}},2p=\frac45\implies a=\frac 15,p=\frac 25,d=\frac45$$

Therefore:

$$\boxed{\frac12 \sqrt{\pi} \text{erf}\left (\frac{x-2}{\sqrt{10}}\right) + \frac12 \sqrt{\pi} \text{erf }\left(\frac{x+2}{\sqrt{10}}\right) = \frac25 \sqrt{\pi}\implies x=5\,\text{chndtrix}\left(\frac25,\frac45,1\right)}$$

which is 5 times the the value of $x$ such that a Non-central Chi Squared Distribution, with non centrality parameter of $1$ and $\frac45$ degrees of freedom, has a probability of $\frac25$.