Interpolation of Gaussian function - minimize relative error

Question

I have been trying to interpolate the function $e^{-x^2}$ on interval [-15,15] using standard methods like Lagrange or Newton interpolation for over a month. The goal is for it to be bound by $-\epsilon$ and $\epsilon$.

No result reached. So I switched to other type of nodes distribution : tried Chebyshev nodes, skewed Chebyshev nodes (see e.g. this link), cosinus nodes... No result reached.

Maybe the skew coefficient was not high enough? I tried to increase the power to which I raise the initial nodes, so they were skewed to the edges even more. No result with different parameters.

Aware of the fact that the higher number of nodes (well, I used about 80-100 nodes in the previous attempts) doesn't necessary lead to higher precision ( because of the Runge's phenomenon), I decided to split the intervals which I work on, and interpolate each of them separately. The edge points were included in the nodes list, because I needed the interpolation polynomial to remain continuous. With the help of this method, I managed to reach the desired error, but absolute, not relative, and also was forced to increase the number of nodes (80 - 45 - 80).

And right now I am stuck. Could anyone please help, provide an explanation of what I am doing wrong or give me links on any other algorithms (maybe a good explanation of Hermite interpolation or so), so that I am able to implement those methods in Maple?

P.S. Important note : I am allowed to use only cubic spline. Of course the desired result was reached easily with the spline degree of 4 or 5.

Answer 1

Comment (but probably not solution) too long for the comment field: The series of the exponential gives an alternating series for the given function

$$e^{-x^2}=1-x^2+\tfrac12 x^2-\tfrac16x^4+\tfrac1{24}x^8+\dots+\tfrac{1}{k!}(-x^2)^k+\dots$$

By the quantitative statements of the Leibniz criterion,

$$e_{2n-1}(-x^2)\le e_{2n+1}(-x^2)\le e^{-x^2}\le e_{2n}(-x^2),$$

where $e_n$ is the $n$-th partial sum of the exponential series and the estimate is true for $x^2<2n$. This gives an absolute error of

$$\frac1{(2n)!}x^{4n}\left(1-\tfrac{x^2}{2n+1}\right)\le\left|e_{2n-1}(-x^2)- e^{-x^2}\right|\le\frac1{(2n)!}x^{4n}$$

which shows that in terms of polynomial approximation, one can not do much better than the partials of the exponential series. To get a small relative error, one now needs to chose $n$ large enough to get

$$\frac{e^{15^2}(15)^{4n}}{(2n)!}<ϵ$$

or about

$$2n(\ln(2n)-\ln(225e))>225+|\ln(ϵ)|$$

Which means that you need $n=406,...,419$ to get $ϵ=10^{-4},...,10^{-18}$

---

Using arithmetic rules of the exponential, one notices that $e^{-x^2}=\left(e^{-(\tfrac{x}{m})^2}\right)^{m^2}$, for instance with $m=15$, and the relative error $ϵ$ is guaranteed if the partial sum $e_n$ is taken such that the relative error of $e_n(-x^2)$ to $e^{-x^2}$ on $[-\tfrac{15}m,\tfrac{15}m]$ is smaller than $\tfrac{ϵ}{2m^2}$. With $m=15$ this requires the easier to control inequality

$$\frac{e^1}{(2n)!}<\tfrac{ϵ}{450}$$

which gives the usual precisions for $n=6,...,12$.

$e_{811}(-15^2)$ on Magma has a gross error, instead of ...e-98 it shows 8e73.

$e_{21}(-\tfrac{x^2}{225})^{225}$ has relative error 5e-19 at $x=15$.

---

If only basic operations are allowed, then we go back to basic divide and conquer or half-and-square in this case. Take in the above calculus $m$ a nice power of $2$, $m=16$ could work, but let's take $m=32$. Then the error estimate is still largest at the interval end, and the condition at the even larger $x=16$ is

$$\frac{e^{\frac14} \left(\frac14\right)^{2n}}{(2n)!}<\tfrac{ϵ}{2048} \iff \frac{\sqrt[4]e}{2^{4n-11}(2n)!}<ϵ$$

Trying out some small values leads to $n=5$ as a likely candidate, and indeed, $e_9(-\tfrac{15^2}{1024})\cdot\exp(15^2)-1\approx -1e-10$

$n=4$ with $e_7$ works as well with an actual error of about $2e-7$, with error bound of actually $1e-6$.

With less than 30 operations, replace 7 by 9 for higher accuracy:

xx=-x*x/1024;
a=xx;
e=1+a;
for k=2 to 7 do a*=xx/k; e+=a; end for;
for k=1 to 10 do e=e*e; end for;
return e

or with a division

xx=x*x/1024;
a=xx;
e=1+a;
for k=2 to 7 do a*=xx/k; e+=a; end for;
e=1/e;
for k=1 to 10 do e=e*e; end for;
return e