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?

  • $\begingroup$ You don't. After the symbolic part, you use another package, scipy.integrate, for example. $\endgroup$ – nicoguaro Jan 4 '18 at 14:52
  • $\begingroup$ I'm not quite sure scipy.integrate is able to handle a SimPy object... $\endgroup$ – Matt Jan 4 '18 at 16:20
  • $\begingroup$ See @Wrzlprmft answer on how to do that. I have used ``lambdify``` and it works. $\endgroup$ – nicoguaro Jan 4 '18 at 16:24
  • $\begingroup$ OK, cheers mate. $\endgroup$ – Matt Jan 4 '18 at 16:37

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.

  • $\begingroup$ Thank you very much for your detailed answer! I hope that lambdify will be enough to tackle this. I accept your answer. $\endgroup$ – Matt Jan 4 '18 at 16:40

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.

  • $\begingroup$ Wah, cheers mate, the library seems impressive, I'm gonna try that too. For the time being, I'm not sure I will spend too much time on optimization, since I mainly have to compute each curve once. $\endgroup$ – Matt Jan 4 '18 at 18:45
  • $\begingroup$ Not my library, but spreading the word because it is pretty awesome. $\endgroup$ – Chris Rackauckas Jan 4 '18 at 19:16
  • $\begingroup$ @Matt I'm the author of pyodesys, let me know on github if you run into any trouble. The current design is not optimized for discretized PDEs (though a smaller number of variables should still work of course). You could also try to extend it by subclassing (I'd be happy to answer questions too if you go down that road). For performance you can either generate C++ or have SymEngine jit-compile expressions. $\endgroup$ – Bjoern Dahlgren Jan 9 '18 at 22:31
  • $\begingroup$ @Bjoern Dahlgren Thank you very much for help offered! I therefore will be contacting you if I have some troubles running what I need to. :) $\endgroup$ – Matt Jan 11 '18 at 11:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.