# Bell-shaped polynomial over a limited domain

The function $f(x) = e^{-x^2}$ has a bell-shaped peak at $x=0$ and then approaches an asymptote at $y=0$. I need to achieve a similar result, but with a polynomial function.

I can use a series approximation of $e^{-x^2}$: $$1 - x^2 + \frac{x^4}2 - \frac{x^6}6 + \cdots,$$ but I need my function to be accurate over the domain of $[-20000, 20000]$, which for this method requires a polynomial of a very high degree.

How can I create a polynomial function of the lowest degree where the area under the curve from $x=-2$ to $2$ accounts for 99% of the total area on the domain of $x= -20000$ to $20000$?

My intention is to investigate the feasibility of a type of additive sound synthesis wherein the sine waves are not stored as discrete waves, but instead as a spectrogram that represents the amplitude of each partial within the signal and converting that directly into the resulting waveform, which allows for evaluating the partials by only calculating the sine of one angle, no matter how many partials there are.

So the first thing I need is a function that peaks at a specified input value and quickly drops to near-zero, which I'll call $R(f)$. Each partial will also be capable of changing amplitude and frequency over time, so each partial, p will be represented by $A_p(t)R(f_p(t))$, where $A_p(t)$ represents the amplitude of the partial at time t and $f_p(t)$ represents its frequency.

Now what we do is add several of these together to get a function where each x-value corresponds to a frequency and the y-values correspond to its amplitude: $y(f) = \sum_{p=1}^n R_p(f)$

Finally, to obtain the waveform, multiply by $\sin(2\pi f t)$ and integrate with respect to f and sample to result at distinct values of t to recreate the waveform: $\int_0^{20000} y(f) sin(2\pi f t) df$ (20000Hz is about the limit of human hearing)

Some interesting things happen if $R(f)$ can be represented as a polynomial of f. Namely, y(f) will also be a polynomial of f of the same degree. Then if we instead approximate its multiplication by $sin(2\pi f t)$ with another polynomial function, evaluating the integral is trivial. There are some optimizations that can be applied to make this not so computationally demanding. For example, if we assume that the amplitude and frequency of the partials remain fairly constant, the integral can instead be evaluated by integrating by parts and only recalculating the value given by the part involving y(f) every few-hundred samples.

Now if you're wondering why I need the function to be valid from -20000 to 20000, that's because a function peaking at +20000 must still be near-zero at x=0, 20000 units to the left and when centered at x=0, in should be near-zero 20000 units to the right. So the function should extend 20000 units on either side.

Also, the requirement for 99% of the area to be within $[-2, 2]$ comes from the fact that each partial should ideally be just 1 frequency. In the case of it actually being spread across a full 4Hz, it creates an interference pattern with itself that causes its amplitude to be modulated in undesirable ways, although it can be somewhat countered by an intelligent $A_p(t)$

-
I don't think a polynomial function will approach $0$ for $x\to \infty$ – Patrick Li Oct 8 '12 at 3:33
@Patrick you are correct as far as I know - polynomials cannot have asymptotes. But I'm only evaluating the function from x=-20000 to 20000, not all the way out to infinity. So I only need it to be nearly zero on that domain, except for a sharp peak at x=0. Sorry if the wording about the asymptote in the beginning was misleading. – Wallacoloo Oct 8 '12 at 3:46
I think there is nothing useful. However, one can do very well with a rational function, indeed one of type $\frac{1}{P(x)}$, where $P(x)$ has not too large degree. – André Nicolas Oct 8 '12 at 3:47
@André Good suggestion - I'll look into it. I had only been considering polynomials because this is part of a larger scenario where I need to add curves like this, with various origins and amplitudes to each other, and then multiply it all by an approximation of $\sin(x)$ and integrate the result in a way that's optimized computationally for several thousands of these curves. $\frac{1}{P_1(x)} + \frac{1}{P_2(x)}$ is simple enough, but I'll have to see if it can be integrated trivially. – Wallacoloo Oct 8 '12 at 4:08
@Wallacoloo: Definitely not trivially! I mentioned rational functions because there is a quite large literature for those. – André Nicolas Oct 8 '12 at 4:11

Look up Chebyshev approximation. For example at http://en.wikipedia.org/wiki/Chebyshev_approximation#Chebyshev_approximation.

It's purpose is to solve exactly this sort of problem.

To compute the best approximation, you need to use something like the Remez algorithm. If you're willing to accept something somewhat less than optimal, interpolation at the zeros of Chebyshev polynomials gives a very good result with a lot less effort.

There is a Matlab add-on called Chebfun that specializes in computing high-degree polynomial approximations. See http://www2.maths.ox.ac.uk/chebfun/.

Using the series expansion you suggested is not likely to work very well. It will give you an approximation that is very good when $x$ is close to zero, but very poor when $x$ is large. Looks like you want an approximation that is equally good (in some sense) over the entire interval. This is what the Chebyshev approach gives you.

Reading your question again, it seems that you are involved with signal processing. Design of filters in DSP is often based on Chebyshev approximations. The signal processing guys refer to the Remez algorithm as the "Parks-McClellan" method.

Interpolating at equally-spaced points is definitely not the right thing to do. A famous example due to Runge shows how badly things can go wrong (and actually, its shape is quite similar to your function). You need the points to be more densely spaced at the ends of your interval, and more widely spaced in the middle. A common choice is the set of zeros of a Chebyshev polynomial. For degree $n$, these points are:

$$x_n = \cos\left( \frac{(2k-1)\pi}{2n} \right) \quad (k = 1,2,\ldots,n)$$

These points are suitable for use on the interval $[-1,1]$. You have to shift/scale them to use them on different intervals. See http://en.wikipedia.org/wiki/Chebyshev_nodes for more details.

-
Thank you for your answer. This does appear to be a much better solution. A quick test of using the Remez algorithm using points spaced 0.5 units apart from $[-8, 8]$ performs much better than a series expansion of the same degree ( i.imgur.com/CT4jW.png ), although I have a feeling that there may be a more optimal way of spacing the points than simply splitting the domain into equal intervals. Could you explain what you mean by interpolating at the zeros of the Chebyshev polynomial? – Wallacoloo Oct 15 '12 at 3:47
It's been a week since this thread's been active, so I went ahead and accepted your answer even though I don't think there's technically a single "best" answer for this question. If I try to increase the domain of the function past about +/-12 units from the origin and add more terms, I get erratic behavior. But I think that's due to the limitations of 64-bit floating point numbers and I'm in the process of implementing a fix. – Wallacoloo Oct 23 '12 at 5:25
I recommend trying Chebfun. It routinely uses polynomials with degrees in the hundreds and thousands, and sometimes even in the millions. If Chebfun can produce a good approximation, then you know that your problems are numerical, rather than theoretical. – bubba Oct 23 '12 at 23:23
There is a "best" answer for most reasonable definitions of "best". If "best" means minimizing the maximum error, then a best solution is guaranteed by the theory and given by the Remes algorithm. Chebyshev approximation (as in Chebfun) gives an answer that is provably not far from best. – bubba Oct 23 '12 at 23:26
I would have tried Chebfun long ago if I had Matlab. I signed up for the 30 day trial and it says a representative will contact me within 3 business days. I suppose you are right about there being a "best" answer though. – Wallacoloo Oct 24 '12 at 0:38

I'm not sure this question really makes sense. Perhaps I misunderstand, but I think you are solving the wrong problem.

Out at $x=\pm 19000$, the value of the function, and therefore the putative polynomial, will be vanishingly small, approximately zero. Presumably a value like "$0.000000\ldots0093$", with 150 million zeroes, is of no use, and you would rather have something that looks like "$9.3\times10^{-150000000}$", say in the form $10^n$. But to calculate the result in this form, you don't need a polynomial that approximates $e^{-x^2}$. The value of the exponent is just $$\log_{10} e^{-x^2} \\ = -x^2 \cdot \log_{10} e.$$

So all you need to calculate is $x^2$, which is trivial. You also need to have the constant $\log_{10} e$ stored ahead of time to your desired degree of precision> Then you can easily and very quickly produce the answer in the form $m\cdot10^x$ with high accuracy.

For example, consider $x=19000$.

$19000^2 = 361000000$, and $\log_{10} e \approx 0.43429448190325182765$, so \begin{align} \log_{10} e^{-19000^2} &\approx -361000000\times 0.43429448190325182765 \\& = -156780307.96707390978165 \end{align}

so therefore the value you are trying to calculate is

\begin{align} e^{-19000^2} &\approx 10^{-156780307.96707390978165} \\ & = 10^{0.96707390978165}\times 10^{-156780307} \\ & \approx 9.26987568001943595963 \times 10^{-156780307} \end{align}

No enormous polynomial approximations were required.

-
Thanks for the humorous ending there :) I added a more in-depth reason as to why I think I need a polynomial defined over that large domain into my original question. – Wallacoloo Oct 8 '12 at 4:53