Brian,
Can you take a look at what odeint returns? Specifically, at:
‘mused’ a vector of method indicators for each successful time step: 1:
adams (nonstiff), 2: bdf (stiff)
I suspect it's using Adams all the way, which means it's doesn't need a
Jacobian.
Emil
On 12/10/15 1:51 PM, Barry Smith wrote:
Brian,
I see two distinct issues here
1) being apply to apply your right hand side efficiently and
2) what type of ODE integrators, if any, can work well for your problem with
its funky, possibly discontinuous right hand side?
1) Looking at the simplicity of your data structure and function evaluation I think you
should just write your right hand side functions in C. The code would be no more
complicated than it is now, from what you showed me. Even with "just in time"
compilation likely you lose at least a factor of 10 by coding in python, maybe 100? I
don't know. It looks to me like an easy translation of your current routines to C.
2) This is tricky. You a) cannot compute derivatives? and b) the function and
derivatives may not be smooth?
If the function was well behaved you could use finite differences to
compute the Jacobian (reasonably efficiently, the cost is just some number of
function evaluations) but with a flaky right hand side function the finite
differences can produce garbage.
You could use explicit schemes with a small enough timestep to satisfy any
stability condition and forget about using implicit schemes. Then it becomes
crucial to have a very fast right hand side function (since you will need many
time steps). If you are lucky you can use something like a 4th order RK scheme
(but again maybe with a non smooth right hand side may you can't).
I am no expert. Perhaps Emil who is far more knowledgable about these
things has questions?
Barry
On Dec 10, 2015, at 12:52 PM, Brian Merchant <bhmerch...@gmail.com> wrote:
Hi Barry,
Here's some non-trivial example code:
https://gist.github.com/bmer/2af429f88b0b696648a8
I have still made some simplifications by removing some phase variables,
expanding on variable names in general, and so on.
The rhs function itself is defined on line 578. The functions referred to
within it should be all defined above, so you can have a peek at them as
necessary.
Starting from line 634 I show how I use the rhs function. In particular, note the "disjointed"
evaluation of the integral -- I don't just evaluate from 0 to t all at one go, but rather evaluate the
integral in pieces (let's call the time spent between the end of one integral evaluation, and the start of
the next integral evaluation a "pause"). This is so that if there were multiple amoebas, during the
"pause", I can take into account changes in some of the parameters due to contact between one
amoeba and another -- poor man's simplification.
Please let me know if this is what you were looking for. I wouldn't be
surprised if it wasn't, but instead would be happy to try to rework what I've
got so it's more in line with what would be meaningful to you.
Kind regards,
Brian
On Wed, Dec 9, 2015 at 2:18 PM, Barry Smith <bsm...@mcs.anl.gov> wrote:
I prefer the actual code, not the mathematics or the explanation
On Dec 9, 2015, at 3:42 PM, Brian Merchant <bhmerch...@gmail.com> wrote:
Hi Barry,
Could send an example of your "rhs" function; not a totally trivial example
Sure thing! Although, did you check out the exam I tried to build up in this
stackexchange question, along with a picture:
http://scicomp.stackexchange.com/questions/21501/is-it-worth-switching-to-timesteppers-provided-by-petsc-if-i-cant-write-down-a
I ask because that's probably the best I can do without using as little math as
possible.
Otherwise, what I'll do is take a couple of days to carefully look at my work,
and write up a non-trivial example of a difficult-to-differentiate RHS, that
still is a simplification of the whole mess -- expect a one or two page PDF?
Kind regards,
Brian
On Mon, Dec 7, 2015 at 12:45 PM, Barry Smith <bsm...@mcs.anl.gov> wrote:
Brian,
Could send an example of your "rhs" function; not a totally trivial example
Barry
On Dec 7, 2015, at 11:21 AM, Brian Merchant <bhmerch...@gmail.com> wrote:
Hi all,
I am considering using petsc4py instead of scipy.integrate.odeint (which is a
wrapper for Fortran solvers) for a problem involving the solution of a system
of ODEs. The problem has the potential to be stiff. Writing down its Jacobian
is very hard.
So far, I have been able to produce reasonable speed gains by writing the RHS functions
in "something like C" (using either numba or Cython). I'd like to get even more
performance out, hence my consideration of PETSc.
Due to the large number of equations involved, it is already tedious to think
about writing down a Jacobian. Even worse though, is that some of the functions
governing a particular interaction do not have neat analytical forms (let alone
whether or not their derivatives have neat analytical forms), so we might have
a mess of piecewise functions needed to approximate them if we were to go about
still trying to produce a Jacobian...
All the toy examples I see of PETSc time stepping problems have Jacobians
defined, so I wonder if I would even get a speed gain going from switching to
it, if perhaps one of the reasons why I have a high computational cost is due
to not being able to provide a Jacobian function?
I described the sort of problem I am working with in more detail in this
scicomp.stackexchange question, which is where most of this question is copied
from, except it also comes with a toy version of the problem I am dealing with:
http://scicomp.stackexchange.com/questions/21501/is-it-worth-switching-to-timesteppers-provided-by-petsc-if-i-cant-write-down-a
All your advice would be most helpful :)
Kind regards,Brian
