Spectral Methods for SPDEs in Python and MATLAB

Spectral methods are a class of techniques used to numerically solve differential equations. They are based on expanding the solution of a differential equation in terms of a series of basis functions, often trigonometric functions (Fourier series) or orthogonal polynomials (Chebyshev, Legendre, etc.).

Stochastic partial differential equations in surface growth models | ViaDean亚图跨际 on Notion

🌵Python snippet

Spectral methods are a powerful approach for solving Stochastic Partial Differential Equations (SPDEs), leveraging the efficiency of spectral decomposition to solve differential equations with high accuracy. Python provides robust tools and libraries to implement spectral methods for SPDEs.

Overview of Spectral Methods

• Spectral Decomposition: Expands the solution of the SPDE in terms of basis functions (e.g., Fourier, Chebyshev, or Legendre polynomials).

• Stochastic Galerkin Projection: Projects the SPDE onto a reduced space using basis functions, reducing it to a system of ordinary differential equations (ODEs) or deterministic PDEs.

• Efficiency: Spectral methods achieve exponential convergence rates for smooth problems.

Libraries for Spectral Methods in Python

1. NumPy: Provides efficient array operations and Fourier transforms.

2. SciPy: Includes solvers for deterministic PDEs, useful in conjunction with spectral methods.

3. SymPy: Useful for symbolic computation, including the generation of basis functions.

4. Dedalus: A Python framework specifically designed for solving PDEs and SPDEs using spectral methods.

---

Example Workflow: Solving an SPDE using Spectral Methods

Problem Setup

We consider a 1D stochastic heat equation:

$$\frac{\partial u(x, t, \omega)}{\partial t} = \nu \frac{\partial^2 u(x, t, \omega)}{\partial x^2} + f(x, \omega),$$

where ( $\omega$) represents the stochastic component (random variable), and ( $f(x, \omega)$ ) is a random forcing term.

Steps

1. Discretize the Spatial Domain: Use a spectral basis, such as Fourier or Chebyshev polynomials, to approximate ( $u(x, t, \omega)$).

2. Expand Stochastic Term: Represent the stochastic term ( $f(x, \omega)$ ) using a series expansion, such as Karhunen-Loève expansion.

3. Galerkin Projection: Apply stochastic Galerkin projection to reduce the SPDE into a deterministic PDE system.

4. Time Integration: Use ODE solvers for temporal evolution.

---

Python Implementation

 import numpy as np
 from scipy.fftpack import fft, ifft
 from scipy.integrate import solve_ivp
 ​
 # Problem Parameters
 L = 1.0            # Length of domain
 N = 64             # Number of spectral modes
 nu = 0.1           # Diffusivity
 t_max = 1.0        # Simulation time
 dt = 0.01          # Time step
 ​
 # Spatial Grid
 x = np.linspace(0, L, N, endpoint=False)
 k = 2 * np.pi * np.fft.fftfreq(N, d=L / N)  # Wavenumbers
 ​
 # Initial Condition (example: random perturbation)
 u0 = np.random.normal(0, 1, N)
 ​
 # Random Forcing (stochastic term, e.g., white noise)
 np.random.seed(42)
 f = lambda: np.random.normal(0, 1, N)
 ​
 # RHS of the SPDE in Fourier Space
 def rhs(t, u_hat):
     u = np.real(ifft(u_hat))
     forcing = f()
     forcing_hat = fft(forcing)
     return -nu * (k**2) * u_hat + forcing_hat
 ​
 # Time Integration
 u_hat0 = fft(u0)
 t_eval = np.arange(0, t_max, dt)
 solution = solve_ivp(rhs, [0, t_max], u_hat0, t_eval=t_eval, method="RK45")
 ​
 # Transform Back to Physical Space
 u_t = np.real(ifft(solution.y.T, axis=1))
 ​
 # Visualization
 import matplotlib.pyplot as plt
 from matplotlib.animation import FuncAnimation
 ​
 fig, ax = plt.subplots()
 line, = ax.plot(x, u_t[0])
 ​
 def update(frame):
     line.set_ydata(u_t[frame])
     return line,
 ​
 ani = FuncAnimation(fig, update, frames=len(t_eval), interval=50)
 plt.show()

---

Key Concepts in the Code

1. Fourier Transform:

   • The spatial domain is represented in Fourier space.

   • fft and ifft are used to switch between physical and spectral spaces.

2. Stochastic Forcing:

   • The function f generates random forcing terms at each time step.

3. Spectral Derivatives:

   • Derivatives are efficiently computed in Fourier space using the wavenumber ( k ).

4. Integration:

   • solve_ivp integrates the time evolution of the spectral coefficients.

---

Extensions

1. Higher Dimensions:

   • Extend the approach to 2D or 3D using multidimensional FFTs (np.fft.fftn and np.fft.ifftn).

2. Other Basis Functions:

   • Replace Fourier with Chebyshev or Legendre polynomials for non-periodic domains.

3. Karhunen-Loève Expansion:

   • Use KL expansion to represent complex stochastic processes.

4. Frameworks:

   • Use Dedalus for more complex SPDEs with built-in spectral capabilities.

---

References

1. Dedalus Documentation: Dedalus Project

2. Numerical Recipes in Python: Excellent resource for spectral methods.

3. Scientific Papers: Explore SPDE-related publications for more advanced models.