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$?
Some additional background info:
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)$
Best Answer
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.
See also: http://mathdl.maa.org/images/upload_library/4/vol6/Sarra/Chebyshev.html