0

I'm working on a system of differential equation that I solve with scipy.integrate.solve_ivp. For each variation term of this system I combine different intermediate variables. I'm looking for a way to store those intermediate variables and not only the resulting variables.

I've thinking to use a list or some other kind of iterable to store the intermediate values directly in the func passed to solve_ivp but it doesn't seem to be the best approach and doesn't allow me to use the interpolants, when using a dense_output. I could also recompute those intermediate variables by using the results of solve_ivp, but as the computation are somehow complicated it should be way easier to store directly the values and only compute them only once.

My func function looks like that :

def dY(t, Y):
    C1, C2 = Y

    dC1 = A(C1, C2) + B(C1, C2)
    dC2 = C(C1, C2) + B(C1, C2)

    return [dC1, dC2]

I would like to store the values of A, B, C and D at each step, as are the values of C1 and C2 (for those to whom it speaks, it's a chemical engineering problem and A, B, C and D are different sink and source terms)

Thanks for any helps !

Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51
Rowin
  • 435
  • 2
  • 9

1 Answers1

2

Use the values of the solutions and call the functions after the fact. Do not use intermediate values that get stored inside the dY function. Stiff ODE integrators will call the function multiple times during a time step, you will end up with meaningless data getting stored. Nonstiff ODE integrators might also attempt multiple time steps within one try which will have the same problem.

Here is an example of how to reconstruct A and B after the fact:

from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

tspan = (0, 10)
y0 = [1, 2]

def A(C1, C2):
    return 2*C1 + 3*C2 - C1**2
def B(C1, C2):
    return 1*C1 + 4*C2
def C(C1, C2):
    return 1*C1 + 3*C2 - C2**2
def D(C1, C2):
    return 2*C1 + 5*C2
def dY(t, Y):
    C1, C2 = Y

    dC1 = A(C1, C2) + B(C1, C2)
    dC2 = C(C1, C2) + B(C1, C2)

    return [dC1, dC2]
sol = solve_ivp(dY, tspan, y0)

fig,ax = plt.subplots()
ax.plot(sol.t, sol.y[0], label='C1')
ax.plot(sol.t, sol.y[1], label='C2')
ax.legend()
fig2,ax2 = plt.subplots()
ax2.plot(sol.t, A(sol.y[0], sol.y[1]), label='A')
ax2.plot(sol.t, B(sol.y[0], sol.y[1]), label='B')
ax2.legend()
Matthew Flamm
  • 448
  • 3
  • 13
  • I assumed it's the better way to do it. I then refactored my code to make it easier to reconstruct A and B after the resolution. – Rowin Oct 04 '19 at 08:39