In [95]:
%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 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 [96]:
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])
"""
h=x[1]-x[0]
y=np.ndarray(len(x))
y[0]=y0
for n in range(len(x)-1):
y[n+1]=y[n]+h*derivs(y[n],x[n])
return y
In [97]:
solve_euler(lambda y, x: 1, 0, [0,1,2])
Out[97]:
In [98]:
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 [99]:
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.ndarray(len(x))
y[0]=y0
for n in range(len(x)-1):
y[n+1]=y[n]+h*derivs(y[n]+(h/2)*derivs(y[n],x[n]),x[n]+h/2)
return y
In [100]:
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 [101]:
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.ndarray(len(x))
for n in range(len(x)):
y[x[n]]=0.25*np.exp(2*x[n])-0.5*x[n]-0.25
return y
In [102]:
solve_exact(np.array([0,1,2]))
Out[102]:
In [103]:
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:
odeintHere are the details:
derivs function for the above differential equation.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:
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 [ ]:
In [105]:
x=np.array(range(0,11,1))/10
def derivs(y,x):
return x+2*y
u=odeint(derivs,0,x)
e=(u - solve_exact(x))
w=(solve_midpoint(derivs,0,x) - solve_exact(x))
q=(solve_euler(derivs,0,x) - solve_exact(x))
plt.figure(figsize=(10,10))
plt.subplot(211)
plt.plot(x,solve_euler(derivs,0,x))
plt.plot(x,solve_midpoint(derivs,0,x))
plt.plot(x,odeint(derivs,0,x))
plt.plot(x,solve_exact(x))
plt.subplot(212)
plt.plot(x,abs(q))
plt.plot(x,abs(w))
plt.plot(x,abs(e))
Out[105]:
In [67]:
assert True # leave this for grading the plots
In [ ]: