I wish to implement a function in python to numerically solve what we might call an intermediate value problem. By an intermediate value problem I mean the union of a final value problem (fvp) and an initial value problem (ivp) in the following sense. Let $f:[a,b] \times \Omega \to \mathbb{R}$, with $\Omega \subseteq \mathbb{R}^n$, let $c \in (a,b)$ and let $y_0 \in \Omega$.

Consider the following fvp in the function $y:[a,c] \to \Omega$: $(\mathrm{fvp}): dy/dt = f(t,y) \ \mathrm{for} \ t \in (a,c), \ y(c) = y_0$ and let the solution be denoted $y_\mathrm{fvp}$.

Similarly, consider the following ivp in the function $y:[c,b] \to \Omega$: $(\mathrm{ivp}): dy/dt = f(t,y) \ \mathrm{for} \ t \in (c,b), \ y(c) = y_0$ and let the solution be denoted $y_\mathrm{ivp}$.

Then define the solution $y:[a,b] \to \Omega$ of the intermediate value problem by $y(t) = y_\mathrm{ivp} \ \mathrm{if} \ t \geq c, \ y(t) = y_\mathrm{fvp} \ \mathrm{if} \ t < c.$

To solve this problem numerically, we can convert the fvp to an ivp by a simple transformation, as discussed e.g. here.

I have done so in the following example code, which seems to work.

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

a = -2
b = 1
c = 3

def f(t, y): return 0.7 * y  
def f_left(s,y): return -f(b-s,y)

right = solve_ivp(f, [b, c], [1], dense_output = True)
left = solve_ivp(f_left, [0, b-a], [1], dense_output = True)

def my_sol(t):
    if t >= b:
        return right.sol(t)
        return left.sol(b-t)

sol_vect = np.vectorize(my_sol,signature='()->(n)')

t_right = right.t
y_right = right.y[0]

t_left = b-left.t[-1::-1]
y_left = left.y[0]
y_left = y_left[-1::-1]

t = np.linspace(a,c)
y = sol_vect(t)

plt.plot(t_right,y_right,'ro',label = "right y")
plt.plot(t_left,y_left,'bo',label = "left y")
plt.plot(t,y,'g',label = "sol_vect")

As you can see, I have defined the solution of the intermediate value problem using numpy.vectorize. However, my concern is that this is somehow inefficient. I am also concerned that there could be some sort of issue with array shapes that I fail to handle in a more general case.

A cleaner solution would be to simply concatenate the dense outputs from the fvp and the ivp. The issue is that the dense outputs for the fvp would have to be reversed according to $y(t) \rightarrow y(b-t)$. Is there a clean and efficient way to do this? Is there overall a better approach to solving the problem?

  • $\begingroup$ Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking. $\endgroup$
    – Community Bot
    Sep 9 at 17:14
  • $\begingroup$ There is no need for a special ODE function, solve_ivp is fully capable to integrate backwards, left = solve_ivp(f, [b,a], [1], dense_output = True). You could also implement sol_vect directly by splitting the input array at b and returning the joined result arrays. $\endgroup$ Sep 10 at 8:46
  • 1
    $\begingroup$ @LutzLehmann Thank you very much for pointing this out! Ok, that seems to work. In case anyone else is reading this, here's how you can get the joined solution as an instance of OdeSolution: joined = OdeSolution(t,interpolants), where t = np.concatenate((left.t[-1:0:-1],right.t)) and where interpolants = left.sol.interpolants[-1::-1] + right.sol.interpolants (list concatenation). $\endgroup$ Sep 10 at 15:32


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.