Lecture 2: Numerical Methods

Jan Benáček

Institute for Physics and Astronomy, University of Potsdam

October 24, 2024

Table of contents

• Lecture overview
• 1 Plasma properties
• 2 Plasma simulations
• 3 Hands-on
• 4 Bibiography and references
• 5 Introduction to numerical methods: ODEs
• 6 Partial differential equations
• 7 References

Lecture overview

Table of Contents 1/2

• 17.10. Introduction + numerical methods summary
• 24.10. Numerical methods of differential equations — lecture + hands-on
• 31.10. Test particle approach — lecture
• 07.11. Test particle approach — hands-on
• 14.11. PIC method — lecture
• 21.11. PIC method — hands-on
• 28.11. Fluid and MHD — lecture

Table of Contents 2/2

• 05.12. Fluid and MHD — hands-on
• 12.12. Canceled
• 19.12. Canceled
• 09.01. Radiative transfer — lecture + hands-on
• 16.01. HPC computing — lecture + hands-on
• 23.01. Advanced — Smooth particle hydrodynamics method — lecture
• 30.01. Advanced — Hybrid, Gyrokinetics — lecture
• 06.02. Advanced — Vlasov and Linear Vlasov dispersion solvers — lecture

1 Plasma properties

1.1 Plasma parameters

• Temperature
• Number density
• Magnetic field (not in the figure)
  
  
• Cold plasmas
• Hot plasmas

1.2 Plasma parameters

ICF = Inertial Confinement Fusion

(Boyd and Sanderson 2003)

2 Plasma simulations

2.1 Hierarchy of plasma physics models

Kinetic description

Microscopic properties, it uses the velocity distribution function \(f(\vec{x}, \vec{v}, t)\).

Fluid description

Uses a few macroscopic quantities, averages of the distribution function (mean velocity \(v(\vec{x},t)\), pressure/temperature). Valid for exact or near thermodynamic equilibrium.

Hierarchy of plasma physics models

2.2 Validity of plasma models

Range of validity of different plasma codes based on typical magnetospheric parameters: \(n=50cm^{-3}\), \(B=50 nT\), \(T_e=T_i=100 eV\) (Winske and Omidi 1996).

2.3 A bunch of rock

https://xkcd.com/505/

3 Hands-on

3.1 Testing the computer setup

1. Login to local computers using your student accounts
2. Open command line
3. Run Jupyter Notebook/Lab: $ jupyter lab
4. Load basic Python libraries Numpy and Matplotlib
5. Plot an arbitrary analytical function (e.g., parabola, sin, cos,…)

If Jupyter Lab does not exist, run:

$ pip3 install --upgrade pip
$ pip3 install --user jupyterlab numpy scipy matplotlib
$ export PATH=$HOME/.local/bin:$PATH
$ jupyter lab

Perhaps, you want to add the PATH into your ~/.bashrc file to be loaded automatically when you log in.

OR

• Instead of (1.–2.) Install and run Jupyter Lab on your computer.

4 Bibiography and references

4.1 Bibliography

• Büchner (2023)
• Büchner, Dum, and Scholer (2003)
• Matsumoto and Omura (1993)
• Jardin (2010)
• Winske and Omidi (1996)
• Lapenta (2012)
• Umeda (2012)
• Benáček et al. (2024)

5 Introduction to numerical methods: ODEs

5.1 Introduction

• We want to solve the ordinary differential equation (ODE): \[ \frac{dy}{dt}=f(t,y), \] with the initial conditions \(y(t_0)=y_0\).

  \(\Rightarrow\) to find the value of \(y\) at a given \(t_n\): \(y_n:=y(t_n)\) (discretization).

• Idea

  1. To divide the interval \([t_0,t_N]\) in \(N\) intervals (\(N+1\) points), such as \(\Delta t=(t_N-t_0)/N\) and \(t_i=i\Delta t\) (\(n=0,1,2,\cdots,N\))
  2. To find a recursive expression relating \(y_n\) with \(\{y_{n-1},y_{n-2},\cdots\}\)

5.2 General form of an ODE

• Any ODE of order \(N\) can be written as a system of \(N\) first-order ODEs

\[ \frac{dy^{(1)}}{dt}=f^{(1)}\left(t,\{y^{(i)}\} \right)\\ \frac{dy^{(2)}}{dt}=f^{(2)}\left(t,\{y^{(i)}\}\right)\\ \cdots \qquad \cdots \\ \frac{dy^{(N)}}{dt}=f^{(N)}\left(t,\{y^{(i)}\}\right)\\ \]

\[ \] \[ \Leftrightarrow \]

\[ \qquad\frac{d\vec{y}}{dt}=\vec{f}(t,\vec{y})\qquad\qquad\qquad \nonumber\\\ \vec{y}=\begin{pmatrix} y^{(1)} \\ y^{(2)}\\\vdots \\ y^{(N)} \end{pmatrix},\qquad \vec{f}=\begin{pmatrix} f^{(1)}(t,\vec{y})\\ f^{(2)}(t,\vec{y})\\ \vdots \\ f^{(N)}(t,\vec{y}) \end{pmatrix} \qquad(1)\]

• Note that the RHSs do not have explicit derivatives.

5.3 Example

• Newton’s second law is a 2nd-order ODE that can be written as: \[ \frac{d^2x}{dt^2}=\frac{1}{m}F\left(t,\frac{dx}{dt},x\right) \]
• Let \[ y^{(1)} \equiv x, \quad y^{(2)} \equiv \frac{dx}{dt}=\frac{dy^{(1)}}{dt}. \] Then we have a system of two first-order ODEs in the form given by Equation 1: \[ \vec{y} \equiv \begin{pmatrix} x \\ \frac{dx}{dt} \end{pmatrix},\qquad \vec{f} \equiv \begin{pmatrix} y^{(2)}(t) \\ \frac{1}{m}F\left(t,y^{(1)},y^{(2)} \right) \end{pmatrix} \]

5.5 Euler’s method

• Simplest method to solve Equation 1: forward finite difference centered at \((t_i,y_i)\)

  (extrapolating the initial slope).

Euler’s method/Forward difference

\[ y_{i+1}=y_i + \Delta t f(t_i,y_i)+\mathcal{O}(\Delta t) \]

• The accuracy of this method is very low.
• More precise methods are necessary

Explicit methods

They give \(y_{i+1}\) in terms of already known values of \(y_i\) — e.g., the Euler method.

Implicit methods

They require to solve (invert) another equation (matrix) to obtain \(y_{i+1}\).

5.6 Accuracy and precision

5.7 Improved Euler’s method

• We can improve the naive Euler’s method by averaging the slopes at \((t_i,y_i)\) and \((t_{i+1},y_{i+1})\): \[ y_{i+1}^*=y_i+\Delta t\, f(t_i,y_i) \\ y_{i+1} =y_i+\Delta t\frac{f(t_i,y_i)+f(t_{i+1},y_{i+1}^*)}{2}\\ \]

• This method predicts \(y_{i+1}\) by the direct Euler’s method and then corrects it.
• Basis for the predictor-corrector schemes.

5.8 Taylor expansion and finite differences

All finite difference expressions are based on

Taylor expansion

\[ f_{i+1}=f(x_i + \Delta x)=f_i + \Delta x f'(x_i) + \frac{\Delta x^2}{2!}f''(x_i) + \frac{\Delta x^3} {3!}f'''(x_i)+\cdots \]

Better accuracy can be obtained by considering more points \[ f_{i+1}-f_{i-1}= 2\Delta xf' + \frac{\Delta x^3}{3}f'''+\cdots \]

Central differences or 3 points formula

\[ f'(x_i)=\frac{1}{2\Delta x}\left[f_{i+1}-f_{i-1}\right]+\mathcal{O}(\Delta x^2) \]

And similarly, using \(f_{i+2}-f_{i-2}\),

5 points formula

\[ f'(x_i)=\frac{1}{12\Delta x}\left[f_{i-2}-8f_{i-1}+8f_{i+1}-f_{i+2}\right]+\mathcal{O}(\Delta x^4) \]

5.9 Finites difference: forward vs centered

5.10 Second order derivatives with finite differences

By using (only) \(f_{i+1}\) and \(f_{i-1}\), we can get ,

Second order derivatives

\[ f''(x_i)=\frac{f_{i+1} - 2f_i + f_{i-1}}{\Delta x^2}+\mathcal{O}(\Delta x^2) \qquad(2)\]

We can also get this by using 1st-order central differences \[ f''(x_i)=\frac{d}{dx}\left(\frac{df}{dx}\right)=\lim_{\Delta x\to 0}\frac{f'_{i+1/2}-f'_{i-1/2}}{\Delta x}\approx \frac{(f_{i+1}-f_i)/\Delta x-(f_{i}-f_{i-1})/\Delta x}{\Delta x} \]

It is possible to systematically obtain higher order derivatives with arbitrary accuracy. See a list in http://en.wikipedia.org/wiki/Finite_difference_coefficients

5.11 Finite differences in two variables

For 2D domains (e.g., in PDEs), we can generalize the previous procedure by using, e.g, a five-point stencil \[ \left(\frac{\partial^2 f}{\partial x \partial y}\right)_{i,j}=\frac{\left(\frac{\partial f}{\partial y}\right)_{i+1,j}-\left(\frac{\partial f}{\partial y}\right)_{i-1,j}}{2\Delta x} \approx \frac{f_{i+1,j+1}-f_{i+1,j-1}-f_{i-1,j+1}+f_{i-1,j-1}}{4\Delta x\Delta y} \]

5.12 Your turn:

Task 1

1. Select a differential method and numerically solve the differential equation \[ \frac{dx}{dt} = f(t) = cos(t) \qquad(3)\] with \(\Delta t = 0.01\), and initial condition x(t=0) = 0, and interval \(t = [0,6\)\(]\).

2. Compare various methods with analytical solution of Equation 3

3. Analyze impact of the time step \(\Delta t\) on the quality of the solution

4. How much run-time differs when changing the size of the time step?

   Tip: For timing the computation use “%%timeit” command at the beginning of the notebook cells.

Task 2

Analytically prove Equation 2, using the formula for the first derivative

6 Partial differential equations

6.1 Characteristics of PDEs

• Curves in the space of independent variables of a PDE, along which the PDE has only total differentials.

  By integrating along those curves one can find the solutions of such a PDE.

• Example: \[ a\frac{\partial u}{\partial t} + b\frac{\partial u}{\partial x} = c \] with solutions \(u=u(x(t),t)\).

• Total derivative of \(u\): \[ \frac{du}{dt}=\frac{\partial u}{\partial x}\frac{dx}{dt} + \frac{\partial u}{\partial t} \]

• The characteristics are the solutions of: \[ \frac{dx}{dt}=\frac{b}{a},\;\; \frac{du}{dt}=\frac{c}{a} \] which means, for \(u(x_0,0)=u_0\), \[ x = x_0 + \frac{b}{a}t\\ u = u_0 + \frac{c}{a}t\\ \]

6.2 Characteristics of PDEs

• Curves in the space of independent variables of a PDE, along which the PDE has only total differentials
• Solution/characteristics: \[ x = x_0 + \frac{b}{a}t\\ u = u_0 + \frac{c}{a}t\\ \]
• Note that if \(c=0\), \(u\) is constant on the characteristics.
• For non constant \(c\), \(u\) grows linearly in time along a characteristic \(\Rightarrow\) the initial profile \(u(x,0)\) already determines the solution along the characteristics

6.3 Characteristics of PDEs

• For the second order PDE: \[ A\frac{\partial^2 f}{\partial x^2} + B\frac{\partial^2 f}{\partial x \partial y} + C\frac{\partial^2 f}{\partial y^2} + G = 0 \]

The characteristics are the solutions of: \[ \frac{dy}{dx} = \frac{B}{2A}\pm \frac{1}{2A}\sqrt{B^2 - 4AC} \]

This equation is:

• elliptic: \(B^2 - 4AC < 0\). No real solutions.
• parabolic if \(B^2 - 4AC = 0\). One real solutions.
• hyperbolic if \(B^2 - 4AC > 0\). Two real solutions.

6.4 Systems \(n\) of 1st order PDEs

• It is possible to write such a system of \(n\) PDEs in two dependent variables \(x\) and \(y\) as: \[ \left(\bar{A}^T\frac{dy}{dx} - \bar{B}^T \right)\cdot \vec{L} = 0 \] where \(\vec{L}\) is a vector of dimension \(n\), \(\bar{A}^T\), \(\bar{B}^T\), are \(n\times n\) matrices.

  The determinant leads to a \(n\)th order equation for \(dy/dx\), which can be classified as:

  • elliptic for no real roots
  • parabolic for \(1\) to \(n-1\) real roots and no complex roots
  • hyperbolic for \(n\) real roots

6.5 Classification

6.6 Linear Second Order PDEs in 2 independent variables

• The most general form of a linear second order PDE with independent variables \(x,y\) and dependent variable \(f(x,y)\) is: \[ A\frac{\partial^2 f}{\partial x^2} + B\frac{\partial^2 f}{\partial x \partial y} + C\frac{\partial^2 f}{\partial y^2} + D\frac{\partial f}{\partial x} + E\frac{\partial f}{\partial y} + G = 0 \]

  \(A,\dots, G\) are constants.

This equation is:

• elliptic if \(B^2 - 4AC < 0\)

• parabolic if \(B^2 - 4AC = 0\)

• hyperbolic if \(B^2 - 4AC > 0\)

6.7 Linear Second Order PDEs in 2 independent variables

\[ A\frac{\partial^2 f}{\partial x^2} + B\frac{\partial^2 f}{\partial x \partial y} + C\frac{\partial^2 f}{\partial y^2} + D\frac{\partial f}{\partial x} + E\frac{\partial f}{\partial y} + G = 0 \]

• Elliptic (\(B^2 - 4AC < 0\)): Steady state problems (gravity, electrostatic). Laplace equation (2D): \[ \nabla ^2\phi = \frac{\partial^2 \phi}{\partial x^2} + \frac{\partial^2 \phi}{\partial y^2} =0 \]
• Parabolic (\(B^2 - 4AC = 0\)): Time-dependent fluid equations where conduction/diffusion are important. Heat conduction equation \[ \frac{\partial T}{\partial t}=k\frac{\partial^2 T}{\partial x^2} \]
• Hyperbolic ( \(B^2 - 4AC > 0\)): Time-dependent fluid eqs where convection is important. Vibration, propagation. Wave equation \[ \frac{\partial^2U}{\partial t^2}=c^2\frac{\partial^2 U}{\partial x^2} \]

6.8 Multiple independent variables

• Let us consider a single dependent variable \(f\) but \(m\) independent variables: \[ \sum_{j=1}^{m} \sum_{k=1}^{m} a_{jk} \frac{\partial^2 f}{\partial x_j \partial x_k} + H=0 \] The character of this eq is determined by the eigenvalues (EV) of the matrix \(A=(a_{jk})\).
• This equation is
  • elliptic if all the EVs \(\neq0\) and all the same sign.
  • parabolic if any EV is 0
  • hyperbolic if all the EVs \(\neq0\) and all but one have the same sign.

6.9 Your turn

• Which kind of differential equation is:
  • Lorentz force
  • Schrödinger equation
  • Maxwell equations

6.10 Elliptic PDEs

• Typical for steady state problems (time-independent limit of hyperbolic or parabolic eqs), or for potentials (Poisson’s equation).
• Example: 2D Poisson equation (Or Laplace’s equation if the RHS=0) \[ \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} = \rho(x,y) \qquad(4)\]
• Solution for 4 boundary conditions: \[ u(x,y) = \sin(\pi x)\exp(-\pi y) \]
• No real characteristics
• Both the maximum and minimum are assumed on the boundary \(\partial R\).
• Any disturbance on the interior influences the entire computational domain.

6.11 Parabolic PDEs

• Example: heat conduction, magnetic flux diffusion, evolution of a velocity distribution function in velocity space due to particle collisions. \[ \frac{\partial u}{\partial t}= \sigma \frac{\partial^2 u}{\partial x^2} \qquad(5)\]
• Solution for one initial condition (typically Dirichlet) and two boundary conditions: \[ u(x,t) = \sin(\pi x)\exp(-\pi^2 t) \]
• Only one characteristic
• Any local perturbation influences the entire domain, but with decreasing magnitude for increasing distance from the source of the perturbation.

6.12 Hyperbolic PDEs

• Example: wave equation: \[ \frac{\partial^2 u}{\partial t^2} - c^2\frac{\partial^2 u}{\partial x^2} = 0 \qquad(6)\]
• Solution: Two sine waves, traveling into the positive (\(+x\)) and negative (\(-x\)) directions.
• Typical property: propagation of information through a system via waves (e.g. sound waves) along the characteristics of the equation.
• Examples: various time dependent forms of fluid equations without dissipation: heat conduction, viscosity, resistivity
• The number of initial conditions (this case: dependent variables and their derivatives) has to correspond with the number of unknowns or waves.
• Any point source along \(x\) will influence its vicinity along the corresponding characteristics (two straight lines for Eq. Equation 6.

6.13 Hyperbolic PDEs

Characteristics of a hyperbolic wave equation. Left: Domain of dependence. Right: Domain of influence (Jardin (2010)).

6.15 Boundary conditions

Let us define a domain, its boundary \(\partial R\), the coordinates \(n\) (outward normal to the boundary) and \(s\) (along the boundary), and the functions \(f\) and \(g\) on the boundary. The three types of BCs are:

1. Dirichlet conditions: \(u=f\) on \(\partial R\)
2. Neumann conditions: \(\frac{\partial u}{\partial n} =f\) or \(\frac{\partial u}{\partial s} =g\) on \(\partial R\)
3. Mixed (Robin) conditions: \(\frac{\partial u}{\partial n} + ku=f, \;\; k>0\) on \(\partial R\).

6.16 Initial vs Boundary conditions

(a) Initial (Cauchy) vs (b) Boundary value problems (Press et al. (2007)). In (a) the arrows indicate time.

6.17 Stability

6.18 Von Neumann stability analysis

• Based on a finite Fourier transform to convert spatial differences into an expression with the \(\xi\), such as the independent solutions (eigenmodes) are: \[ u_j^n=\xi^n(k) e^{ikj\Delta x} \qquad(7)\] \(j\) is the space index, \(n\) is the time index, \(i\) the imaginary unit, \(k\) the wavenumber, \(\Delta x\) the grid cell spacing.
• The difference equation will be unstable if \(|\xi(k)|>1\) for some \(k\). This means exponential growing modes.
• von Neumann analysis is local.

6.19 Von Neumann stability analysis 2

• Example: numerically unstable hyperbolic equation: wave equation written as a system of 1st order equations (each one is sometimes called “advective equations”). \[ \frac{\partial u}{\partial t} = -v \frac{\partial u}{\partial x} \] The (forward time centered space) FTCS representation is: \[ \frac{u_j^{n+1} - u_j^{n}}{\Delta t} = -v \left(\frac{u_{j+1}^{n} - u_{j-1}^{n}}{2\Delta x} \right) \qquad(8)\]
• Replacing Equation 7 in Equation 8 we get: \[ \xi(k)=1-i\frac{v\Delta t}{\Delta x}\sin(k\Delta x) \] which is \(>1\) for all \(k\), so that the scheme is unstable.

6.20 Von Neumann stability analysis 3

• Example: numerically stable hyperbolic equation: let us replace the time derivative in Equation 8 by its average: \[ u_j^{n+1} = \frac{1}{2}\left(u_{j+1}^{n} + u_{j-1}^{n}\right) -\frac{v\Delta t}{2\Delta x} \left(u_{j+1}^{n} - u_{j-1}^{n} \right) \qquad(9)\] (this is known as the Lax method), which gives the amplification factor: \[ \xi(k)=\cos(k\Delta x)-i\frac{v\Delta t}{\Delta x}\sin(k\Delta x) \] and it leads to the stability condition:

Courant-Friedrichs-Lewy

Courant-Friedrichs-Lewy stability condition (CFL, Courant condition): \[ \frac{v\Delta t}{\Delta x}\leq 1 \] Example: 2nd order hyperbolic wave equation (derived from the determinant of a matrix).

6.21 CFL condition

CFL condition (Jardin (2010)).

6.22 Von Neumann stability analysis 4

• It can be proven that Equation 9 is equivalent to a FTCS (forward time-centered space) representation of \[ \frac{\partial u}{\partial t} = -v \frac{\partial u}{\partial x} + \frac{(\Delta x)^2}{2\Delta t}\nabla^2 u \] where the extra term is a diffusion/dissipative term. Then, the Lax scheme is said to have numerical dissipation or viscosity.

6.23 Methods

6.24 Solution methods for parabolic equations 1

• Diffusion equation in 1D \[ \frac{\partial u}{\partial t}= D \frac{\partial^2 u }{\partial x^2} \qquad(10)\] where \(D\) is the diffuction parameter.

  This can also be written as a flux-conservative equation, with F=-\(D \frac{\partial u}{\partial x}\).

• \(D\geq 1\) in order for Eq. Equation 10 to be stable.

• Applying a FTCS scheme we get: \[ \frac{u_j^{n+1} - u_j^{n}}{\Delta t} = D \left(\frac{u_{j+1}^{n} -2u_j^n + u_{j-1}^{n}}{(\Delta x)^2} \right) \qquad(11)\]

6.25 Solution methods for parabolic equations 2

• The amplification factor is: \[ \xi(k)=1 -\frac{4D \Delta t}{(\Delta x)^2}\sin(k\Delta x/2) \] leading to the stability criterion: \[ \frac{2D\Delta t}{\Delta x}\leq 1 \qquad(12)\] the maximum allowed timestep is approx. the diffusion time across a cell of width \(\Delta x\).
• In general, the diffusion time across a spatial scale of size \(\lambda\) is \(\tau\sim \frac{\lambda^2}{D}\), which leads to a number of timesteps (based on Equation 12) of order \(\frac{\lambda^2}{\Delta x^2}\) \(\rightarrow\) too many.

6.26 Solution methods for parabolic equations 3

• One solution to change the small CFL timestep of the diffusion eq. is the Crank-Nicholson scheme, second-order accurate in time.
• Average of explicit and implicit FTCS schemes: \[ \frac{u_j^{n+1} - u_j^{n}}{\Delta t} = \frac{D}{2}\left(\frac{(u_{j+1}^{n+1} -2u_j^{n+1} + u_{j-1}^{n+1}) +(u_{j+1}^{n} -2u_j^{n} + u_{j-1}^{n}) }{(\Delta x)^2} \right) \qquad(13)\]
• The amplification factor predicts that the method is stable for any \(\Delta t\)
• Crank-Nicholson also works well for Schrödinger-like eqs. \[ i\frac{\partial \psi}{\partial t} = -\frac{\partial ^2\psi}{\partial x^2} + V(x)\psi \]

6.27 Solution methods for elliptic PDEs

• 2D Poisson equation \[ \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} = R(x,y) \qquad(14)\] with Dirichlet boundary conditions: \[ u(0,y), u(L,y) \qquad 0<y<L \] and \[ u(x,0), u(x,L) \qquad 0<x<L \]

The 2nd order finite difference equation for \(j,k\)-element (derived by the Taylor expansion) is: \[ u_{j+1,k} + u_{j-1,k} + u_{j,k+1} + u_{j,k-1} - 4u_{j,k} = \Delta x^2 R_{j,k}. \qquad(15)\] Let the \((N+1)^2\) vector of unknowns be denoted by \(\vec{x}\), which contains \(u_{j,k}\) elements.

The matrix form of Equation 15 is \(\hat{A}\cdot\vec{X}=\vec{b}\)

6.28 Solution methods for elliptic PDEs

Sparse matrix \(\hat{A}\):

• Must be inverted to solve for the unknowns (although the actual inverse is not normally calculated). Some methods are:
  • Direct inversion using Gauss or LU decomposition
  • Direct inversion using block tri-diagonal methods
  • Iterative methods (multigrid, Krylov space methods). A physical approach involves adding time derivatives terms (converting it to parabolic/hyperbolic equations), and advancing it in time until a steady state is obtained.
  • Direct methods based on transform techniques (e.g., Fourier transformation)

6.29 Multigrid methods

Schematics of multigrid method

• Iterative methods benefit greatly from a good initial guess: A coarse grid solution can be used to initialize a finer grid iteration
• On a given grid, the short wavelengths components of the error decay much faster than the long wavelengths. So, the coarse mesh is more efficient in handling long components of the error but the fine mesh is needed for the short wavelength components.
• The error satisfies the same matrix eqs. as the unknown.
• Two basic steps: nested iteration (for the initial guess), coarse grid correction.

6.30 Fast Fourier methods

• Fourier (spectral) methods can be applied to solve boundary problems (e.g., elliptic equations) that lead to large sparse linear systems of the form \(\hat{A}\cdot\vec{u}=\vec{b}\)
• The algorithm is based on a Fast Fourier transform (FFT), a clever trick to avoid the calculation of \(N^2\) multiplications and \(N(N-1)\) that are required for a direct finite Fourier transform. Instead, if \(N\) is equal to a power of 2, it leads to a product of \(log_2(N)\) sparse matrices.
• It works best when the PDE has constant coefficients and the boundaries coincide with the coordinate lines.

7 References

Benáček, Jan, Andrey Timokhin, Patricio A. Muñoz, Axel Jessner, Tatiana Rievajová, Martin Pohl, and Jörg Büchner. 2024. “Poynting flux transport channels formed in polar cap regions of neutron star magnetospheres.” arXiv e-Prints, May, arXiv:2405.20866. https://doi.org/10.48550/arXiv.2405.20866.

Boyd, T. J. M., and J. J. Sanderson. 2003. The Physics of Plasmas.

Büchner, Jörg. 2023. Space and Astrophysical Plasma Simulation. Methods, Algorithms, and Applications. https://doi.org/10.1007/978-3-031-11870-8.

Büchner, Jörg, C Dum, and M Scholer. 2003. Space Plasma Simulation. Berlin Heidelberg: Springer-Verlag.

Jardin, S. 2010. Computational Methods in Plasma Physics. CRC Press.

Lapenta, Giovanni. 2012. “Particle simulations of space weather.” J. Comput. Phys. 231 (3): 795–821. https://doi.org/10.1016/j.jcp.2011.03.035.

Matsumoto, H, and Y Omura. 1993. Computer Space Plasma Physics : Simulation Techniques and Software. Terra Scientific Publishing Company. https://www.terrapub.co.jp/e-library/cspp/.

Press, William H., Saul A. Teukolsky, William T. Vetterling, and Brian P. Flannery. 2007. Numerical Recipes: The Art of Scientific Computing. Third Edition. Cambridge University Press. https://doi.org/10.1137/1031025.

Umeda, Takayuki. 2012. “Simulation of Collisionless Plasma With the Vlasov Method.” In Int. Conf. Simul. Technol., edited by Brian S. Doherty and Amy N. Molloy, 315–32. Nova Science Publishers. http://www.jsst.jp/e/JSST2012/extended_abstract/pdf/12.pdf.

Winske, D., and N. Omidi. 1996. “A nonspecialist’s guide to kinetic simulations of space plasmas.” J. Geophys. Res. 101 (A8): 17287. https://doi.org/10.1029/96JA00982.