# Interpolation of Gaussian function - minimize relative error

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.

• @Amzoti, seen the page before. Did I explain the problem in a bad way? Cause the link is completely irrelevant. I have already tried the method suggested. The fact that there is a mention of my function there doesn't mean that this is what I need. – petajamaja Jan 4 '14 at 4:53
• @Amzoti, or did you mean the "unless the points are very close to each other" part? So, I should increase the number of points? – petajamaja Jan 4 '14 at 4:56
• So it is relative error from the start on, that is, the interpolation $p$ should satisfy $(1-ϵ)e^{-x^2}\le p(x)\le(1+ϵ)e^{-x^2}$? Resp. $e^{x^2}p(x)-1\in(-ϵ,ϵ)$. This gives a pretty hefty magnification factor that would explain your problem in finding the approximation. – LutzL Jan 4 '14 at 8:30
• @LutzL, yes, it is relative error. The formula is $-\epsilon$ ≤ (phi/f-1)(x) ≤ $\epsilon$, where phi(x) is the function given by interpolation, and f(x) is the function I needed to interpolate. – petajamaja Jan 4 '14 at 11:33
• I mean, phi(x) is the same thing as p(x) in your example – petajamaja Jan 4 '14 at 11:49

## 1 Answer

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

• Maple crushes every time I try this... – petajamaja Jan 4 '14 at 18:02
• try what? Verifying that the relative error at 15 is really that small, constructing the actual polynomial,...? If you have it, try the text interface. It may still be faster and more stable than the graphical interface. – LutzL Jan 4 '14 at 18:46
• Basically, I have all the methods already implemented in a standard worksheet (.mw), and I have just tried to modify the variable containing the number of nodes (to be in the range you proved). It keeps evaluating forever when I set the number higher than 150. – petajamaja Jan 4 '14 at 19:01
• Could you explain what you need the polynomial for? Numerical stability would be greatly increased if you were to compute the approximation via the inverse, i.e., approximations of the denominator in $1/e^{x^2}$. I added another variant that may or may not be suitable for your purposes. – LutzL Jan 4 '14 at 19:24
• Ahm... Basically, this is our university assignment, and we are mainly restricted in the methods we should use. I am, however, very much interested in your other algorithm, would you please add it to the answer? P.S. The modification you suggested is awesome! I'll definitely try it and report the results. – petajamaja Jan 4 '14 at 19:42