Can we modify the solution vector between integrations steps with scipy.integrate.ode, using VODE?

user_12296951 0
Viewed 78

I am trying to get a solution for a stiff ODE problem where at each integration step, i have to modify the solution vector before continuing on the integration. For that, i am using scipy.integrate.ode, with the integrator VODE, in bdf mode. Here is a simplified version of the code i am using. The function is much more complex than that and involve the use of CANTERA.

from scipy.integrate import ode
import numpy as np
import matplotlib.pyplot as plt

def yprime(t,y):
    return y

vode = ode(yprime)
vode.set_integrator('vode', method='bdf', with_jacobian=True)

y0 = np.array([1.0])
vode.set_initial_value(y0, 0.0)
y_list = np.array([])
t_list = np.array([])
while vode.t<5.0 and vode.successful:
    vode.integrate(vode.t+1e-3,step=True)
    y_list = np.append(y_list,vode.y)
    t_list = np.append(t_list,vode.t)

plt.plot(t_list,y_list)

Output:
Plot generated

So far so good.

Now, the problem is that within each step, I would like to modify y after it has been integrated by VODE. Naturally, i want VODE to keep on integrating with the modified solution. This is what i have tried so far :

while vode.t<5.0 and vode.successful:
    vode.integrate(vode.t+1e-3,step=True)
    vode.y[0] += 1  # Will change the solution until vode.integrate is called again
    vode._y[0] += 1 # Same here.

I also have tried looking at vode._integrator, but it seems that everything is kept inside the fortran instance of the solver. For quick reference, here is the source code of scipy.integrate.ode, and here is the pyf interface scipy is using for VODE.

Has anyone tried something similar ? I could also change the solver and / or the wrapper i am using, but i would like to keep on using python for that.

Thank you very much !

ode scipy python
1 Answers

For those getting the same problem, the issue lies in the Fortran wrapper from Scipy.

My solution was to change the package used, from ode to solve_ivp. The difference is that solve_ivp is entirely made with Python, and you will be able to hack your way through the implementation. Note that the code will run slowly compared to the vode link that the other package used, even though the code is very well written and use numpy (basically, C level of performances whenever possible).

Here are the few steps you will have to follow.

First, to reproduce the already working code :

from scipy.integrate import _ivp # Not supposed to be used directly. Be careful.
import numpy as np
import matplotlib.pyplot as plt

def yprime(t,y):
    return y

y0 = np.array([1.0])
t0 = 0.0
t1 = 5.0

# WITHOUT IN-BETWEEN MODIFICATION
bdf = _ivp.BDF(yprime,t0,y0,t1)

y_list = np.array([])
t_list = np.array([])
while bdf.t<t1:
    bdf.step()
    y_list = np.append(y_list,bdf.y)
    t_list = np.append(t_list,bdf.t)

plt.plot(t_list,y_list)

Output :

y(t) without modifications between steps

Now, to implement a way to modify the values of y between integration steps.

# WITH IN-BETWEEN MODIFICATION
bdf = _ivp.BDF(yprime,t0,y0,t1)

y_list = np.array([])
t_list = np.array([])
while bdf.t<t1:
    bdf.step()
    bdf.D[0] -= 0.1 # The first column of the D matrix is the value of your vector y.
    # By modifying the first column, you modify the solution at this instant.
    y_list = np.append(y_list,bdf.y)
    t_list = np.append(t_list,bdf.t)

plt.plot(t_list,y_list)

Gives the plot :

y(t) with modification between steps

This does not have any physical sense for this problem, unfortunately, but it works for the moment.

Note : It is entirely possible that the solver become unstable. It has to do with the Jacobian not being updated at the right time, and so one would have to recalculate it again, which is performance heavy most of the time. The good solution to that would be to rewrite the class BDF to implement the modification before the Jacobian Matrix is updated. Source code here.

user_122969510
Related