Integrating an ODE using odeint

consider the toy ODE boundary value problem:

$$\frac{d^2\Psi}{dx^2} + [100-\beta]\Psi = 0, x\in[-1,1]$$

where we have the boundary conditions of $$\Psi(1) = 0\\ \Psi(-1) = 0$$

I will adopt the notation $\Psi_{XX}$ to denote the second derivative through the rest of this notebook.

This is a second order system (second derivative), but we only know how to solve first order ODEs via Runge Kutta or such. Where 'we' = all mathematicians in the world. It is the one trick we know how to do on computers.

Hence, let's reduce this to a system of first order equations, and then solve them.

Define some new variables:

$$y_1 = \Psi \\ y_2 = \Psi_X$$

Now take their derivative to get first order equations:

$$\begin{aligned}y'_1 &= y2\\ y'_2 &= \Psi_{XX} \end{aligned}$$

Now substitute for $Y_{XX}$ using the original equation so that we have eliminated the second order terms from the equations:

$$\begin{aligned}y'_1 &= y2\\ y'_2 &= [\beta-100]\Psi\end{aligned}$$

Now that is done, we can pass this system of first order equations into scipy's odeint function to solve these equations.

Next problem: odeint just uses Runge Kutta (or a similar integrator) to iterate from the first value to the last value. The boundary conditions give our start point, but we have no way to specify the ending point.

Let's look at an example:

That's a fair amount of code to throw at you. The heart of the problem is that you have to write an evaluation function that takes a tuple of the current value of the system equations, time, and any additional arguments that you might need. You evaluate the system equations for that time step, and return them in a tuple. We do this in the function f().

def f(ys, t, beta):
    y1 = ys[0]
    y2 = ys[1]
    return [y2, (beta-100)*y1]

The rest of the code is mostly boilerplate. This gives us a range of values for x. I use 100 points just so we can plot the differential equation with a fair amount of fidelity; it does not affect the precision of the ode integration.

xs = np.linspace(-1, 1, 100)

We need to define the values for the system of equations. We have two equations, so we will need a tuple with two elements in it. We have no idea what the correct values are, so we fairly arbitrarily choose 0 and 1.

system_eq = (0,1)

Finally, we have to choose $\beta$. Our equation is $$\Psi_{XX} = [\beta-100]\Psi$$. This is one of the very few ODE's we can solve analytically, and we know the answer is $e^{(100-B)X}$. The power will be positive if $\beta >100$, and that will never reach our second boundary condition of $\Psi(1) = 0$. Therefore we set $\beta$ to an arbitrary value less than 100.

Finally, we integrate using the function scipy.integrate.odeint. The first parameter is the name of the function that performs the evaluation on the system of equations. The second parameter is the initial conditions for the system of equations. The third parameter is the range of values for $X$. Our function needs to take $\beta$ as an input, so we use the optional args parameter to provide beta to the function.

z = integrate.odeint(f, system_eq, xs, args=(beta,))

Shooting Method to find second boundary condition

So we have integrated our ODE, but we can see that we did not achieve the boundary condition $\Psi(1)=0$. Instead, it looks like $\Psi(1)\approx 0.9$.

All that is left is to vary beta to find the value for which $\Psi(1)=0$. We will define an initial step size for $\beta$

d_beta = 1

and then write a for loop that searches for the correct value.

This is a really naive implementation because the ODE has multiple solutions, and we only find the first one, but the point was just to document the use of odeint, not to make an industrial strength solver. For example, this is what the equation looks like with $\beta=40$