How to start Partial Differential Equation (PDE) Solving in Julia for modified Wave Equation

Here’s an example which uses ApproxFun to discretize space for Burger’s equation:

What it’s doing and why it’s written like that is described in the notes/video. Basically, pick a space, then take every derivative and turn it into an operator, and write out the operator ODE on some finite coefficient size version (to keep it simple). The whole Fourier discretization in this form is:

S = Fourier()
n = 512
x = points(S, n)
D  = (Derivative(S) → S)[1:n,1:n]
T = ApproxFun.plan_transform(S, n)
Ti = ApproxFun.plan_itransform(S, n)

û₀ = T*cos.(cos.(x.-0.1))
tmp = similar(û₀)
p = (D,T,Ti,tmp,similar(tmp))
function burgers_nl(dû,û,p,t)
    D,T,Ti,u,tmp = p
    mul!(tmp, D, û)
    mul!(u, Ti, tmp)
    mul!(tmp, Ti, û)
    @. tmp = tmp*u
    mul!(u, T, tmp)
    @. dû = - u

Here this is defining the ODE in phase space (with a hat). The initial condition is the initial condition transformed into phase space. The RHS of the ODE is:

  1. Apply a derivative to u (in phase space)
  2. Transform the current function back to point space
  3. Transform du/dx back to point space
  4. Do u * du/dx in point space
  5. Transform that result back to phase space (and -)

This then gives du/dt = - u * du/dx written as an ODE in phase space that transforms back to do the nolinear portion.

You can take this idea and apply it directly to your function. Just choose a function space appropriate for R, or transform R -> [-1,1] using something like y = (x / (1-x^2), in which case Chebyshev might be a good choice, where you do exactly this procedure to specify the PDE except change the definition of S and change which operators you compute.


Let’s keep it civil here, folks. PDEs are a tough topic, and multiple users in this thread have considerable domain expertise in different PDE subfields. Best to cool down, give OP time to get up to speed on the topic, and try what’s already been suggested.


Thanks Chris!
I will read & try it out.
Still have to watch your video from before.

Hi again, after doing a ridiculous amount of reading (and watching Chris’ extremely helpful video), most of the reply posts on this topic are starting to make sense to me.

It is now clear that the best way to address my equation is with the Method of Lines for the resulting ODE diff eq’s in time (t), after first discretizing in space (x). The best way here to make x discrete is to use what you’re calling Finite Differences with DiffEqOperators. A possible 2nd choice for backup/verification is discretizing x with ApproxFun, but for my G(t,x) function that should only work using the Hermite basis.
Also, I must use “totally absorbing” boundary conditions, which add a couple of “ghost cells” on either side (in x) that keep the spatial derivatives zero at the edges, to kill off any edge reflections or aliased waves from outside the calc domain.

I spent a lot of time reading about Finite Volume methods, but those are bad for hyperbolic eq’s unless one uses something like Lax-Wendroff. But that causes unacceptable phase shifts here. (Solvable with flux limiters, but then the method gets too complicated and mysterious.) And Finite Element methods (and its spinoffs) aren’t called for in this kind of problem.

Within a day or two, I will post a couple of rough draft code snippets outlining my attempted approach; feedback/corrections would be most welcome, if anyone has time. For now, I just have a couple of general questions about this method:

(i) With DiffEqOperators, getting uxx with CenteredDifference(2, order, Δx, N), the usual setting is order=2. But why not crank it up to order=3,4,5…? Any reason not to do to that (e.g., accuracy side effects), other than computational burden?

(ii) My G(t,x) varies so much in magnitude that I may have to use a LARGE number of points in discrete x, ending up with LOTS of ODE’s. (100? 1000? 5000??) For my ODE work in June, the best solvers I found were Vern9 (non-stiff), RadauIIA5 and Rodas5 (stiff). Chris’ talk also mentions Rodas5 as good for the PDE problem. (KenCarp4 gave me less accuracy on ODE’s, and CVODE_BDF has other warnings.) But the ODE solvers documentation says for large numbers of ODE’s, to use LSODA (or VCABM). Is Rodas5 good enough, or are LSODA or VCABM necessary?

(iii) The global methods with ApproxFun mix things together non-locally (by design). For my problem, there will be regions with zero wave energy, that should numerically stay zero. Will ApproxFun (with Hermite) be fairly bad at that, putting real wave energy where physically there shoudn’t be any? (Also, would ApproxFun mess up my attempted totally-absorbing boundary conditions?)

(iv) Unlikely, but is it possible to make Δx of the spatial grid variable, e.g., finer spaced where the variations are stronger? (Even more unlikely, can the program sense this and distribute points appropriately; and most unlikely, maybe even adapt evolve the Δx grid over time?)

Lots of info & questions, I know, but thanks for any feedback!!