# Ordinary Differential Equations Exercise 1

## Imports



In [1]:

%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




:0: FutureWarning: IPython widgets are experimental and may change in the future.



## 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 [ ]:




In [2]:

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])
"""
# YOUR CODE HERE
h=x[1]-x[0]
y=np.zeros((len(x)))
y[0]=y0
for i in range(1,len(x)):
n=i-1
y[i]=y[i-1]+h*derivs(x[n],y[n])

return y




In [3]:

solve_euler(lambda y, x: 1, 0, [0,1,2])




Out[3]:

array([ 0.,  1.,  2.])




In [4]:

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 [5]:

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])
"""
h=x[1]-x[0]
y=np.zeros((len(x)))
y[0]=y0

for i in range(1,len(x)):
n=i-1
y[i]=y[n]+h*derivs((y[n]+h/2*derivs(y[n],x[n])),x[n]+h/2)
return y




In [6]:

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 [7]:

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=np.zeros((len(x)))
for i in range(1,len(x)):
y[i]=.25*np.exp(2*x[i])-.5*x[i]-0.25
return y




In [8]:

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:

1. Euler's method
2. Midpoint method
3. odeint
4. 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:

1. Plot the $y(x)$ versus $x$ for each of the 4 approaches.
2. 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 [44]:

x=np.linspace(0,1,11)
def derivs(y,x): #is this right gotta check
dy=x+2*y
return dy
y4=solve_exact(x)
y1=solve_euler(derivs,0, x)
y2=solve_midpoint(derivs,0,x)
y3=odeint(derivs,0,x)




In [28]:

plt.plot?




In [45]:

plt.figure(figsize=(8,11))

plt.subplot(2,1,1)
plt.plot(x,y1,color='b',label='Euler')
plt.title('Approximations')
plt.xlabel('x')
plt.ylabel('y(x)')
plt.plot(x,y2,color='r',label='Midpoint')
plt.plot(x,y3,color='g',label='odeint',linestyle='--', )
plt.title('y(x) Calculations')
plt.plot(x,y4,color='k', linestyle='',marker='o',label='Exact')
plt.grid(False)
plt.legend(loc=4)

plt.subplot(2,1,2)
plt.plot(x,abs(y1-y4),color='b',label='Euler Error',linestyle='', marker='o')
plt.plot(x,abs(y2-y4),color='g',linestyle='',marker='^',label='Midpoint Error')
# plt.plot(x,abs(y3-y4),color='k',linestyle='',marker='>',label='Odeint vs Exact') #makes weird graph dont understand
plt.grid(False)
plt.title('Method Error compared to Exact')
plt.ylabel('Difference from Exact')
plt.xlabel('x')

plt.legend(loc=2)




Out[45]:

<matplotlib.legend.Legend at 0x7f0b8034a7b8>




In [ ]:

assert True # leave this for grading the plots