Solving PDEs in Julia

JuliaCon 2018 workshop

Chris Rackauckas

How are finite elements, multigrid methods, ODE solvers, etc. all the same topic?

Teach a man to fish: we won't be going over pre-built domain specific PDE solvers, instead we will be going over the tools which are used to build PDE solvers.

While the basics of numerically solving PDEs is usually taught in mathematics courses, the way it is described is not suitable for high performance scientific computing. Instead, we will describe the field in a very practical "I want to compute things fast and accurately" style.

What is a PDE?

A partial differential equation (PDE) is a differential equation which has partial derivatives. Let's unpack that.

A differential equation describes a value (function) by how it changes. u' = f(u,p,t) gives you a solution u(t) by you inputing / describing its derivative. Scientific equations are encoded in differential equations since experiments uncovers laws about what happens when entities change.

A partial differential equation describes a value by how it changes in multiple directions: how it changes in the x vs y physical directions, and how it changes in time.

Thus spatial properties, like the heat in a 3D room at a given time, u(x,y,z,t) are described by physical equations which are PDEs. Values like how a drug is distributed throughout the liver, or how stress propogates through an airplane hull, are examples of phonomena described by PDEs.

You are interested in PDEs.

Part 1: Representations of PDEs as other mathematical problems: Learning by example: the Poisson Equation

This will either be an overview where I will reframe the most common method of solution, or your first walkthrough of a PDE solver!

The best way to solve a PDE is...

By converting it into another problem!

Generally, PDEs are converted into:

  • Linear systems: Ax = b find x.

  • Nonlinear systems: G(x) = 0 find x.

  • ODEs: u' = f(u,p,t), find u.

There are others:

  • S(tochastic) DEs: du = f(u,p,t)dt + g(u,p,t)dW, find mean(u).

... Yes experts, there are more but we will stick to the usual stuff.

Let's introduce some shorthand: ux=uxu_x = \frac{\partial u}{\partial x}. The Poisson Equation is the PDE:

for x[0,1]x \in [0,1]. In one dimension:

b(x)b(x) is some known constant function (known as "the data"). Given the data (the second derivative), find uu.

How do we solve this PDE?

First Choice: Computational Representation of a Function

First we have to choose how to computationally represent our continuous function u(x)u(x) by finitely many numbers. Let's start with the most basic way (and we'll revisit the others later!). Let Δx\Delta x be some constant discretization size and let xi=iΔxx_i = i\Delta x. Then we represent our function u(x)ui=u(xi)u(x) \approx {u_i} = {u(x_i)}, i.e. we represent a continuous function by values on a grid (and we can assume some interpolation, like linear interpolation)

Δx = 0.1
x = 0:Δx:1
u(x) = sin(2π*x)
using Plots
Julia DiffEq (Julia)

Second Choice: Discretization of Derivatives

Forward difference: uxδ+u=ui+1uiΔxu_x \approx \delta^+ u = \frac{u_{i+1} - u_{i}}{\Delta x}

Backward difference: uxδu=uiui1Δxu_x \approx \delta^- u = \frac{u_{i} - u_{i-1}}{\Delta x}

Central difference for 2nd derivative: uxxδux=δ+uxδuxΔxu_{xx} \approx \delta u_{x} = \frac{\delta^+ u_x - \delta^- u_x}{\Delta x}

This gives the well-known central difference formula:

Quick Recap

We just made two choices:

  • Represent our function by an evenly spaced grid of points

  • Represent the derivative by the central difference formula

Given these two choices, how can we re-write our PDE?

The Representation of Our PDE

Remember we want to solve Δu=b \Delta u = bon x[0,1] x \in [0,1]

  • u(x)u(x)is now a vector of points ui=u(xi)u_i = u(x_i)

  • b(x)b(x)is now a vector of points bi=b(xi)b_i = b(x_i)$b_i = b(x_i)$

  • The second derivative is now the function 2ux(xi)=ui+12ui+ui1Δx2\frac{\partial^2 u}{\partial x}(x_i) = \frac{u_{i+1} - 2u_i + u_{i-1}}{\Delta x^2}

Thus we have a system of ii equations:

But Wait...

What happens at i=0i=0?

Translating back to ui=u(xi)=u(iΔx)u_i = u(x_i) = u(i\Delta x):

The last point is out of the domain! In order to solve this problem we have to impose boundary conditions. For example, let's add the following condition to our problem: u(0)=u(1)=0u(0) = u(1) = 0. Then the 00th point is determined: u0=0u_0 = 0, and the 1st point is:

so we can solve it!

The Linear Representation of Our Derivative

Notice that if U=[u1 u2  uN1 uN], U=\left[\begin{array}{c} u_{1}\ u_{2}\ \vdots\ u_{N-1}\ u_{N} \end{array}\right], then

This is a linear equation!

We know AA, BB is given to us, find UU.

Let's walk through a concrete version with some Julia code now. Let's solve:

Δx = 0.1
x = Δx:Δx:1-Δx # Solve only for the interior: the endpoints are known to be zero!
N = length(x)
B = sin.(2π*x)
A = zeros(N,N)
for i in 1:N, j in 1:N
  abs(i-j)<=1 && (A[i,j]+=1)
  i==j && (A[i,j]-=3)
A = A/(Δx^2)
Julia DiffEq (Julia)
9×9 Array{Float64,2}: -200.0 100.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 100.0 -200.0 100.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 100.0 -200.0 100.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 100.0 -200.0 100.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 100.0 -200.0 100.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 100.0 -200.0 100.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 100.0 -200.0 100.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 100.0 -200.0 100.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 100.0 -200.0
# Now we want to solve AU=B, so we use backslash:
U = A\B
Julia DiffEq (Julia)

Did we do that correctly?

This equation is simple enough we can check via the analytical solution.

Integrate it twice:

# Now we want to solve AU=B, so we use backslash:
plot!([0;x;1],-sin.(2π*[0;x;1])/4(π^2),label="Analytical Solution")
Julia DiffEq (Julia)


We solved the PDE Δu=b\Delta u = b by transforming our functions into vectors of numbers UU and BB, transforming second derivative into a linear operator (a matrix) AA, and solving AU=BAU = B using backslash.

Does this method generally apply? Pretty much.

Because derivatives are linear, when you discretize a function, the derivative operators become linear operators = matrices

Semilinear Poisson Equation

Now the right hand side is dependent on uu! Let's choose the same discretization and the same representation of the derivatives. Then once again Δu=AU\Delta u = AU for the same matrix AA. Now we get the equation:

Find the vector of UU which satisfy this nonlinear system! If we redefine:

then we are looking for the vector UU which causes G(U)=0G(U) = 0.

tl;dr: Semilinear equations convert into nonlinear rootfinding problems

Semilinear Heat Equation

Discretize the function the same way as before. This once again makes uxx=AUu_{xx} = AU. Thus letting UtU_t be the time derivative of each point in the vector, we get:

but since there's now only one coordinate, let the derivative be by time. Then we can write this as:

This is an ODE!

tl;dr: Time-dependent PDEs (can) convert into ODE problems!

Loose End: Higher Dimensions Do The Same Thing

Now let's look at Δu=uxx+uyy=b(x,y) \Delta u = u_{xx} + u_{yy} = b(x,y) on x[0,1],y[0,1]x \in [0,1], y \in [0,1].

In this case, you can let ui,j=u(xi,yi)u_{i,j} = u(x_i,y_i) where xi=iΔxx_i = i\Delta x and yi=iΔyy_i = i\Delta y. You can list out all of the ui,ju_{i,j} into a vector UU by lexicographic ordering:

for some A. Now solve AU=BAU=B

Part 1 Summary

To solve a PDE,

  • You choose a way to represent functions.

  • You choose a way to represent your derivative (and on the function representation, your derivative representation is a matrix!)

Then when you write out your PDE, you get one of the following problems:

  • Linear systems: Ax = b find x.

  • Nonlinear systems: G(x) = 0 find x.

  • ODEs: u' = f(u,p,t), find u.

Part 2: The Many Ways to Discretize

There are thus 4 types of packages in the PDE solver pipeline:

  • Packages with ways to represent functions as vectors of numbers and their derivatives as matrices

  • Packages which solve linear systems

  • Packages which solve nonlinear rootfinding problems

  • Packages which solve ODEs

In this part we will look at the many ways you can discretize a PDE.

Part 2.1: Packages to Represent Functions and Derivatives

There are four main ways to represet functions and derivatives as vectors:

  • Finite difference method (FDM): functions are represented on a grid. Packages: DiffEqOperators.jl (still developing)

  • Finite volume method (FVM): functions are represented by a discretization of its integral. Currently no strong generic package support.

  • Finite element method (FEM): functions are represented by a local basis. FEniCS.jl, JuliaFEM and JuAFEM.jl

  • Spectral methods: functions are represented by a global basis. FFTW.jl and ApproxFun.jl

Finite Difference method: DiffEqOperators.jl

DiffEqOperators.jl is part of the JuliaDiffEq ecosystem. It automatically develops lazy operators for finite difference discretizations (functions represented on a grid). For example, to represent uxxu_{xx}, we'd do:

using Pkg
pkg"add DiffEqOperators BandedMatrices"
Julia DiffEq (Julia)
using DiffEqOperators
# Second order approximation to the second derivative
order = 2
deriv = 2
Δx = 0.1
N = 9
A = CenteredDifference(deriv, order, Δx, N)
Julia DiffEq (Julia)
DerivativeOperator{Float64,1,false,Float64,SArray{Tuple{3},Float64,1,3},SArray{Tuple{0},SArray{Tuple{4},Float64,1,4},1,0},Nothing,Nothing}(2, 2, 0.1, 9, 3, [100.0, -200.0, 100.0], 4, 0, StaticArrays.SArray{Tuple{4},Float64,1,4}[], StaticArrays.SArray{Tuple{4},Float64,1,4}[], nothing, nothing)

This A is lazy: it acts A*u like it was a matrix but without ever building the matrix by overloading * and directly computing the coefficients. This makes it efficient, using O(1)\mathcal{O}(1) memory while not having the overhead of sparse matrices!

This package also makes it easy to generate the matrices without much work. For example, let's get a 2nd order discretization of uxxxxu_{xxxx}:

using BandedMatrices
BandedMatrix(CenteredDifference(4, 2, Δx, N))
Julia DiffEq (Julia)
9×11 BandedMatrix{Float64,Array{Float64,2},OneTo{Int64}}: 20000.0 -90000.0 160000.0 -140000.0 … ⋅ ⋅ ⋅ 10000.0 -40000.0 60000.0 -40000.0 ⋅ ⋅ ⋅ 0.0 10000.0 -40000.0 60000.0 0.0 ⋅ ⋅ 0.0 0.0 10000.0 -40000.0 0.0 0.0 ⋅ 0.0 0.0 0.0 10000.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 … 10000.0 0.0 0.0 0.0 0.0 0.0 0.0 -40000.0 10000.0 0.0 ⋅ 0.0 0.0 0.0 60000.0 -40000.0 10000.0 ⋅ ⋅ 0.0 0.0 160000.0 -90000.0 20000.0

Brief brief brief overview of finite element methods

Represent your function u(x)=iciφi(x)u(x) = \sum_i c_i \varphi_i(x) with some chosen basis φi(x)\varphi_i(x).

"Matrix Assembly" = calculate the matrix representations of the derivatives in this function representation. The core of an FEM package is its matrix assembly tools.

Finite difference method is good if your domain is a square / hypercube. Finite element methods can solve PDEs on more complicated domains.

FEM Package 1: FEniCS.jl

FEniCS is a well-known finite element package for Python (fellow NumFOCUS project!). It lets you describe the kind of PDE you want to solve and what elements (basis functions) you want to discretize with in a DSL, and it handles the rest.

FEniCS.jl is a wrapper over FEniCS which is maintained by the JuliaDiffEq organization.

  • Pro: very full featured (since it's wrapping an existing package). Linear solvers are built in.

  • Con: not Julia-based, so it's missing a lot of the fancy Julia features (generic programming, arbitrary number types, etc.). Python loading is tricky.

FEM Package 2: JuliaFEM

JuliaFEM is an organization with a suite of packages for performing finite element discretizations. It focuses on FEM discretizations of physical PDEs and integrates with Julia's linear solver, nonlinear rootfinding, and DifferentialEquations.jl libraries to ease the full PDE solving process.

FEM Package 3: JuAFEM.jl

JuAFEM.jl is a FEM toolbox. It gives you functionality that makes it easier to write matrix assembly routines. The below is the heat equation example.

Spectral Methods

Like finite element methods, spectral methods represent a function in a basis:u(x)=iciφi(x)u(x) = \sum_i c_i \varphi_i(x) with some chosen basis φi(x)\varphi_i(x)

"Spectral" usually refers to global basis functions. For example, the Fourier basis of sines and cosines.

Spectral Packages 1: FFTW.jl

Okay, this is not necessarily a spectral discretization package. However, it is a package to change from a pointwise representation of a function to a Fourier representation via a Fast Fourier Transform (FFT). Thus if you want to find the coefficients ckic_ki for

you'd use:

using Pkg
pkg"add FFTW"
Julia DiffEq (Julia)
using FFTW
x = range(0,2π,length=100)
u(x) = sin(x)
freqs = fft(u.(x))[1:length(x)÷2 + 1]
c = 2*abs.(freqs/length(x))'
Julia DiffEq (Julia)
1×51 Adjoint{Float64,Array{Float64,1}}: 5.27356e-18 0.994805 0.0135642 0.00760374 … 0.000635192 0.000634878

This corresponds to basically saying that sin(x)sin(x) is represented in the Fourier basis as:

It almost got it, and you can see the slight discretization error. This goes away as you add more points. But now you can represent any periodic function!

Spectral Packages 2: ApproxFun.jl

ApproxFun.jl is a package for easily approximating functions and their derivatives in a given basis, making it an ideal package for building spectral discretizations of PDEs. It utilizes lazy representations of infinite matrices to be efficient and save memory. It uses a type system to make the same code easy to translate between different basis choices.

Representing a Function and its Derivative in the Fourier Basis

Let's represent

in the Fourier basis, and build a matrix representation of the second derivative in this basis:

using Pkg
pkg"add ApproxFun"
Julia DiffEq (Julia)
using ApproxFun
S = Fourier()
n = 100
T = ApproxFun.plan_transform(S, n)
Ti = ApproxFun.plan_itransform(S, n)
x = points(S, n)
r = (T*cos.(cos.(x.-0.1)))'
Julia DiffEq (Julia)
1×100 Adjoint{Float64,Array{Float64,1}}: 0.765198 2.3331e-18 -2.36654e-17 -0.0456556 … -1.77238e-17 -1.9984e-17
D2 = Derivative(S,2);
Julia DiffEq (Julia)

Now let's change to the Chebyshev basis:

This basis is kckTk(x) \sum_k c_k T_k(x) for TkT_k the Chebyschev polynomials (1,x,2x21,4x33x,)(1, x, 2x^2-1, 4x^3 - 3x, \ldots)

S = Chebyshev()
n = 100
T = ApproxFun.plan_transform(S, n)
Ti = ApproxFun.plan_itransform(S, n)
x = points(S, n)
r = (T*cos.(cos.(x.-0.1)))'
Julia DiffEq (Julia)
1×100 Adjoint{Float64,Array{Float64,1}}: 0.712955 -0.0529031 0.155596 0.01012 … -1.22125e-17 -6.66134e-18
D2 = Derivative(S,2);
Julia DiffEq (Julia)

Part 2 Summary

Using these packages, you can easily translate your PDE functions to coefficient vectors and your derivatives to matrices:

  • Use spectral methods or finite difference methods for cases with "simple enough" boundary conditions and on simple (square) domains

  • Use finite element packages to discretize complex spatial domains

There are a lot more factors. Sometimes given discretizations are better/worse on given PDEs. There are mathematicians who spend their entire life researching the differences between discretization methods!

But now...

Our PDE functions are now vectors of coefficients. Our PDE derivatives are now matrices. We are left with one of the following to solve:

  • Linear systems: Ax = b find x.

  • Nonlinear systems: G(x) = 0 find x.

  • ODEs: u' = f(u,p,t), find u.

Each of these problems is a specific area of research in and of itself!

Part 3: Equation Solving

In order to do equation solving correctly, we must make note of one important fact:

The matrices which come out of a derivative discretization are often very sparse!

  • The finite difference matrix [1 -2 1] only has 3 non-zero values per row no matter what N is! This type of matrix is called tridiagonal.

  • For matrices which have more general semi-diagonal structures which are only non-zero near the center, they are called banded matrices.

  • Derivatives in Fourier space are diagonal!

  • Generic sparse matrices arise from FEM discretizations (since the derivatives only depend on nearby basis elements which are very few!).

Thus solving PDEs is more specifically equation solving with large numbers of equations and high sparsity with structure

Part 3.1: Solving Large Sparse Linear Systems

When solving AU=BAU=B where AA arises from the derivative term in the PDE discretization, the structure of AA determines what type of method and what pacakage to use.

  • If A is tridiagonal, then Julia's \ will use a specialized fast method.

  • If A is small enough, then \ is a multithreaded LU or QR factorization which is fine.

  • If A is a banded matrix (many spectral discertizations and finite difference methods), then uses BandedMatrix(A)\B from BandedMatrices.jl utilizes a fast LU/QR for this matrix type.

  • If A is a block banded matrix (blocking arises from lexicographical ordering when 2+ dimensional), then BlockBandedMatrix(A)\B from BlockBandedMatrices.jl utilizes a fast LU/QR for this matrix type.

  • If A is a small enough sparse matrix, then Julia's Base \ will use SuperLUMT, a multithreaded sparse LU which will be fast if it has enough memory (hence small enough).

  • If A is a large sparse matrix (or matrix-free), then iterative methods are required. For this, there is IterativeSolvers.jl and JuliaSmoothOptimizers/Krylov.jl

Iterative Solvers

Iterative solvers repeatedly calculate A*U \ k in order to find a sequence such that UkUwhereAU=BU_k \rightarrow U \\ where \\ AU=B.

Because iterative solvers only use multiplication, "matrix-free" operators, i.e. Julia types which just define *, can be used in these methods!

IterativeSolvers.jl has many methods which are specialized for different forms of sparsity structures. For example, cg is a fast method for symmetric positive definite matrices. The fallback method for a general A is gmres. Let's give it a try on some random AU=B problem.

using Pkg
pkg"add IterativeSolvers"
Julia DiffEq (Julia)
n = 10
A = rand(n, n)
B = rand(n)
# Let's use the gmres method from IterativeSolvers.jl
using IterativeSolvers, LinearAlgebra
# Same thing as A\B
U = gmres(A, B, tol = 1e-8) # Get at least within 1e-8 of the solution!
norm(A*U - B)
Julia DiffEq (Julia)

Using a Matrix-Free Type

Let's create a type that is the second order central difference derivative [1 -2 1] and solve the Poisson equation Δu=b\Delta u = b with our matrix-free operator by using gmres instead of \. We could just use DiffEqOperators.jl's lazy matrix here, but let's build the full example from scratch:

struct SizedStrangMatrix
Base.eltype(A::SizedStrangMatrix) = Float64
Base.size(A::SizedStrangMatrix) = A.size
Base.size(A::SizedStrangMatrix,i::Int) = A.size[i]
A = SizedStrangMatrix((length(B),length(B)))
import LinearAlgebra
function LinearAlgebra.mul!(C,A::SizedStrangMatrix,B)
    for i in 2:length(B)-1
        C[i] = B[i-1] - 2B[i] + B[i+1]
    C[1] = -2B[1] + B[2]
    C[end] = B[end-1] - 2B[end]
Base.:*(A::SizedStrangMatrix,B::AbstractVector) = (C = similar(B); mul!(C,A,B))
using IterativeSolvers, LinearAlgebra
U = gmres(A,B,tol=1e-14)
norm(A*U - B)
Julia DiffEq (Julia)

LinearMaps.jl is a nice package for helping build matrix-free discretizations!


To solve the linear system AU=BAU=B faster, you can use a preconditioners to partial solve the linear problem and thereby present a simpler problem to the iterative method. This is a deep deep area of research and the best preconditioner for your problem is highly dependent on your PDE, so I'm just going to point to some of the best packages in the area:

You take one of these, plop it in as an additional argument to gmres and then it can go faster!

Linear Solver Parallelism Libraries

Because of the ubiquity of solving large sparse linear systems in PDEs, there exists libraries which are dedicated to parallelizing the large scale sparse linear algebra. Julia packages for this are:

These packages allow you to distibute a linear solve computation amongst a whole cluster mixed with GPUs! But fundamentally the techniques are the same.

Part 3.1 Recap

Linear solving requires specializing on the matrix type:

  • If small enough or the structure is known, special methods should be used

  • If large enough or no sparsity structure to explot, Krylov methods like gmres need to be used.

  • If iterative (Krylov) methods are used, then a good preconditioner can heavily speed up the solving.

  • \, Julia's LinearAlgebra special matrix types, BandedMatrices.jl, BlockBandedMatrices.jl, and IterativeSolvers.jl are your friends!

Part 3.2 Nonlinear Solvers

Recall that if your data is a nonlinear function of your unknown then your PDE discretization doesn't produce a linear system but instead a nonlinear system G(x)=0G(x) = 0. There are a few packages for this case:

  • NLSolve.jl

  • Sundials.jl (KINSOL)

  • MINPACK.jl

These have not been benchmarked against each other ( go do it!) but have different characteristics.

But due to features, KINSOL.jl is the recommended one for now for PDEs, with NLsolve features coming soon

Quick Note: Don't use optimization routines to solve rootfinding problems

G(x)=0G(x) = 0 is theoretically the same as finding the minimum of G(x)\Vert G(x) \Vert (which should be zero!). However, this is bad:

  • Optimization methods are generally rootfinding methods on the derivative, since an optima is found when the derivative is zero!

  • Optimization methods generally use higher order derivatives. For example, Newton's method for optimization uses the second derivative (the Hessian), instead of the first derivative (the Jacobian), so it's more error prone and more expensive.

Just don't do it.

Newton's Method

where J(zn)J(z_n) is the Jacobian of G(zn)G(z_n). To solve this without needed inverses, you can rewrite this with dz=zn+1zndz = z_{n+1} - z_n and see that

Thus every step of Newton's method is actually a linear solving problem, where the matrix is the Jacobian of GG. Since GG is from our discretized PDE and all of the nonlinear behavior is local, the Jacobian of GG has the same structure and sparsity pattern as the AA from before! This means that:

Newton's method is essentially nonlinear solving via repeat linear solving. Almost all of the time is spent linear solving

All of the same tricks from before apply: we want to specialize on banded matrices, we want to use Krylov methods, and preconditioners, etc.

Part 3.2 Recap

To solve nonlinear systems G(x)=0G(x) = 0, you either do:

  • Some simple iteration scheme (Picard iteration, Powell iterations, Anderson acceleration). These have a smaller area of convergence and converge slower, but do not require solving linear systems.

  • Some form of Newton's method. This requires solving linear systems.

Usually a Newton method is required. In this case, almost all of the cost comes from solving the linear systems. Right now, KINSOL in Sundials.jl is the best tool because of the flexibility it gives in the linear solver choices, though NLsolve.jl will be adding similar features shortly (PRs are in progress).

Part 3.3 ODE Solvers

Recall that if you leave a derivative un-discretized you get an ODE system. Example:

You could discretize this time derivative the same way as you did space, getting a two-dimensional system and solving a big linear/nonlinear equation. But there can be advantages to leaving it as an ODE:

  • The full linear/nonlinear system is much larger. If you have N points in space and M time points, then you have u as a size NM vector! This grows fast!

  • Newton's method requires a good starting point. For a huge system in time and space, initial values for u(xi,tj)u(x_i,t_j) can be hard to come by, and Newton's method can diverge.

Keeping one part undiscretized and solving the ODE can alieviate these problems.

Features of a PDE-derived ODE

A PDE-derived ODE tends to have the following features:

  • Stiffness. This is due to the CFL condition which does require ΔtΔxk<1\frac{\Delta t}{\Delta x^k} < 1 for stability

  • Large sparse Jacobian. This is because of the matrices from the derivative terms!

What ODE methods are applicable to this case?

DifferentialEquations.jl has a wealth of unique options specifically for reducing the cost of integration on these types of problems.

Quick Note: Order of Integrator

FAQ: I did a second order discretization in space, does it make sense to use a higher than second order integrator in time?

Answer: Yes! The order of the integrator isn't just for reducing error. Higher order integrators can be more efficient! You can always incrase the error by taking larger time steps (increasing the tolernace) and thereby use less steps with a more efficient method.

Implicit and Rosenbrock ODE Integrators

These methods step by repeatedly solving nonlinear systems. Remember that solving nonlinear systems boils down to repeatedly solving linear systems. Thus for the ODE u=f(u,p,t)u'=f(u,p,t) with Jacobian JJ, an implicit ODE solver spends most of its time solving:

Since this JJ has the same characteristics as before, solving this efficiently requires specializing the linear solver to the sparsity pattern (bandedness, etc.) or using a good preconditioner in an iterative solver.

Good choices in DifferentialEquations.jl:


  • KenCarp4

  • Rodas5

Example: Spectral time stepping for Heat Equation in Fourier Space

using Pkg
pkg"add ApproxFun Sundials Plots DifferentialEquations"
Julia DiffEq (Julia)
using ApproxFun, Sundials, Plots; gr()
S = Fourier()
n = 100
x = points(S, n)
T = ApproxFun.plan_transform(S, n)
Ti = ApproxFun.plan_itransform(S, n)
# Convert the initial condition to Fourier space
u₀ = T*cos.(cos.(x.-0.1)) 
D2 = Derivative(S,2)
L = D2[1:n,1:n]
using LinearAlgebra
# The equation is trivial in Fourier space
heat(du,u,L,t) = mul!(du, L, u) 
prob = ODEProblem(heat, u₀, (0.0,10.0),L)
# Specialize the linear solver on the diagonalness of the Jacobian
sol = solve(prob, CVODE_BDF(linear_solver=:Diagonal); reltol=1e-8,abstol=1e-8);
Julia DiffEq (Julia)
# The solution is in Fourier space, so use inverse to transform back
Julia DiffEq (Julia)

Pseudospectral Solving the Reaction-Diffusion Equation

Take the PDE

For some nonlinear uu. In the Fourier basis, uxx=Auu_{xx} = Au for a diagonal matrix uu, so this is a nice way to solve the Heat Equation as above (by solving in Fourier space, and converting the solution back). We can solve a nonlinear equation by, whenever we need to apply the nonlinear ff, transform uu back to the origional coordinates and applying ff on uu before transforming back. This pseudospectral discretization is seen as:

gives an ODE for how Fu\mathcal{F}u evolves, and from that we can recover the true solution via the inverse Fourier transform F1(Fu)\mathcal{F}^{-1}(\mathcal{F}u)

Implicit-Explicit (IMEX) Integrators

But wait, we shouldn't solve this as one whole ODE. If we split the ODE into two parts:

then we notice that the first part is the linear part and has the stiff term, while the second part is nonlinear and (can be) non-stiff. If we only are implicit on the first part, we can use a linear solver instead of a nonlinear solver to significantly reduce the amount of work! Integrators which allow you to split and do two separate methods at the same time are called IMEX integrators.

In DifferentialEquations.jl, IMEX integrators to be aware of are SplitODEProblem

  • CNAB2

  • SBDF2

  • KenCarp3

  • KenCarp4

# By making it Diagonal, DiffEq internally specializes the linear solver
using DiffEqOperators, LinearAlgebra, DifferentialEquations
A = DiffEqArrayOperator(Diagonal(L)) 
function f(dû,,tmp,t)
  # Transform u back to point-space
  # apply nonlinear function 0.75sqrt(u)-u in point-space
  @. tmp = 0.75sqrt(tmp) - tmp 
  mul!(dû,T,tmp) # Transform back to Fourier space
# Define u' = Au + f
prob = SplitODEProblem(A, f, u₀, (0.0,10.0), similar(u₀));
Julia DiffEq (Julia)
# has_Wfact() issue
sol = solve(prob, KenCarp4())
Julia DiffEq (Julia)

Exponential Integrators

A more recent set of integrators are the exponential integrators. These methods do not require solving a linear system, thus avoiding the most costly calculation of the implicit methods. However, they must perform a Krylov version of matrix exponential times vector multiplications, w=exp(γA)vw = exp(\gamma A)v. Recent literature has been showing large performance gains of exponential integrators over traditional implicit methods. DifferentialEquations.jl is the first open source library in a high performance language to include exponential integrators, and a brand new (adaptive) Krylov-based approach for large sparse systems was just released as part of GSoC 2018. These methods will be continued to be improved but are ready for widespread use and benchmarking is to come.

sol = solve(prob, ETDRK4(), dt=0.1)
Julia DiffEq (Julia)


Additional methods in DifferentialEquations.jl should be explored. The Runge-Kutta Chebyshev (RKC) methods are stabilized explicit methods which can solve stiff equations without linear solving. Options allow for swapping out the internal Newton method to Picard and Anderson iteration to not require linear solving. Passing a linsolve allows you to take control over the linear solver technique: make it an AMG-preconditioned GMRES. Please see the documentation for a full feature list.

SSP Methods

Strong-Stability Presurving Runge-Kutta methods are "more stable" explicit methods for PDEs. Some PDEs, like hyperbolic PDEs, are not amenable to implicit solvers and thus require explicit solvers. SSPRK methods can increase the allowable stepsize to increase the efficiency!

Part 3.3 Recap

PDEs can also be solved by leaving an axis undiscretized and using an ODE solver. This can increase efficiency, stability, and reduce memory consumption. Classes of methods to utilize are:

  • Implicit and Rosenbrock(-W) Methods

  • IMEX Methods

  • Exponential Integrators

  • Runge-Kutta Chebyshev Methods

  • SSPRK Methods

Why Julia?

  • IterativeSolvers.jl is compatible with matrix-free types via * overloads able to be used in the nonlinear and ODE solvers

  • BandedMatrices.jl and BlockBandedMatrices.jl are packages for specialization of linear solvers on common PDE matrix types

  • KINSOL from Sundials.jl, and soon NLsolve.jl.

  • DifferentialEquations.jl

  • Julia is a language which allows for zero-cost abstractions, letting one compile code specific to the application and compose packages without extra overhead

The native Julia methods of DifferentialEquations.jl also let you swap in all of the Julia tools for linear solvers and allows type-genericity (which allows the use of GPUs for example), allowing you to utilize the full arsenal of tools with these unique implementations.

Conclusion: Solving PDEs Takes Layers of Tools

You need discretization tooling, linear solvers, nonlinear solvers, and finally ODE solvers to build an efficient PDE solver.

We showed how you can use simple loops to write simple PDE solvers, but the efficient methods require using packages in order to get the latest and most efficient methods.

Many of the latest and most efficient methods only have implementations in Julia, and the JuliaDiffEq organization is committed to continuing the development of such toolling.

Near Future

In the near future we with to provide a set of problem types and algorithms for common PDEs. Example: Reaction-Diffusion Equation

# Define the PDE problem to solve
prob = ReactionDiffusionProblem(domain,discretization,f)
# Solve it with a pseudospectral IMEX method
# Or an EPIRK method

We are now very close.


I am deeply indebted to every JuliaDiffEq contributor. I would like to especially those who have specifically been involved in developing and funding the large range of tools which have been demonstrated in this talk:

  • Yingbo Ma (@YingboMa)

  • Shivin Srivastava (@shivin9)

  • Shubham Maddhashiya (@sipah00)

  • Xingjian Guo (@MSeeker1340)

  • Yiannis Simillides (@ysimillides)

  • Sheehan Olver (@dlfivefifty)

  • Patrick Kofod Mogensen (@pkofod)

  • Harmen Stoppels (@haampie)

  • @dextorious

  • Bart Janssens (@barche)

  • David Widmann (@devmotion)

  • Hendrik Ranocha (@ranocha)

  • Alexey Stukalov (@alyst)

  • Spencer Lyon (@sglyon)

  • Jiahao Chen (@jiahao)

  • Viral Shah (@ViralBShah)

  • Vijay Ivaturi (@vjd)

  • Jesse Perla (@jlperla)

  • Kristoffer Carlsson (@KristofferC)

  • Antoine Levitt (@antoine-levitt)

  • Jukka Aho (@ahojukka5)

Runtimes (1)