Ordinary Differential Equations Exercise 1

Imports



In [117]:

    
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns
from scipy.integrate import odeint
from IPython.html.widgets import interact, fixed

Euler's method

Euler's method is the simplest numerical approach for solving a first order ODE numerically. Given the differential equation

$$ \frac{dy}{dx} = f(y(x), x) $$

with the initial condition:

$$ y(x_0)=y_0 $$

Euler's method performs updates using the equations:

$$ y_{n+1} = y_n + h f(y_n,x_n) $$$$ h = x_{n+1} - x_n $$

Write a function solve_euler that implements the Euler method for a 1d ODE and follows the specification described in the docstring:



In [126]:

    
def solve_euler(derivs, y0, x):
    """Solve a 1d ODE using Euler's method.
    
    Parameters
    ----------
    derivs : function
        The derivative of the diff-eq with the signature deriv(y,x) where
        y and x are floats.
    y0 : float
        The initial condition y[0] = y(x[0]).
    x : np.ndarray, list, tuple
        The array of times at which of solve the diff-eq.
    
    Returns
    -------
    y : np.ndarray
        Array of solutions y[i] = y(x[i])
    """
    y= np.ndarray(len(x))
    y[0]=y0
    h=x[1]-x[0]
    
    
    for i in range(len(x)-1):
        y[i+1]=y[i]+h*derivs(y[i],x[i])
        
        
    return y



In [127]:

    
assert np.allclose(solve_euler(lambda y, x: 1, 0, [0,1,2]), [0,1,2])

The midpoint method is another numerical method for solving the above differential equation. In general it is more accurate than the Euler method. It uses the update equation:

$$ y_{n+1} = y_n + h f\left(y_n+\frac{h}{2}f(y_n,x_n),x_n+\frac{h}{2}\right) $$

Write a function solve_midpoint that implements the midpoint method for a 1d ODE and follows the specification described in the docstring:



In [120]:

    
def solve_midpoint(derivs, y0, x):
    """Solve a 1d ODE using the Midpoint method.
    
    Parameters
    ----------
    derivs : function
        The derivative of the diff-eq with the signature deriv(y,x) where y
        and x are floats.
    y0 : float
        The initial condition y[0] = y(x[0]).
    x : np.ndarray, list, tuple
        The array of times at which of solve the diff-eq.
    
    Returns
    -------
    y : np.ndarray
        Array of solutions y[i] = y(x[i])
    """
    y= np.ndarray(len(x))
    h= x[1]-x[0]
    y[0]=y0
    for i in range(len(x)-1):
        y[i+1]= y[i]+h*derivs((y[i]+(h/2)*derivs(y[i],x[i])),x[i]+(h/2))
    return y



In [121]:

    
assert np.allclose(solve_midpoint(lambda y, x: 1, 0, [0,1,2]), [0,1,2])

You are now going to solve the following differential equation:

$$ \frac{dy}{dx} = x + 2y $$

which has the analytical solution:

$$ y(x) = 0.25 e^{2x} - 0.5 x - 0.25 $$

First, write a solve_exact function that compute the exact solution and follows the specification described in the docstring:



In [122]:

    
def solve_exact(x):
    """compute the exact solution to dy/dx = x + 2y.
    
    Parameters
    ----------
    x : np.ndarray
        Array of x values to compute the solution at.
    
    Returns
    -------
    y : np.ndarray
        Array of solutions at y[i] = y(x[i]).
    """
    y= .25*np.exp(2*x)-.5*x-.25
    return y



In [123]:

    
assert np.allclose(solve_exact(np.array([0,1,2])),np.array([0., 1.09726402, 12.39953751]))

In the following cell you are going to solve the above ODE using four different algorithms:

Euler's method
Midpoint method
odeint
Exact

Here are the details:

Generate an array of x values with $N=11$ points over the interval $[0,1]$ ($h=0.1$).
Define the derivs function for the above differential equation.
Using the solve_euler, solve_midpoint, odeint and solve_exact functions to compute the solutions using the 4 approaches.

Visualize the solutions on a sigle figure with two subplots:

Plot the $y(x)$ versus $x$ for each of the 4 approaches.
Plot $\left|y(x)-y_{exact}(x)\right|$ versus $x$ for each of the 3 numerical approaches.

Your visualization should have legends, labeled axes, titles and be customized for beauty and effectiveness.

While your final plot will use $N=10$ points, first try making $N$ larger and smaller to see how that affects the errors of the different approaches.



In [153]:

    
N= 11
h=0.1
x=np.linspace(0,1,N)


def derivs(x,y):
    return x+2*y

plt.figure(figsize=(10,15))

plt.subplot(2,1,1) #EULER's
plt.plot(x,solve_euler(derivs,0,x), color='r',label="Euler's") 

plt.subplot(2,1,1) # Midpoint
plt.plot(x,solve_midpoint(derivs,0,x), label ="Midpoint", marker='^')

plt.subplot(2,1,1) #ODEINT
plt.plot(x,odeint(derivs,0,x), color ='g', label= "ODEINT")

plt.subplot(2,1,1) #Exact
plt.plot(x,solve_exact(x),color='black',label="Exact")
plt.title('ODE solutions')
plt.xlabel('x')
plt.ylabel('y(x)')
plt.legend()

plt.subplot(2,1,2) # Euler- Exact
plt.plot(x,abs((solve_euler(derivs,0,x))-solve_exact(x)),color='r', label='Euler-Exact')


plt.subplot(2,1,2) #Midpoint - Exact
plt.plot(x,abs((solve_midpoint(derivs,0,x)-solve_exact(x))),color='b', label='Midpoint-Exact')

plt.subplot(2,1,2) #ODEINT-Exact
plt.plot(x,abs((odeint(derivs,0,x)-solve_exact(x))), color='g', label='ODEINT-Exact')

plt.title('Differences Between Methods and the Exact')
plt.xlabel('x')
plt.ylabel('y(x)')
plt.legend()

plt.tight_layout()

I HAVE NO IDEA WHAT IS WRONG WITH MY ODEINT-EXACT GRAPH (the green one)



In [125]:

    
assert True # leave this for grading the plots



In [ ]: