Efficiently Solving Differential Equations with FFT Techniques
Written on
Chapter 1: Introduction to Fourier Transform in Differential Equations
In our previous discussion, we explored how Fourier methods enable the replacement of differentiation with multiplication. This principle can be harnessed to numerically solve partial differential equations with relative ease. Today, we will delve deeper into this concept.
Daily Dose of Scientific Python
Here, we tackle both common and unique challenges using libraries like Numpy, Sympy, SciPy, and matplotlib.
The Concept
Let's consider the task of solving a differential equation defined as follows:
for a function ( f ) over the domain ( [a,b] ). We will impose the boundary conditions ( f(pm a)=0 ). While this is technically an ordinary differential equation, the methodology can be extended to multidimensional cases.
By expressing the equation as:
and similarly for ( f(x) ), the differential equation transforms into:
The use of complex exponentials provides an orthogonal basis for function space, allowing us to compare coefficients:
From here, we can derive the Fourier component of ( f ):
The essential takeaway is that with this formula, we can perform the inverse (discrete) Fourier transformation to obtain our solution ( f(x) ).
However, caution is needed! If ( f=0 ), what does this imply? Indeed, ( f ) can be zero, which indicates a constant term in the Fourier series:
Therefore, we can assign ( c_0 ) a value of our choice. After executing the inverse Fourier transform, the resultant function may be offset by a constant value, which we can easily adjust. It's important to note that while a second-order differential equation allows for two boundary conditions, we only appear to have one choice here: the offset. However, by applying the Fourier technique, we implicitly impose periodicity, which means selecting the offset effectively determines both boundary conditions!
Now, let's examine a specific coding example. We will select:
How did I generate this example? I constructed it by inserting the function ( f(x)=exp(-x^2) ) into the differential equation, leading to ( f(x) ). I then pretended that I was unaware of the solution’s form and set out to solve the equation. This is a common and effective strategy for creating examples with analytical solutions, useful for verifying code functionality.
The Implementation
Implementing this in Python with Numpy is straightforward:
import numpy as np
import matplotlib.pyplot as plt
from numpy.fft import fft, fftfreq, ifft
N = 201 # number of grid points
a = 10
x = np.linspace(-a, a, N)
# The exact known solution for later comparison:
f_exact = np.exp(-x**2)
# The right-hand side of the ODE
rho = (4*x**2 - 2)*np.exp(-x**2)
# The right-hand side in Fourier space:
rho_k = fft(rho)
# The coordinate axis in Fourier space:
k = fftfreq(N) * 2*np.pi * N / (x[-1] - x[0])
# Prepare the solution in Fourier space. We are
# setting the constant term in the Fourier series to 0:
f_k = np.zeros_like(rho_k)
# Avoid division by zero:
mask = k**2 != 0
f_k[mask] = -rho_k[mask] / (k**2)[mask]
# Perform the inverse FFT transformation:
f = ifft(f_k)
The only new element compared to our last article is the utilization of the fftfreq function in Numpy, which simplifies the construction of the "frequency" axis.
The resulting function ( f ) is real again. Here’s how it appears:
Clearly, this solution does not yet meet the boundary conditions ( f(pm a)=0 ). This is expected because we set the constant term of the Fourier series arbitrarily to zero. The constant term represents the mean value of the function over the interval. We simply need to adjust our result to ensure it meets the boundary condition:
f -= f[0]
And that concludes our discussion for today. In the next installment, we will explore how to extend this methodology to multiple dimensions, where the true strength of this technique becomes apparent.
Chapter 2: Practical Applications of FFT in PDEs
This video titled "Solving PDEs with the FFT [Matlab]" offers a practical demonstration of applying FFT techniques to solve partial differential equations. It serves as a valuable resource for understanding the concepts discussed.
In the second video, "Solving PDEs with the FFT [Python]," viewers can see a Python implementation of FFT methods in action, complementing our previous examples and providing further insights into the topic.