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

Specific Domains Numerics
41

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

https://benchmarks.sciml.ai/html/MOLPDE/burgers_spectral_wpd.html

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
end

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.

6 Likes
42

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.

8 Likes
43

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

44

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!!

2 Likes
45

Yeah, why not? It increases the band size of the matrix, but it can be helpful in reducing the number of terms required. You do need to be careful with numerical errors when using higher order operators, but 4th order should be fine. I think the 4th order central difference should be used more, it’s just difficult to, by hand, work out the boundary conditions since you no longer have a single term to the left and right. But yes, when a library takes care of that for the user, why not?

Rodas5 won’t scale as well to tons of ODEs, but it is L-stable while other methods like LSODA and CVODE_BDF are not, so Rodas5 will have a much better time with highly oscillatory solves. RodauIIA5 and KenCarp4 should have similar properties. The newer KenCarp47 is worth a try too.

You’d have to analytically look at the problem a bit, but properties like dissipation will occur numerically. The error has to go somewhere! You can mix in the ManifoldProjection callback to correct for this though. Or, if you have a mode that you know should be analytically zero, you can just have a callback that projects that back to zero after each step (this retains the order of the solution BTW so it’s perfectly safe!)

Yes, DiffEqOperators (now with fixed higher dimensional BC handling) allows for a non-uniform grid. You could use that to do an adaptive grid, and that could be something fun to investigate.

2 Likes
46

Hi Chris, thanks so much for your detailed reply. Been trying everything out.

The higher orders library for finite difference is extremely useful; is there some
theoretical resource or documentation explaining how they actually calculate
d^2/dx^2 (spatial 2nd derivative) for the n>2 orders?

Thanks for the commentary on the L-stability, Rodas5 seems like the way to go.
KenCarp47 sounds good to try too, if I get the chance.

Still working on figuring out how to use ManifoldProjection (now Manifolds?) and
DiffEqOperators for those purposes (energy cancelling and non-uniform grids).

Thanks again!

47

It’s just Fornberg’s algorithm: https://web.njit.edu/~jiang/math712/fornberg.pdf

5 Likes
48

Thanks, that reference is perfect!

49

Actually, going through the paper, they also give one-sided finite difference formulas too; but implementing them by hand my results are strange.

In analogy with CenteredDifference(...) (and BandedMatrix(CenteredDifference(...))), does Julia have commands like BackwardDifference(...) and ForwardDifference(...) (or LeftDifference(...) and RightDifference(...) ) ?

Thanks!

50

The UpwindOperators are forward and reverse differencing.

51

Thanks again Chris, I found your web article on the topic.