First load the package for plotting.
using Pkg
Pkg.activate("../../.")
using PyPlot, FFTW
Create the physical grid $x$. For a domain of length $L_x$ with $n_x$ grid points the grid spacing is $\delta=L_x/n_x$. Our grid points then are: $x = 0,\delta, 2\delta,\dots,L_x-\delta$. Notice, that the final grid point does not correspond to $x=L_x$ but to $x=L_x-\delta$. This is because of periodicity: $f(L_x)=f(0)$ and, therefore, the grid point $x=L_x$ is the same as $x=0$.|
Thus, our grid points for an array with elements $x_j = (j-1)\delta$, where $j=1,\dots,n_x$.
nx = 32
Lx = 2.0*pi
dx = Lx/nx
x = 0:dx:Lx-dx;
Let's define a function, e.g., $f(x) = \sin(2x)+\frac1{2}\cos(5x)$. The approximation of this function on our $x$-grid is the array with elements $$f_j \equiv f(x_j),\text{ for }j=1,\dots,n_x.$$
Let us see how well the array $f_j$ approximates function $f(x)$.
f = @. sin(2*x) + 0.5*cos(5*x) # a test function f(x)
x_cont = range(0, Lx, length=1001) # a much finer x-grid used for plotting the "actual" f(x)
f_cont = @. sin(2*x_cont) + 0.5*cos(5*x_cont) # the same function f(x) defined on the finer grid
# plot f(x) and f_j; for f(x) we just plot f on the finer grid, f_cont
figure(1, figsize=(8, 3) )
plot(x_cont, f_cont, linestyle="solid", label=L"f(x)") # plots the abs value of fft of f(x)
plot(x, f, marker="*", color="r", linestyle="none", label=L"f_j") # plots the abs value of fft of f(x)
legend()
xlim(0, Lx)
xlabel(L"x")
ylabel(L"f(x)");
Remark: Do $n_x=32$ grid points capture the behavior of $f(x)$ good enough? How would a grid with $16$ points look? A rule-of-thumb is that we need about $4$ grid points per wavelength of the highest harmonic. In our case $f(x)$ includes a $\sin(5x)$ term, oscillates 5 times within our domain. Therefore our rule-of-thumb would argue that we capture the behavior of $f(x)$ if we use at least $20$ grid points.
The derivative of the above function is $df/dx = 2\cos(2x)-\frac5{2}\sin(5x)$. Let us pretend that we do not know the derivative and instead try to compute it numerically.
(When we are coding something new we should always try to start with something we know what the answer should be so we can compare our numerical results. This step may seem obvious, but, nonetheless, it is often overlooked by most people...)
The simplest approximation for the array of the derivative $(df/dx)_j$ is the, so called, second-order centered finite differences: $$\left(\frac{df}{dx}\right)_j = \frac{f_{j+1}-f_{j-1}}{2\delta},\text{ for }j=1,\dots,n_x,$$ where it is understood that $f_{0}=f_{n_x}$ and $f_{n_x+1}=f_1$ due to periodic boundary conditions.
Let's see how this approximation compares with the actual derivative.
# the analytical expression for df/dx
dfdx_an = 2*cos.(2*x) - 2.5*sin.(5*x)
# the analytical expression for df/dx computed on the finer grid (for plotting purposes)
dfdx_an_cont = 2*cos.(2*x_cont) - 2.5*sin.(5*x_cont)
# the finite differences approximation
dfdx_fd = zeros(nx)
for j in 2:nx-1
dfdx_fd[j] = (f[j+1]- f[j-1])/(2*dx)
end
# for the j=1 and j=nx grid points we use the fact that f(x) is periodic!
dfdx_fd[1] = (f[2] - f[nx]) /(2*dx)
dfdx_fd[nx] = (f[1] - f[nx-1])/(2*dx)
# Plot (df/dx)_j and compare with the actual df/dx
fig, axs = subplots(nrows=2, ncols=1, figsize=(8, 6) )
sca(axs[1])
plot(x_cont, dfdx_an_cont, label="analytical")
plot(x, dfdx_fd, marker="d", color="g", linestyle="none", label="numerical")
xlabel(L"x")
ylabel(L"df / dx")
legend(loc="upper right", fancybox="true")
xlim(0, Lx-dx)
sca(axs[2])
plot(x, abs.(dfdx_an - dfdx_fd), marker="d", label="|numerical-analytical|")
legend(loc="upper left", fancybox="true")
xlabel(L"x")
ylabel("absolute error")
xlim(0, Lx-dx);
We see that the error at each grid-point is of order $10^{-1}$. If we increase the number of grid-points, of course, the error will decrease.
For second-order finite differences the error should decrease as $n_x^2$. This means that if we increase the numbers of grid-points by a factor of $10$ then the error should decrease by a factor of $10^2$. This is OK, but we can do much better using spectral methods.
Exercise: Confirm the above assertion. Show (merely by plotting) that the error for the finite-difference approximation of the derivative decreases as $n_x^2$?
To use spectral methods we need to remind ourselves what a Fourier series is.
Any periodic and square-integrable function defined on the domain $0\le x\le L_x$ can be expanded as a Fourier series: $$ f(x) = \sum_{k} \hat{\;f}_{k}\,e^{i k x},$$ where the wavenumbers $k = \frac{2\pi}{L_x}[0,\pm1,\pm2,\dots]$. The coefficents $\hat{\;f}_{k}$ are $$\hat{\;f}_{k} = \frac1{L_x}\int_{0}^{L_x} f(x)\,e^{-i k x}\,dx.$$
By "square-integrable" we mean that $\int_0^{L_x} \left|f(x)\right|^2 dx <\infty$. Note that $f$ does not have to be differentiable. In fact, it can even be discontinuous within our domain. However, in this tutorial we will only use differentiable functions for our examples.
The derivative is trivially expressed as a Fourier series as: $$ \frac{df}{dx} = \sum_{k} i k \hat{\;f}_{k}\,e^{i k x}.$$ (Just differentiate the right-hand-side of the Fourier series for $f(x)$.)
Our function $f$ is, conveniently, already expanded as a Fourier series so we can read off the coefficients $\hat{\;f}_{k}$: $$\hat{\;f}_{k}=\begin{cases} \frac1{2i} \quad\text{if $k=4\pi/L_x$,} \\ -\frac1{2i} \quad\text{if $k=-4\pi/L_x$,} \\ \frac1{4} \quad\text{if $k=10\pi/L_x$,} \\ \frac1{4} \quad\text{if $k=-10\pi/L_x$,} \\ 0 \quad\text{otherwise.} \\ \end{cases}$$ This mean that for our example the Fourier series is not an infinite sum but rather it terminates since after some wavenumber all coefficients are zero.
Exercise: Use the coefficients above and compute the Fourier series for $\sum_{k} i k \hat{\;f}_{k}\,e^{i k x}$ (i.e. sum all the terms). Do you recover the analytical expression for $df/dx$?
The wavenumbers that take part in the sum for the Fourier series of a function $f(x)$ form an array: $k_m = (2\pi/L_x) m$, with $m=0,\pm1,\pm2,\dots$. However, if our function $f(x)$ is restricted only on the values of the $x$-grid, $x_j=(j-1)\delta=(j-1)L_x/n_x$, $j=1,\dots,n_x$, then only $n_x$ of all wavenumbers are independent. Thus the Fourier series truncates from an infinite sum to a finite sum. To see first that consider: $$ f(x_j) = \sum_{m=-\infty}^{+\infty} \hat{\;f}_{k} \exp{\Big[ i \underbrace{\frac{2\pi}{L_x} m}_{= k_m} \; \underbrace{(j-1)\frac{L_x}{n_x}}_{=x_j}\Big]}$$ But the exponential inside the sum above remains unchanged if we replace $m\mapsto m+n_x$, i.e., $$ \exp{\left[ 2\pi i \frac{(m+n_x) (j-1)}{n_x}\right]} = \underbrace{\exp{\left[ 2\pi i (j-1)\right]}}_{=1\text{ since }j\in\mathbb{Z}}\exp{\left[ 2\pi i\frac{m (j-1)}{n_x}\right]} = \exp{\left[ 2\pi i \frac{m (j-1)}{n_x}\right]}. $$
Therefore we only need to take $n_x$ different wavenumbers in the sum. It does not really matter which ones we choose to have, given that we include $k_x=0$. Conventionally, we take $n_x$ to be even and choose: $$k = \frac{2\pi}{n_x}\left[0,1,2,\dots,\frac{n_x}{2},-\frac{n_x}{2}+1,-\frac{n_x}{2}+2,\dots,-2,-1\right].$$
If our example function $f$ was not as simple we couldn't not just read off the values of the coefficients $\hat{\;f}_{k}$ from $f$. In that case we could go back to the definition and compute the coefficients $\hat{\;f}_{k}$ through the integral above (or by approximating the integral numerically). Fortunately though, even in those cases we don't have to do that. There exist the Fast Fourier Transform (fft
) algorithm (which is included in Matlab, Python, Julia, etc). This algorithm is very well optimized, i.e., it is much faster than actually approximating the integral numerically. We can obtain the coefficents $\hat{\;f}_{k}$ using fft
as:
$$
\mathtt{fft}(f) = \hat{\;f}_{k}\,n_x.
$$
Above by $\mathtt{fft}(f)$ we mean performing fft
on the array with elements $f_j$. The multiplication with $n_x$ on the right-hand-side above is due to the normalization that fft
algorithm uses. Furthermore, with the inverse, ifft
we can compute $f$ from $\hat{\;f}_{k}$.
$$
\mathtt{ifft}(\hat{\;f}_{k}) = f\big/n_x.
$$
With the coefficients $\hat{\;f}_{k}$ in hand we compute $df/dx$ as $$\frac{df}{dx} = \mathtt{ifft}\left[ i k \mathtt{fft}(f)\right] .$$
Let us see how well the numerical approximation of the derivative using spectral methods compares with the actual analytical expression for $df/dx$.
# the k_x wavenumbers in the order that the FFT algorithm assumes
k = 2.0*pi/Lx * cat(0:nx/2, -nx/2+1:-1, dims=1);
fh = fft(f)/nx # the FFT of f(x)
dfdx_num = real(ifft(im*k.*fh)*nx); # the numerical approximation for df/dx
Note that after we used ifft
we also took the real part. We do that to remove any imaginary parts that the algoritm might introduce (which should be of order $10^{-16}$).
First we plot the coefficients $\hat{\;f}_{k}$ obtained using the fft
algorithm. Do we see power at the wavenumbers we were expecting?
figure(1, figsize=(8, 3) )
plot(k, abs.(fh), marker="d", markersize=10, linestyle=":", linewidth=0.5) # plots the abs value of fft of f(x)
xlim(-12, 12)
xlabel(L"k_x")
ylabel(L"|\hat{f}_{k_x}|");
Now let's compare the numerical approximation of $df/dx$ with the analytical one.
fig, axs = subplots(nrows=2, ncols=1, figsize=(8, 6) )
sca(axs[1])
plot(x_cont, dfdx_an_cont, label="analytical")
plot(x, dfdx_num, marker="*", color="r", linestyle="none", label="numerical")
xlabel(L"x")
ylabel(L"df / dx")
legend(loc="upper right", fancybox="true")
xlim(0, Lx-dx)
sca(axs[2])
plot(x, abs.(dfdx_an - dfdx_num), marker="d", label="|numerical-analytical|")
legend(loc="upper left", fancybox="true")
xlabel(L"x")
ylabel("absolute error")
xlim(0, Lx-dx);
Wow!
See how much better accuracy we got now? We got an error of (nearly) machine precission. And we obtained this just with only $n_x=32$ points. Try doing that with finite differences...
(The error with finite differences using $n_x=10^{20}\approx 1000000$ points is still about $10^{-10}$, i.e., $10^4$ larger compared to spectral methods!)
Exercise: Can we perhaps push spectral methods even further? What about if we had $n_x=16$ points? How far can we decreasee the resolution and still evaluate $df/dx$ with such good accuracy?
Let's apply the same principles on a two-dimensional domain for periodic functions $f(x,y)$, i.e., function with the property $f(x+L_x,y)=f(x,y+L_y)=f(x,y)$. Same principle apply; now the Fourier series is: $$ f(x,y) = \sum_k\sum_l \hat{\;f}_{k,l} \, e^{i(k x+l y)},$$ with $k=(2\pi/L_x)[0,\pm 1,\pm2,\dots]$ and $l=(2\pi/L_y)[0,\pm 1,\pm2,\dots]$.
First, we create the physical grid $x$, $y$. It is instructive to take different domain lengths, $L_x$ and $L_y$ and different grid points in each direction, $n_x$, $n_y$. This way make sure that we treat each dimension correctly.
nx, ny = 64, 32 # number of grid points
Lx, Ly = 2.0*pi, 1.0*pi # size of the domain in each direction
# constructing the physical grid (x,y)
dx, dy = Lx/nx, Ly/ny
x = 0:dx:Lx-dx
y = 0:dy:Ly-dy
X = zeros(nx,ny)
Y = zeros(nx,ny)
for j in 1:ny, i in 1:nx
X[i,j], Y[i, j] = x[i], y[j]
end
# constructing the wavenumber grid (k,l)
k = 2.0*pi/Lx * cat(0:nx/2, -nx/2+1:-1, dims=1);
l = 2.0*pi/Ly * cat(0:ny/2, -ny/2+1:-1, dims=1);
K = zeros(nx, ny)
L = zeros(nx, ny)
for j in 1:ny, i in 1:nx
K[i, j] = k[i]
L[i, j] = l[j]
end
Now let's define a two-dimensional function: $$f(x,y) = \cos\left(3\frac{2\pi}{L_x} x\right)\sin\left(2\frac{2\pi}{L_y} y\right).$$ Let's compute $\partial f/\partial y$ numerically and compare with the analytical expression.
k0 = 2.0*pi/Lx # the fundamental x-wavenumber
l0 = 2.0*pi/Ly # the fundamental y-wavenumber
mx, my = 3, 2
f = cos.(mx*k0*X).*sin.(my*l0*Y)
fh = fft(f)
dfdy_an = my*l0*cos.(mx*k0*X).*cos.(my*l0*Y)
dfdy_num = real(ifft(im*L.*fh));
Exercise: Plot the two-dimensional power spectra of $f$: $|\hat{\;f}_{k,l}|$. Is it how you expexted it to be? (Regarding normalization, the coefficients $\hat{\;f}_{k,l} = \mathtt{fft}(f)\big/(n_x n_y)$.)
Now, as we did before, let's plot the numerical and analytical $\partial f/\partial y$ and compare. If all is good, we should expect machine precision error...
levs = -4:1:4 # the contour levels for the plots of \partial f/\partial y
# contour plot of the analytical \partial f/\partial y
figure(figsize=(8, 3))
contourf(X, Y, dfdy_an, levs) # notice that contour levels are prescribed
xlabel(L"$x$")
ylabel(L"$y$")
title(L"analytical $\partial f/\partial y$")
xlim(0, Lx-dx)
ylim(0, Ly-dy)
clim(-4,4)
colorbar(ticks=-4:2:4)
# contour plot of the numerical approcimation to \partial f/\partial y
figure(figsize=(8, 3))
contourf(X, Y, dfdy_num, levs)
xlabel(L"$x$")
ylabel(L"$y$")
title(L"numerical $\partial f/\partial y$")
xlim(0, Lx-dx)
ylim(0, Ly-dy)
clim(-4,4)
colorbar(ticks=-4:2:4)
# contour plot of the absolute error
figure(figsize=(8, 3))
contourf(X, Y, abs.(dfdy_an - dfdy_num))
xlabel(L"$x$")
ylabel(L"$y$")
title("absolute error = |analytical - numerical|")
xlim(0, Lx-dx)
ylim(0, Ly-dy)
colorbar();
What's your favorite number? Mine is $10^{-16}$; this is the machine precision with 64-bit arithmetic. When I see that probably I know I did everything correct. (Often if I see just plain $0$ the usually something's going wrong.)
Do we get machine precision if we calculate the relative error that we made with our approximation?
rel_error = abs.(maximum(abs.(dfdy_an)) - maximum(abs.(dfdy_num)))/maximum(abs.(dfdy_an))
1.1102230246251565e-15
Sure we do...