CODEWITHSUNDEEP
pythonscipyodeodeint
Eiem Lupus
Eiem Lupus

Reputation: 23

How to go through exact points with solve_ivp?

I have a system of ODE, and I want to change the value of the variables when the solver reaches an exact point in time. What I'm trying to do is similar to this example in Julia

What I've tried is to do is using else and if to check if the time t has been reached. However, this doesn't work because the array that contains the time doesn't pass through that exact point. Also tried to use t_eval without success.

Using the same problem showed in the julia example we have:

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

def f(t, u):
    if t == 4 and u[0] < 4:
        u[0] += 10
    du = np.empty(1)
    du[0] = -u[0]
    return du

u0=[10.0]
V = 1
t = [0,10]
u0 = [10.0]
sol = solve_ivp(f,t,u0)

plt.plot(sol.t, sol.y[0])

Is there a way to replicate the julia callback function?

Upvotes: 2

Views: 1404

Answers (1)

Warren Weckesser
Warren Weckesser

Reputation: 114811

Instead of trying to exactly replicate that code with a callback, it seems simpler to solve the equations in two stages:

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


def f(t, u):
    du = -u
    return du


dose = 10.0

u0 = dose
V = 4

t0 = 0
t1 = 4
t2 = 10

t = np.linspace(t0, t1, 100)
sol1 = solve_ivp(f, [t0, t1], [u0], dense_output=True, t_eval=t)

u0 = sol1.y[0, -1]
if u0 / V < 4:
    u0 += dose

t = np.linspace(t1, t2, 150)
sol2 = solve_ivp(f, [t1, t2], [u0], dense_output=True, t_eval=t)

plt.plot(sol1.t, sol1.y[0], 'b')
plt.plot([sol1.t[-1], sol2.t[0]], [sol1.y[0, -1], sol2.y[0, 0]], 'b--')
plt.plot(sol2.t, sol2.y[0], 'b')
plt.grid()
plt.xlabel('t')
plt.show()

plot

Upvotes: 1

Related Questions