DifferentialEquations.jl: split differential equation with explicit propagators

Specific Domains Numerics
1

I have a PDE that is (not exactly but close to) du/dt + v \cdot du/dx = A(x)u, where x is 2D, u is a vector with two components and A is a 2x2 matrix. The initial data possibly has shocks (discontinuities). Right now I use finite differences + DifferentialEquations, but the numerical diffusion is killing me. I thought about using those fancy hyperbolic schemes, but they look complicated. Alternatively, I can use a splitting method (I know how to integrate without v, and I know how to integrate without A). Can I fit this into DifferentialEquations.jl to take care of the higher-order/timestepping for me? I did not find it in the doc, I just ask in case it’s hidden somewhere.

2

This is something we’re still working on. Operator splitting methods are in a very rudimentary state:

http://docs.juliadiffeq.org/latest/solvers/split_ode_solve.html

It’s one of the next big sets of algorithms we are working on though. Until we have those kinds of methods, I would recommend just using the SSP methods

http://docs.juliadiffeq.org/latest/solvers/ode_solve.html#Explicit-Strong-Stability-Preserving-Runge-Kutta-Methods-for-Hyperbolic-PDEs-(Conservation-Laws)-1

with a sufficient limiter. The other types of methods which would do well are the IMEX schemes which we don’t have yet:

https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/issues/191

So yeah, probably not satisfactory but it’s what we got. You can also just try throwing it into Sundials’ CVODE_BDF which might do well enough.

3

Another thing that helps with numerical differentiation is to use a higher order upwinding scheme. You can check out DiffEqOperators.jl for some of the operators for that:

The documentation there is still lacking though. But if you really need to get rid of numerical diffusion, you may need to use a better spatial discretization.

4

Actually I thought I knew how to solve without the A on a grid, but of course that’s bullshit: if v*delta_t is not a multiple of delta_x there’s no obvious way to transport the solution. So as you say I probably need to rethink the discretization. Still since in my case v is constant, maybe I can impose that v*delta_t = delta_x and solve exactly. I wonder if that works (and if I’m rediscovering Lattice Boltzmann or something like that). In any case I definitely need to read some literature. Thanks for the input!