5

I'm using sympy to solve a simple linear system of equations.

It's a coupled ODE, there are time-derivatives of variables and I need to solve the system of equations for the highest derivatives. Since sympy doesn't allow me to solve for statements like phi_1.diff(t), I've replaced all derivatives with placeholder symbols.

For example:

phi.diff(t).diff(t) + phi(t) =0

becomes

ddphi + phi(t) = 0

This works fine. The solutions are correct and I can simulate the system - it's a pendulum: https://youtu.be/Gc_V2FussNk

The problem is that solving the system of equations (with linsolve) takes very long.

For just 2 equations, it takes 2 seconds. For 3 equations, it's still calculating (after over 10 minutes).

EDIT: @asmeurer advised me to try out solve instead. For n=3, linsolve took about 34 minutes - I only made one measurement. solve takes 31 seconds (averages over 3 runs).

Still, I believe that a linear 3x3 system should be solved in fractions of a second.

And for n=4 solve becomes unbearably slow, too (still calculating)

I've formatted the code and created an iPython notebook: http://nbviewer.jupyter.org/gist/lhk/bec52b222d1d8d28e0d1baf77d545ec5 If you scroll down a little, you can see the formatted output of the system of equations and directly below that the call to linsolve

The equations are rather long but strictly linear in the second derivatives. I'm sure that this system can be solved. All I need to do is solve a 3x3 system of linear equations where the coefficients might be symbols.

Is there a more performant way to do this ?

lhk
  • 27,458
  • 30
  • 122
  • 201
  • You should also provide the jacobian to **odeint** if you want the solution to be more performant. – mtadd Jun 21 '16 at 21:09
  • Thanks for the tip. But for the 2-pendulum odeint wasnt the bottleneck. So far I havent even reached that line of code with 3 pendulums. – lhk Jun 21 '16 at 21:11
  • 1
    Your differential equations don't look linear to me. Some of the first order derivative terms of phi are squared in your set of equations. – mtadd Jun 21 '16 at 21:17
  • @mtadd sure that's right. But I'm not solving for the first order derivatives, only for the second order. The solver should be able to handle them like constants. – lhk Jun 21 '16 at 21:19
  • 1
    Looks like `solve` is faster than `linsolve`. – asmeurer Jun 22 '16 at 04:38
  • @asmeurer Oh, you're right. And by quite a bit. I actually waited for linsolve to complete: took 35min. Solve on the other hand only needs 31 seconds. This is only a solution for n=3 though. I've updated the question accordingly – lhk Jun 24 '16 at 07:12
  • 1
    The problem is that `linsolve` calls `simplify`, which is slow, and there's no way to disable it. `solve` does have a flag, though, `simplify=False`. Does that help any? – asmeurer Jun 24 '16 at 16:51
  • @asmeurer thanks for that amazing tip. Just tried that flag and it makes a gigantic difference. 3 pendulum in 1.7s , 4 pendulum in 4.6s, 5 pendulum in 30s. This is a huge progress. Still I think that a 5x5 system of equations shouldn't take 30s to solve. Do you think there's still room for improvement ? Honestly, I've got basically no experience with sympy and don't know its limitations – lhk Jun 24 '16 at 20:00
  • 1
    You could play with the other flags to [solve](http://docs.sympy.org/latest/modules/solvers/solvers.html#sympy.solvers.solvers.solve). I'm not an expert in it. – asmeurer Jun 24 '16 at 20:27
  • @asmeurer I did that. There is a flag which disables the conversion of floats to rationals - this brings another huge boost. Now I can simulate pendulums up to order 7 and 8. The integration of the ODE is becoming the bottleneck. Thank you very much. If you convert this into an answer: "Tuning the flags will improve performance" I'll accept it. – lhk Jun 24 '16 at 20:53

2 Answers2

8

solve (not linsolve) has some flags that you can set which can make it faster:

  • simplify=False: disables simplification of the result.
  • rational=False: disables the automatic conversion of floats to rationals.

There is a warning in the solve docstring that rational=False can lead to some equations not being solvable due to issues in the polys, so be aware that that's a potential issue.

asmeurer
  • 86,894
  • 26
  • 169
  • 240
2

I have found that solve can be very slow in jupyter notebook if you have run sp.init_printing() before your equations. I have a module "equations" where I write my equations and solve them.

This is faster: 

import sympy as sp
import equations
sp.init_printing()

Than this:

import sympy as sp
sp.init_printing()
import equations
Martin Alexandersson
  • 1,269
  • 10
  • 12