How to compute the Fourier coefficients of $\exp(f(x))$ given the Fourier coefficients of $f(x)$

Question

In an attempt to solve a certain differential equation numerically with the Galerkin method I bumped into this problem:

Given the real function $f(x) = \sum_{n = -N}^{N} a_n \, \exp(inx) $ compute the Fourier coefficients with indices between $-N$ and $N$ of the function $\exp(f(x))$.

Of course the coefficients can be computed to any required accuracy by brute force numerical methods such as Gaussian Quadrature. In terms of computation time the $O(N^2)$ operations needed by the Gaussian quadrature are however to expensive in this case. What I am looking for is therefore a closed formula or algorithm for the coefficients such that the $2N+1$ coefficients can be computed with $O(N \log N)$ operations or less.

Here is what I have tried:

Let $b_p$ denote the $p$'th Fourier coefficient of $\exp(f(x))$. For notational convinience we will consider the coefficient $b_{-p}$, which can be calculated as

$$b_{-p} = \frac{1}{2 \pi} \int_{0}^{2 \pi} \exp(f(x)) \exp(ipx) \, dx. $$

Inserting the expansion for $f$ into this expression yields

$$b_{-p} = \frac{1}{2 \pi} \int_{0}^{2 \pi} \exp\bigg( \sum_{n = -N}^{N} a_n \exp(inx) \bigg) \exp(ipx) \, dx. $$

We now rewrite this as a contour integral along the unit circle in the complex plane. One parametrization of the unit circle is $\gamma(x) = \exp(ix)$ where $0 \le x \le 2 \pi$ and by definition of a complex contour integral we therefore have

$$b_{-p} = - \frac{i}{2 \pi} \int_{\gamma} \exp\bigg( \sum_{n = -N}^{N} a_n z^n \bigg) \, z^{p-1} dz. $$

By the residue theorem of Cauchy this integral is easily calculated if the residues of the function $\exp\bigg( \sum_{n = -N}^{N} a_n z^n \bigg) \, z^{p-1}$ are known but this is where I get stuck. To get the residues I have tried to rewrite the exponential as

$$ \exp\bigg( \sum_{n = -N}^{N} a_n z^n \bigg) = \exp(a_{-N} z^{-N}) \exp(a_{-N+1} z^{-N+1}) \cdots \exp(a_{N-1} z^{N-1}) \exp(a_{N} z^{N}) $$

and then expanded each of the exponential functions as a Taylor series (the positive powers of $z$) or a Laurent series (the negative powers of $z$). This unfortunately creates a mess and I have not been able to extract a useful formula for $b_{-p}$ from it.

Answer 1

Maybe this helps: If $\exp\bigl(f(x)\bigr)=: g(x)$ then $g'(x)=f'(x)\>g(x)$. This immediately translates into $$b_n={1\over n}\sum_{k+l=n} k\, a_k \,b_l\ .$$

Answer 2

In case anybody should be interested, here is my solution to the problem.


Given the Fourier coefficients of $f(x)$ we can compute the Fourier coefficients of $f(x)^n$ for $n = 0, 1, ..., K$ in $O(N log(N))$ operations (if $K$ is negligible compared to $N$) using the fast Fourier transform. We will in the following let $a_p^{(n)}$ denote the $p$'th Fourier coefficient of $f(x)^n$.


Let us now consider the $-p$'th Fourier coefficient of $\exp(f(x))$ given by $$ b_{-p} = \frac{1}{2 \pi} \int_{0}^{2\pi} \exp(f(x)) \exp(ipx) dx .$$ Taylor expanding the exponential function and taking the integral inside the sum gives $$ b_{-p} = \sum_{n = 0}^{\infty} \frac{1}{n!} \frac{1}{2 \pi} \int_{0}^{2\pi} f(x)^n \exp(ipx) dx $$ The integral divided by $2 \pi$ is just the $-p$'th Fourier coefficient of $f(x)^n$ meaning that we can express $b_{-p}$ as $$ b_{-p} = \sum_{n = 0}^{\infty} \frac{1}{n!} a_{-p}^{(n)}$$ To obtain an approximation for $b_{-p}$ I truncate this series after $K$ terms. In my case it turns out that the series is rapidly convergent and truncating it after $K = 5$ terms gives an accurate enough solution.