Author : Ninad Joshi, email: njoshi@fias.uni-frankfurt.de, webpage: fias.uni-frankfurt.de/~njoshi
Numerical Solutions to ODEs¶
In [1]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import math
Consider following equation
\begin{equation} \frac{dy}{dt} = f (y,t) \end{equation}where,
\begin{equation} f(t,y) = y - t^2 + 1 \end{equation}within interval $[0,2]$ with $f(0) =0.5$
In [2]:
def f( t, y ):
return y - t**2 + 1
The analytical solution is
$y(t) = (t+1) ^2 - 0.5~e^t$
In [3]:
def Exact(t):
return (t+1.0)**2.0-0.5*np.exp(t)
Euler Method:
\begin{equation} y_{n+1} =y_{n} + h ~ f(t_n,y_n) \end{equation}$h$ is step size
In [4]:
def Euler_method( f, a, b, h, y0 ):
n = int( (b-a)/h ) # determine step-size
t = np.arange( a, b+h, h ) # create equidistant time/ x-axis
w = []
t[0] = a
w.append (y0 )
for i in range(1,n+1): # apply Eulers's Method
w. append ( w[i-1] + h * f( t[i-1], w[i-1] ) )
return np.asarray(w)
Heun Method:
\begin{equation} \tilde{y}_{n+1} =y_{n} + h ~ f(t_n,y_n) \\ y_{n+1} =y_{n} + \frac{h}{2} \Big[ f(t_n,y_n)+f(t_{n+1}, \tilde{y}_{n+1} ) \Big] \end{equation}$h$ is step size
In [5]:
def Heun_method( f, a, b, h, y0):
n = int( (b-a)/h ) # determine step-size
t = np.arange( a, b+h, h ) # create equidistant time/ x-axis
w = []
t[0] = a
w.append (y0 )
for i in range(1,n+1): # apply Heun's Method
f1 = f( t[i-1], w[i-1] )
f2 = f( t[i-1] + h, w[i-1] + h * f1 )
w. append ( w[i-1] + h * ( f1 + f2 ) *0.5 )
return np.asarray(w)
In [6]:
# limits: 0.0 <= t <= 2.0
a = 0.0
b = 2.0
# steps
N = 10
# step-size
h = (b-a)/N
# initial value: y(0.0) = 0.5
y0 = 0.5
# equidistant time/ x-axis
t = np.arange( a, b+h, h )
result1 = Exact(t)
result2 = Euler_method( f, a, b, h, y0 )
result3 = Heun_method( f, a, b, h, y0 )
plt.plot( t, result1, label= 'Exact', linestyle='-', color='red' )
plt.plot( t, result2, label='Euler', linestyle='-.', color='green')
plt.plot( t, result3, label='Heun',linestyle=':', color='blue')
plt.xlabel('t')
plt.ylabel('y(t)')
plt.legend(loc=4)
plt.grid()
Euler Method:
\begin{equation} \tilde{y}_{n+1} =y_{n} + h ~ f(t_n,y_n) \\ y_{n+1} =y_{n} + \frac{h}{2} \Big[ f(t_n,y_n)+f(t_{n+1}, \tilde{y}_{n+1} ) \Big] \end{equation}$h$ is step size
In [7]:
def Taylor_method( f, a, b, h, y0 ):
n = int( (b-a)/h ) # determine step-size
t = np.arange( a, b+h, h ) # create equidistant time/ x-axis
w = []
t[0] = a
w.append (y0 )
for i in range(1, n+1): # Apply Taylor Method
T = 0
for j in range(len(f)):
h_factor = h**(j)/float(math.factorial(j+1))
T += h_factor * f[j]( t[i-1], w[i-1] )
w.append (w[i-1] + h * T)
return w
In [8]:
a, b = 0.0, 2.0
N = 10
# step-size
h = (b-a)/N
# initial value: y(0.0) = 0.5
y0 = 0.5
t = np.arange( a, b+h, h )
f_ = lambda t, y: y - t**2.0 + 1
df_ = lambda t, y: y - t**2.0 + 1 - 2*t
ddf_ = lambda t, y: y - t**2.0 - 2*t - 1
result4 = Taylor_method( [ f_, df_ ], a, b, h, y0 )
result5 = Taylor_method( [ f_, df_, ddf_ ], a, b, h, y0 )
plt.plot( t, result4, label='2nd order', linestyle=':', color='blue' )
plt.plot( t, result5, label='3rd order', linestyle='-.', color='green' )
plt.plot( t, Exact(t), label='exact', linestyle='-', color='red' )
plt.title( "Taylor's Method Example, N="+str(N) )
plt.title( "Taylor's Method Example, N="+str(N) )
plt.xlabel('t')
plt.ylabel('y(t)')
plt.legend(loc=4)
plt.grid()
In [9]:
def Midpoint_method( f, a, b, h, y0 ):
n = int( (b-a)/h )
t = np.arange( a, b+h, h )
w = []
t[0] = a
w.append (y0 )
for i in range(1,n+1): # Apply Midpoint Method
w.append ( w[i-1] + h * f( t[i-1] + h/2.0, w[i-1] + h * f( t[i-1], w[i-1] ) / 2.0 ) )
return w
In [10]:
a, b = 0.0, 2.0
N = 10
# step-size
h = (b-a)/N
# initial value: y(0.0) = 0.5
y0 = 0.5
t = np.arange( a, b+h, h )
f_ = lambda t, y: y - t**2.0 + 1
result6 = Midpoint_method( f_, a, b, h, y0 )
plt.plot( t, result6, label='midpoint', linestyle='-.', color='green' )
plt.plot( t, Exact(t), label='exact', linestyle='-', color='red' )
plt.title( 'Midpoint Method' )
plt.xlabel('t')
plt.ylabel('y(t)')
plt.legend(loc=4)
plt.grid()
In [11]:
def Runge_Kutta_4_method( f, a, b, h, y0 ):
n = int( (b-a)/h )
t = np.arange( a, b+h, h )
w = []
t[0] = a
w.append (y0 )
for i in range(1,n+1): # Apply Fourth Order Runge-Kutta Method
k1 = h * f( t[i-1], w[i-1] )
k2 = h * f( t[i-1] + h / 2.0, w[i-1] + k1 / 2.0 )
k3 = h * f( t[i-1] + h / 2.0, w[i-1] + k2 / 2.0 )
k4 = h * f( t[i], w[i-1] + k3 )
w.append ( w[i-1] + ( k1 + 2.0 * k2 + 2.0 * k3 + k4 ) / 6.0 )
return w
In [12]:
a, b = 0.0, 2.0
N = 10
# step-size
h = (b-a)/N
# initial value: y(0.0) = 0.5
y0 = 0.5
t = np.arange( a, b+h, h )
f_ = lambda t, y: y - t**2.0 + 1
result7 = Runge_Kutta_4_method( f_, a, b, h, y0 )
plt.plot( t, result7, label='Runge Kutta 4rth order', linestyle='-.', color='green' )
plt.plot( t, Exact(t), label='exact', linestyle=':', color='red' )
plt.title( 'Midpoint Method' )
plt.xlabel('t')
plt.ylabel('y(t)')
plt.legend(loc=4)
plt.grid()
For adaptive step size, we first calculate the integral with a arbritrary step size $h$ and $h/2$. The difference between area is compared against tolerancence. If the error is larger than the tolerance the step size is decreased to the half of its value. This cycle is then repeated till desired tolarance is achieved.
In [13]:
def AdaptDiff(Interval, y0, h, tol):
a = Interval[0] ;
b = Interval[1] ;
print( 'Initial h=',h)
while True:
b1 = a + 20.* h
if (b1> b):
print('initial interval too large ')
return
y1 = Heun_method(f, a, b1, h, y0 )
y2 = Heun_method(f, a, b1, h/2, y0)
err = np.absolute(y2[-1]-y1[-1])
if np.min(err) < tol:
print( 'Final h=',h)
result = Heun_method(f, a, b, h, y0)
break
if np.min(err) < tol/4:
h = h*3.0
else:
h = h*0.5
return(result , h )
In [14]:
tol =0.0001
a, b = 0.0, 2.0
N = 20
# step-size
h = (b-a)/N
y0 = np.array([0.5])
result8, h = AdaptDiff([0.0,2], y0, h, tol)
tf = np.arange( 0, 2+h,h )
plt.plot( tf, result8, label='Heun Method', linestyle='-.', color='green' )
plt.plot( tf, Exact(tf), label='exact', linestyle=':', color='red' )
plt.title( 'Adaptive Stepsize' )
plt.xlabel('t')
plt.ylabel('y(t)')
plt.legend(loc=4)
plt.grid()
In [15]:
def AdaptDiff_Euler(Interval, y0, h, tol):
a = Interval[0] ;
b = Interval[1] ;
print( 'Initial h=',h)
while True:
b1 = a + 20.* h
if (b1> b):
print('initial interval too large ')
return
y1 = Euler_method(f, a, b1, h, y0 )
y2 = Euler_method(f, a, b1, h/2, y0)
err = np.absolute(y2[-1]-y1[-1])
if np.min(err) < tol:
print( 'Final h=',h)
result = Euler_method(f, a, b, h, y0)
break
if np.min(err) < tol/4:
h = h*3.0
else:
h = h*0.5
return(result , h )
In [16]:
tol =0.0001
a, b = 0.0, 2.0
N = 20
# step-size
h = (b-a)/N
y0 = np.array([0.5])
result9, h = AdaptDiff([0.0,2], y0, h, tol)
tf = np.arange( 0, 2+h,h )
plt.plot( tf, result9, label='Euler Method_adaptive', linestyle='-.', color='blue' )
plt.plot( t, result2, label='Euler Method', linestyle='-.', color='green' )
plt.plot( tf, Exact(tf), label='exact', linestyle=':', color='red' )
plt.title( 'Adaptive Stepsize' )
plt.xlabel('t')
plt.ylabel('y(t)')
plt.legend(loc=4)
plt.grid()
