# Announcements
- No Problem Set this week, Problem Set 4 will be posted on 9/28.
- Stay on at the end of lecture if you want to ask questions about Problem Set 3.

<style>
@import url(https://www.numfys.net/static/css/nbstyle.css);
</style>
<a href="https://www.numfys.net"></a>

# Ordinary Differential Equations - higher order methods

<section class="post-meta">
Based on notes and notebooks by Niels Henrik Aase, Thorvald Ballestad, Vasilis Paschalidis and Jon Andreas Støvneng
</section>


## Algorithms for initial value problem ODEs


Assume we have a first-order differential equation which can be expressed in the form

$$ \frac{dy}{dt} = g(y,t) $$

We will solve this on constant-interval mesh of the independent variable $t$ defined by

$$ t_n = t_0 + n h $$

### Forward-Euler method

In Lecture 10 we derived Euler's method, which simply solves the first-order forward difference approximation to $dy/dt$

$$ \frac{y_{i+1}-y_i}{h} = g(y_i,t_i)$$

as

$$ y_{i+1} = y_i + h g(y_i,t_i) \label{Euler_fwd}\tag{3}$$


```python
# Importing the necessary libraries
import numpy as np # NumPy is used to generate arrays and to perform some mathematical operations
import matplotlib.pyplot as plt # Used for plotting results


```


```python
def forwardEuler_step(t, y, h, g, *P):
    """
    Implements a single step of the forward-Euler finite-difference scheme
    Parameters:
            t: time t
            y: Numerical approximation of y at time t
            h: Step size
            g: RHS of our ODE (RHS = Right hand side). Can be any function with signature g(t,y,*P).
            *P: tuple of parameters, arguments to g
        Returns:
            next_y: Numerical approximation of y at time t+h
    """
    next_y = y + h*g(t, y, *P)
    return next_y


```

We now need some sort of framework which will take this function and do the integration for us. Let's rewrite `full_Euler` from Lecture 10 to be more general:


```python
def odeSolve(t0, y0, tmax, h, g, method, *P):
    """ A full numerical aproximation of an ODE in a set time interval. Performs consecutive steps of `method`
    with step size h from start time until the end time. Also takes into account the initial values of the ODE
    
    Parameters:
            t0: start time
            y0 : Initial condition for y at t = t0
            tmax: The end of the interval where the `method` is integrated, t_N
            h: Step size
            g: RHS of our ODE (RHS = Right hand side). Can be any function with signature g(t,y,*P).
            *P: tuple of parameters, arguments to g
        Returns:
            t_list: Evenly spaced discrete list of time with spacing h. 
                    Starting time = start_t, and end time = end_t 
            y_list: Numerical approximation of y at times t_list
    """    
    # make the t-mesh; guarantees we stop precisely at tmax
    t_list = np.arange(t0,tmax+h,h)
    # allocate space for the solution
    y_list = np.zeros_like(t_list)
    # set the initial condition
    y_list[0] = y0
    
    # find out the size of the t-mesh, and then integrate forward one meshpoint per iteration of the loop
    n, = t_list.shape
    for i in range(0,n-1):
        y_list[i+1] = method(t_list[i], y_list[i], h, g, *P)
        
    # return the solution
    return t_list,y_list
```

Armed with this machinery, let's set up another simple problem and try it out.

Last time, we looked at exponential growth, let's solve exponential decay this time:

$$ \frac{dy}{dt} = - c y, \quad y[0] = 1 $$

First, we provide a function to implement the RHS:


```python
def expRHS(t, y, c):
    """
    Implements the RHS (y'(x)) of the DE
    """
    return -c*y
```

Now we set up the problem to compute and plot the result, along with a plot of the magnitude of the fractional error

### Runge-Kutta Schemes

The idea of the Runge-Kutta schemes is to take advantage of derivative information at the times between $t_i$ and $t_{i+1}$ to increase the order of accuracy.

For example, in the midpoint method, the derivative at the initial time is used to approximate the derivative at the midpoint of the interval, $f(y_i+\frac{1}{2}hf(y_i,t_i), t_i+\frac{1}{2}h)$. The derivative at the midpoint is then used to advance the solution to the next step. The method can be written in two stages $k_i$,

$$ \begin{aligned} \begin{array}{l} k_1 = h f(y_i,t_i)\\ k_2 = h f(y_i+\frac{1}{2}k_1, t_n+\frac{1}{2}h)\\ y_{i+1} = y_i + k_2 \end{array} \end{aligned}\label{RK2}\tag{4} $$

The midpoint method is known as a __2nd-order Runge-Kutta__ formula.


In general, an explicit 2-stage Runge-Kutta method can be written as,
$$ \begin{array}{l} k_1 = h f(y_n,t_n)\\ k_2 = h f(y_n+b_{21}k_1, t_n+a_2h)\ \\ y_{n+1} = y_n + c_1k_1 +c_2k_2 \label{explicitrk2}\tag{5}\end{array} $$

The scheme is said to be *explicit* since a given stage does not depend *implicitly* on itself, as in the backward Euler method , or on a later stage.

Other explicit second-order schemes can be derived by comparing Eq.(\ref{explicitrk2}) to other second-order expansions and matching terms to determine the coefficients $a_2$, $b_{21}$, $c_1$ and $c_2$.

### Explicit Fourth-Order Runge-Kutta Method

Explicit Runge-Kutta methods are popular as each stage can be calculated with one function evaluation. In contrast, implicit Runge-Kutta methods usually involves solving a non-linear system of equations in order to evaluate the stages. As a result, explicit schemes are much less expensive to implement than implicit schemes.

The higher-order Runge-Kutta methods can be derived by in manner similar to the midpoint formula. An s-stage method is compared to a Taylor method and the terms are matched up to the desired order.

As it happens to be, <strong>The Fourth Order Runge-Kutta Method</strong> uses three such test-points and is the most widely used Runge-Kutta Method. You might ask why we don't use five, ten or even more test-points, and the answer is quite simple: It is not computationally free to calculate all these test-points, and the gain in accuracy rapidly decreases beyond the fourth order of the method. That is, if high precision is of such importance that you would require a tenth-order Runge-Kutta, then you're better off reducing the step size $h$, than increasing the order of the method. 

Also, there exists other more sophisticated methods which can be both faster and more accurate for equivalent choices of $h$, but obviously, may be a lot more complicated to implement. See for instance <i>Richardson Extrapolation</i>, <i>the Bulirsch-Stoer method</i>, <i>Multistep methods, Multivalue methods</i> and <i>Predictor-Corrector methods</i>.


The classic fourth-order Runge-Kutta formula is:
$$ \begin{array}{l} k_1 = h f(y_n,t_n)\\ k_2 = h f(y_n+\frac{k_1}{2}, t_n+\frac{h}{2})\\ k_3 = h f(y_n+\frac{k_2}{2}, t_n+\frac{h}{2})\\ k_4 = h f(y_n+k_3, t_n+h)\\ y_{n+1} = y_n + \frac{k_1}{6}+ \frac{k_2}{3}+ \frac{k_3}{3} + \frac{k_4}{6} \label{RK4}\tag{6}\end{array} $$



```python
def RK2_step(t, y, h, g, *P):
    """
    Implements a single step of the second-order, explicit midpoint method
    """
    thalf = t + 0.5*h
    k1 = h * g(t, y, *P)
    k2 = h * g(thalf, y + 0.5*k1, *P)

    return y +k2
```


```python
def RK4_step(t, y, h, g, *P):
    """
    Implements a single step of a fourth-order, explicit Runge-Kutta scheme
    """
    thalf = t + 0.5*h
    k1 = h * g(t, y, *P)
    k2 = h * g(thalf, y + 0.5*k1, *P)
    k3 = h * g(thalf, y + 0.5*k2, *P)
    k4 = h * g(t + h, y + k3, *P)
    return y + (k1 + 2*k2 + 2*k3 + k4)/6
```


```python
# set up problem
c = 1.0
h = 0.5
t0 = 0.0
y0 = 1.0
tmax = 5.0

# call the solver for RK2
t, y = odeSolve(t0, y0, tmax, h, expRHS, RK2_step, c)

# plot the result
fig,ax = plt.subplots(1,2)
ans = np.exp(-c*t)
ax[0].plot(t,ans,'r')
ax[0].set_xlabel('t')
ax[0].set_ylabel('y')

ax[0].plot(t,y,'o','RK2')
err_RK2 = np.abs((ans-y)/ans)

# call the solver for Euler
t, y = odeSolve(t0, y0, tmax, h, expRHS, forwardEuler_step, c)
ax[0].plot(t,y,'o','Euler')
err = np.abs((ans-y)/ans)

# call the solver for RK2
t, y4 = odeSolve(t0, y0, tmax, h, expRHS, RK4_step, c)
ax[0].plot(t,y4,'o','RK4')
err_RK4 = np.abs((ans-y4)/ans)

# along with the errors
err_RK2 = np.abs((ans-y)/ans)
ax[1].plot(t, err_RK2, 'o',label = "RK2")
ax[1].plot(t, err_RK4, 'o',label = "RK4")
ax[1].plot(t, err, 'o',label = "Euler")
ax[1].set_xlabel('t')
ax[1].set_ylabel('fractional error')
ax[1].legend()
# now also overplot the error we calculated for forward-Euler

# this gives better spacing between axes
plt.tight_layout()
plt.show()
```

### Systems of First-Order ODEs

Next, we turn to systems of ODE's. We'll take as our example the Lotke-Volterra equations, a simple model of population dynamics in an ecosystem (with many other uses as well).

Imagine a population of rabbits and of foxes on a small island. The rabbits eat a plentiful supply of grass and
would breed like, well, rabbits, with their population increasing exponentially with time in the absence of preditors. The foxes eat the rabbits, and would die out exponentially in time with no food supply. The rate at which foxes eat rabbits depends upon the product of the fox and rabbit populations.

The equations for the population of the rabbits $R$ and foxes $F$ in this simple model is then

\begin{eqnarray*}
\frac{dR}{dt} &= \alpha R - \beta R F \\
\frac{dF}{dt} &= \delta R F - \gamma F
\end{eqnarray*}

Without the cross terms in $RF$, these are just two decay equations of the form we have used as an example above.

A random set of parameters (I am not a biologist!) might be that a rabbit lives four years, so $\alpha=1/4$ and
a fox lives 10 years, so $\gamma=1/10$. Let's pick the other parameters as $\beta = 1$ and $\delta = 1/4$.

We can express the unknown populations as a vector of length two: $y = (R, F)$. The rate of change of populations then can also be expressed as a vector $dy/dt = (dR/dt, DF/dt)$. With such a definition, we can write the RHS function of our system as


```python
def lvRHS(t, y, *P):
    # Lotke-Volterra system RHS
    
    # unpack the parameters from the array P
    alpha, beta, gamma, delta = P

    # make temporary variables with rabbit and fox populations
    R = y[0]
    F = y[1]
    
    # LV system
    dRdt = alpha * R - beta * R * F
    dFdt = delta * R * F - gamma * F
    
    # return an array of derivatives with same order as input vector
    return np.array([ dRdt, dFdt ])
```

We now have to generalize our odeSolve function to allow more than one equation


```python
def odeSolve(t0, y0, tmax, h, RHS, method, *P):
    """
    ODE driver with constant step-size, allowing systems of ODE's
    """
    # make array of times and find length of array
    t = np.arange(t0,tmax+h,h)
    ntimes,  = t.shape

    # find out if we are solving a scalar ODE or a system of ODEs, and allocate space accordingly
    
    if type(y0) in [int, float]:  # check if primitive type -- means only one eqn
        neqn = 1
        y = np.zeros( ntimes )
    else:                         # otherwise assume a numpy array -- a system of more than one eqn
        neqn, = y0.shape
        y = np.zeros( (ntimes, neqn) )

    # set first element of solution to initial conditions (possibly a vector) 
    y[0] = y0

    # march on...
    for i in range(0,ntimes-1):
        y[i+1] = method(t[i], y[i], h, RHS, *P)

    return t,y

```

Now we can solve our system of two coupled ODEs. Note that the solution is now a vector of 2D vectors... the first index is the solution time, the second the variable:


```python
alpha = 1.0 
beta = 0.025
gamma = 0.4
delta = 0.01

h = 0.2
t0 = 0.0
y0 = np.array([ 30, 10 ])
tmax = 50

# call the solver
t, y = odeSolve(t0, y0, tmax, h, lvRHS, RK4_step, alpha, beta, gamma, delta)

fig,ax = plt.subplots()
ax.plot(t,y[:,0],'b', label='prey')
ax.plot(t,y[:,1],'r', label='preditor')
ax.set_xlabel('time')
ax.set_ylabel('population')
ax.legend()

plt.tight_layout()
plt.show()
```

### Higher Order Derivatives and Sets of 1st order ODEs

The trick to solving ODEs with higher derivatives is turning them into systems of first-order ODEs.

As a simple example, consider the second-order differential equation describing the van der Pol oscillator

$$ \frac{d^2 x}{dt^2} - a (1-x^2) \frac{dx}{dt} + x = 0 $$

We turn this into a pair of first-order ODEs by defining an auxiliary function $v(t) = dx/dt$ and writing the system as

\begin{align}
\begin{split}
\frac{dv}{dt} &= a (1-x^2) v - x\\
\frac{dx}{dt} &= v
\end{split}
\end{align}

Note that there are only functions (and the independent variable) on the RHS; all "differentials" are on the LHS.

Now that we have a system of first-order equations ,we can proceed as above. A function describing
the RHS of this system is


```python
def vdpRHS(t, y, a):
    # we store our function as the array [x, x']
    return np.array([
        
        y[1],                          # dx/dt = v
        a*(1-y[0]**2)*y[1] - y[0]      # dv/dt = a*(1-x**2)*v - x
        
    ])
```


```python
a = 15 # parameter

h = 0.01
t0 = 0.0
y0 = np.array([ 0,  1])
tmax = 50

# call the solver
t, y = odeSolve(t0, y0, tmax, h, vdpRHS, RK4_step, a)

fig,ax = plt.subplots()
ax.plot(t,y[:,0],'b', label='x')
ax.plot(t,y[:,1],'r--', label='v')
ax.set_xlabel('time')
ax.legend()
ax.set_title(f"van der Pol Oscillator for a={a}")

plt.tight_layout()
plt.show()

```

A somewhat more complex example is the Lane-Emden equation, which is really just Poisson's equation in spherical symmetry for the graviational potential of a self-gravitating fluid whose pressure is related to its density as $P\propto\rho^\gamma$. Such a system is called a _polytrope_, and is often used in astrophysics as a simple model for the structure of system such as a a star in which outward pressure and inward gravity are in equilibrium.

Let $\xi$ be the dimensionless radius of the system, and let $\theta$ be related to the density as
$\rho = \rho_c \theta^n$, where $\rho_c$ is the density at the origin and $n = 1/(\gamma-1)$. We then have the dimensionless second-order differential equation

$$ \frac{1}{\xi^2}\frac{d}{d\xi}\left(\xi^2\frac{d\theta}{d\xi}\right) + \theta^n = 0 $$

Note that the first term is just the divergence $\nabla\cdot\theta$ in spherical symmetry.

If we expand out the first term, we have 

$$ \frac{d^2\theta}{d\xi^2} + \frac{2}{\xi}\frac{d\theta}{d\xi} + \theta^n = 0 $$

Defining an auxiliary function $v(\xi) = d\theta/d\xi$, we can then convert this into a system of two first-order ODEs:

\begin{align}
\begin{split}
\frac{dv}{d\xi} &= -\frac{2}{\xi} v - \theta^n \\
\frac{d\theta}{d\xi} &= v
\end{split}
\end{align}

Again, we have "derivatives" only on the LHS and no derivatives on the RHS of our system.

Looking at this expression, one can see right away that at the origin $\xi=0$ we will have a numerical problem; we are dividing by zero.
Analytically, this is not a problem, since $v/\xi\rightarrow0$ as $\xi\rightarrow0$, but here we need to address this numerically.

The first approach is to take care of the problem in our RHS function:


```python
def leRHS(x, y, n):
    
    dthetadx = y[1]
    
    if x==0:
        dvdx = -y[0]**n
    else:
        dvdx = -2/x*y[1] - y[0]**n
    
    return np.array([ dthetadx, dvdx ])
```

This is somewhat clunky, however, and you would first have to convince yourself that in fact $v(\xi)\rightarrow0$ faster than $\xi$ (don't just take my word for it!).

Instead, we could use a more direct RHS function


```python
def leRHS(x, y, n):
    
    dthetadx = y[1]
    
    dvdx = -2/x*y[1] - y[0]**n
    
    return np.array([ dthetadx, dvdx ])
```

and expand the solution in a Taylor series about the origin to get a starting value for our numerical integration at a small distance away from the origin. To do this, write 

$$\theta(\xi) = a_0 + a_1 \xi + a_2 \xi^2 + \dots$$

The first thing to notice is that, by symmetry, only even powers of $\xi$ will appear in the solution.
Thus we will have

$$ \theta(\xi) = a_0 + a_2 \xi^2 + a_4 \xi^4 + \dots$$

By the boundary condition $\theta(0) = 1$, we have immediately that $a_0 = 1$.

Next, substitute $\theta(\xi) = 1 + a_2 \xi^2 + a_4 \xi ^4 + O(\xi^6)$ into the Lane-Emden equation. $\theta$ and its first two derivatives are

\begin{align}
\begin{split}
\theta(\xi) &= 1 + a_2 \xi^2 + a_4 \xi^4 + O(\xi^6)\\
\theta'(\xi) &= 2 a_2 \xi + 4 a_4 \xi^3 + O(\xi^5) \\
\theta''(\xi) &= 2 a_2 + 12 a_4 \xi^2 + O(\xi^4)
\end{split}
\end{align}

Putting these into the Lane-Emden equation, we have
\begin{align}
\begin{split}
2 a_2 + 12 a_4 \xi^2 + O(\xi^4) + \frac{2}{\xi} (2 a_2 x + 4 a_4 \xi^3 + O(\xi^5))  &= -\theta^n \\
6 a_2 + 20 a_4 \xi^2 + O(\xi^4) &= -\theta^n
\end{split}
\end{align}

A boundary condition $\theta(0)=1$, and thus we have $a_2 = -1/6$. Away from zero, then, we have

\begin{align}
\begin{split}
-1 + 20 a_4 \xi^2 + O(\xi^4) &= -\left(1 - 1/6 \xi^2 + a_4 \xi^4 + O(\xi^6)\right)^n
\end{split}
\end{align}

The term on the RHS is $ 1 - n \xi^2/6 + O(\xi^4)$, and so we must have $a_4 = n/120$.

Thus, the series expansion of the solution around the origin is

$$ \theta(\xi) = 1 - \frac{1}{6}\xi^2 + \frac{n}{120} \xi^4 + \dots $$

We can now use this expansion to take a first step slightly away from the origin before beginning our
numerical integration, thus avoiding the divide by zero. Note that this series solution near the origin is $O(h^5)$ and so is a good match for RK4 if we take the same (or smaller) step-size.


```python
n = 3

xi0 = 0.01                            # starting value of xi for our numerical integration
theta0 = 1 - xi0**2/6 + n*xi0**4/120  # Taylor series solution to the DE near zero derived above
theta0p = -xi0/3 + n*xi0**3/30
y0 = np.array([ theta0,  theta0p])    # set IC's for numerical integration
print(f"IC at {xi0:10.5e}: {y0[0]:10.5e}, {y0[1]:10.5e}")

h = 0.1
tmax = 8

# call the solver
t, y = odeSolve(xi0, y0, tmax, h, leRHS, RK4_step, n)

fig,ax = plt.subplots()
ax.plot(t,y[:,0],'b', label=r'$\theta(\xi)$')
ax.plot(t,y[:,1],'r--', label=r'$\frac{d\theta}{d\xi}$')
ax.plot([0,tmax],[0,0],'k')
ax.set_xlabel(r'$\xi$')
ax.set_title(f"Lane Emden Equation for n={n}")
ax.legend()

plt.tight_layout()
plt.show()
```

For values of $n\le5$, the solutions of the Lane Emden equation (the so-called Lane-Emden functions of index $n$) decrease to zero at finite $\xi$. Since this is the radius at which the density goes to zero, we can interpret it as the surface of the self-gravitating body (for example, the radius of the star). Knowing this value $\xi_1$ is thus interesting... Let us see how to determine it numerically.

Cleary, we are looking for the solution to $\theta(\xi_1)=0$; this is just root-finding, which we already know how to do. Instead of using some closed-form function, however, the value of the function $\theta(\xi)$ must in this case be determined numerically. But we have just figured out how to do this!

Let's use the bisection method for our root-finding algorithm; here is a quick version (no error checking!)


```python
def bisection(func, low, high, eps, *P):
    flow = func(low, *P)
    fhigh = func(high, *P)
    mid = 0.5*(low+high)
    fmid = func(mid,*P)
    
    while (high-low)> eps:
        if fmid*flow < 0:
            high = mid
            fhigh = fmid
        else:
            low = mid
            flow = mid
            
        mid = 0.5*(low+high)
        fmid = func(mid,*P)
            
    return low
```

Now let us make a function which returns $\theta(\xi)$, the solution to the Lane-Emden equation at $\xi$


```python
def theta(xi, n):
    h = 1e-4
    xi0 = 1e-4
    theta0 = 1 - xi0**2/6 + n*xi0**4/120
    theta0p = -xi0/3 + n*xi0**3/30
    y0 = np.array([ theta0,  theta0p])
    t, y = odeSolve(xi0, y0, xi, h, leRHS, RK4_step, n)
    
    return y[-1,0]
```

Using these, we can compute the surface radius of the polytrope


```python
n = 3
xi1 = bisection(theta, 6, 8, 1e-5, n)
print(f"xi_1 = {xi1:7.5f}")
```

A more careful treatment gives a value $\xi_1 = 6.89685...$, so we are doing pretty well...


```python

```

30c00abbaaa3111abafe96512375232710a15b33
ipynb
Jupyter Notebook
Lectures/Lecture 12/Lecture12_ODE_part3.ipynb
astroarshn2000/PHYS305S20
18f4ebf0a51ba62fba34672cf76bd119d1db6f1e
MIT
2020-09-10T06:45:46.000Z
2020-10-20T13:50:11.000Z
Lectures/Lecture 12/Lecture12_ODE_part3.ipynb
astroarshn2000/PHYS305S20
18f4ebf0a51ba62fba34672cf76bd119d1db6f1e
MIT
Lectures/Lecture 12/Lecture12_ODE_part3.ipynb
astroarshn2000/PHYS305S20
18f4ebf0a51ba62fba34672cf76bd119d1db6f1e
MIT
Qwen/Qwen-72B
1. YES
2. YES