1

I have the following problem:

I have to use an ode-solver to solve a chemical reaction equation. The rate constants are functions of time and can suddenly change (puls from electric discharge).

One way to solve this is to keep the stepsize very small hmax < dt. This results in a high comp. affort --> time consuming. My question is: Is there an efficient way to make this work? I thought about to def hmax(puls_ON) with plus_ON=True within the puls and plus_ON=False between. However, since dt is increasing in time, it may dose not even recognize the puls, because the time interval is growing hmax=hmax(t).

A time-grid would be the best option I thin, but I don't think this is possible with odeint?

Or is it possible to somehow force the solver to integrate at a specific point in time (e.g. t0 ->(hmax=False)->tpuls_1_start->(hmax=dt)->tpuls_1_end->(hmax=False)->puls_2_start.....)?

thx

Matthias K
  • 93
  • 9

1 Answers1

0

There is an optional parameter tcrit for the odeint that you could try:

Vector of critical points (e.g. singularities) where integration care should be taken.

I don't know what it actually does but it may help to not simply step over the pulse.

If that does not work you can of course manually split your integration into different intervals. Integrate until your tpuls_1_start. Then restart the integration using the results from the previous one as initial values.

astoeriko
  • 730
  • 8
  • 20
  • Thx for the suggestion with tcrit. However, odeint seems to be not the best method for this problem. I recognized on the scipy page they now suggest to use scipy.integrate.solve_ivp instead of odeint? Finally, I think that VODE should be the best option for such ODE problems, since other plasma chemistry solver also use it. integrate.ode features VODE, but is it also fast? (comp. to odenit). And do you/anyone knows some examples, how to use it. It seems to be more tricky (have to calc jacobi...)? My odeint needs a lot of args, since df(t) depends on many time dependent parameter! – Matthias K Apr 30 '19 at 14:48
  • If you need the Jacobian you might want to have a look at how to symbolically represent your ODE and compute the Jacobian using [sympy](https://www.sympy.org/en/index.html). The symbolic expression can then be converted to a python function. [Here](https://www.sympy.org/scipy-2017-codegen-tutorial/notebooks/32-chemical-kinetics-symbolic-construction.html) is an example (with odeint, but it should be similar with other solvers). – astoeriko May 02 '19 at 01:14