Specify different timepoints for boundary conditions of Scipy's "odeint"?

Question

I am trying to numerically solve a system of two nonlinear ODEs. I am using Scipy's odeint function. odeint requires an argument y0 which specify the initial conditions. However, it seems to assume that the initial conditions of y0 begin at the same timepoint (i.e both conditions are at t=0). In my case, I want to specify two different boundary conditions that are specified for different times (i.e. omega(t=0) = 0, theta(t=100) = 0). I can't seem to figure out how to do this and would greatly appreciate any help!

Some example code below:

from scipy.integrate import odeint

def pend(y, t, b, c):
    theta, omega = y
    dydt = [omega, -b*omega - c*np.sin(theta)]
    return dydt

b = 0.25
c = 5.0
t = np.linspace(0, 100, 101)

# I want to make these initial conditions specified at different times
y0 = [0, 0]

sol = odeint(pend, y0, t, args=(b, c))

Answer 1

odeint solves the initial value problem. The problem that you describe is a two-point boundary value problem. For that, you can use scipy.integrate.solve_bvp

You could also take a look at scikits.bvp1lg and scikits.bvp_solver, although it looks like bvp_solver hasn't been updated in a long time.

For example, here is how you could use scipy.integrate.solve_bvp. I changed the parameters so the solution does not decay so fast and has a lower frequency. With b = 0.25, the decay is rapid enough that θ(100) ≈ 0 for all solutions where ω(0) = 0 and |θ(0)| is on the order of 1.

The function bc will be passed the values of [θ(t), ω(t)] at t=0 and t=100. It must return two values that are the "residuals" of the boundary conditions. That just means it must compute values that must be 0. In your case, just return y0[1] (which is ω(0)) and y1[0] (which is θ(100)). (If the boundary condition at t=0 had been, say ω(0) = 1, the first element of the return value of bc would be y0[1] - 1.)

import numpy as np
from scipy.integrate import solve_bvp, odeint
import matplotlib.pyplot as plt


def pend(t, y, b, c):
    theta, omega = y
    dydt = [omega, -b*omega - c*np.sin(theta)]
    return dydt


def bc(y0, y1, b, c):
    # Values at t=0:
    theta0, omega0 = y0

    # Values at t=100:  
    theta1, omega1 = y1

    # These return values are what we want to be 0:
    return [omega0, theta1]


b = 0.02
c = 0.08

t = np.linspace(0, 100, 201)

# Use the solution to the initial value problem as the initial guess
# for the BVP solver. (This is probably not necessary!  Other, simpler
# guesses might also work.)
ystart = odeint(pend, [1, 0], t, args=(b, c,), tfirst=True)


result = solve_bvp(lambda t, y: pend(t, y, b=b, c=c),
                   lambda y0, y1: bc(y0, y1, b=b, c=c),
                   t, ystart.T)


plt.figure(figsize=(6.5, 3.5))
plt.plot(result.x, result.y[0], label=r'$\theta(t)$')
plt.plot(result.x, result.y[1], '--', label=r'$\omega(t)$')
plt.xlabel('t')
plt.grid()
plt.legend(framealpha=1, shadow=True)
plt.tight_layout()

plt.show()

Here's the plot of the result, where you can see that ω(0) = 0 and θ(100) = 0.

Note that the solution to the boundary value problem is not unique. If we modify the creation ystart to

ystart = odeint(pend, [np.pi, 0], t, args=(b, c,), tfirst=True)

a different solution is found, as seen in the following plot:

In this solution, the pendulum starts out almost at the inverted position (result.y[0, 0] = 3.141592653578858). It starts to fall very slowly; gradually it falls faster, and it reaches the straight down position at t = 100.

The trivial solution θ(t) ≡ 0 and ω(t) ≡ 0 also satisfies the boundary conditions.