Developping PDE with Python symbolically and numericaly

Question

I feel like publishing some previous works from my PhD thesis. I was using Mathematica to build a system of 2N partial differential equations for 2N functions by symbolic spatial Taylor expansion, then numerically integrate them with NDSolve with respect to time.

Mathematica had a strange behaviour for some N orders, stopping to some apparent numerical singularities I was not able to manage at that time. Furthermore, I have an access to a supercomputer, which doesn't offer Mathematica. I would therefore like to rebuild my model in Python.

For the symbolic part, I guess SymPy will handle easily the Taylor expansion part (notice that after order 2 or 3, it took far more than a full page to Mathematica to write the PDE system, should I ask it to). But how do you numerically (Runge-Kutta or whatever) solve the objects coming out of SymPy?

Answer 1

But how do you numerically (Runge-Kutta or whatever) solve the objects coming out of SymPy?

That depends a lot on your specific problem and what level of optimization you need.

• Most ODE solvers (such as those included in SciPy) require that you provide them a Python function representing the right-hand side of the ODE. As you do not want things to be horribly slow, you have to take care of vectorisation somehow:

  • Depending on how your PDE is structured, this may work using SymPy’s lambdify with a NumPy backend.
  • An alternative is SymPy’s ufuncify function.
  • Finally, I wrote a module that takes SymPy input, compiles it, and feeds it into a SciPy integrator. This does not require vectorisable structures, but it can also not exploit them very well (as it is made with non-vectorisable problems in mind).

  Either way, his means that you have to take care of spacial discretisations on the SymPy level. You also will not be able to make use of PDE-specific optimisation techniques, in particular those involving GPUs; you may get multi-kernel support though (if that’s how you want to make use of the supercomputer)

• You convert SymPy’s results into whatever input is needed by PDE-specific Python tools. This probably involves meta-programming (code that writes code) of some sort, but could be rather harmless. Note that SymPy features several code printers and similar that may assist you with this.

Answer 2

I want to point out that PyODESys is a Python library that does exactly what you're looking for: send SymPy expressions to ODE solvers. It's nice because it links to more solvers, but @Wrzlprmft's JITCODE is probably better because it compiles the resulting expression and the runtime of the user's function is very important for solving an ODE efficiently.