# Motivation for Implicit Time Integration¶

## CH EN 2450 - Numerical Methods¶

Department of Chemical Engineering
University of Utah

In [1]:
import numpy as np
from numpy import *
%matplotlib inline
%config InlineBackend.figure_format = 'svg'
import matplotlib.pyplot as plt
from scipy.integrate import odeint

In [2]:
def forward_euler(rhs, f0, tend, dt):
''' Computes the forward_euler method '''
nsteps = int(tend/dt)
f = np.zeros(nsteps)
f[0] = f0
time = np.linspace(0,tend,nsteps)
for n in np.arange(nsteps-1):
f[n+1] = f[n] + dt * rhs(f[n], time[n])
return time, f

def be_residual(fnp1, rhs, fn, dt, tnp1):
'''
Nonlinear residual function for the backward Euler implicit time integrator
'''
return fnp1 - fn - dt * rhs(fnp1, tnp1)

def backward_euler(rhs, f0, tend, dt):
'''
Computes the backward euler method
:param rhs: an rhs function
'''
nsteps = int(tend/dt)
f = np.zeros(nsteps)
f[0] = f0
time = np.linspace(0,tend,nsteps)
for n in np.arange(nsteps-1):
fn = f[n]
tnp1 = time[n+1]
fnew = fsolve(be_residual, fn, (rhs, fn, dt, tnp1))
f[n+1] = fnew
return time, f

def cn_residual(fnp1, rhs, fn, dt, tnp1, tn):
'''
Nonlinear residual function for the Crank-Nicolson implicit time integrator
'''
return fnp1 - fn - 0.5 * dt * ( rhs(fnp1, tnp1) + rhs(fn, tn) )

def crank_nicolson(rhs,f0,tend,dt):
nsteps = int(tend/dt)
f = np.zeros(nsteps)
f[0] = f0
time = np.linspace(0,tend,nsteps)
for n in np.arange(nsteps-1):
fn = f[n]
tnp1 = time[n+1]
tn = time[n]
fnew = fsolve(cn_residual, fn, (rhs, fn, dt, tnp1, tn))
f[n+1] = fnew
return time, f


# Solve the following initial value problem using the Forward Euler method¶

Solve the ODE: $$\frac{\text{d}y}{\text{d}t} = -ay;\quad y(0) = y_0$$ The analytical solution is $$y(t) = y_0 e^{-a t}$$

In [3]:
y = lambda a,y0,t : y0 * np.exp(-a*t)
t = linspace(0,0.1,300)
y0 = 1
# plt.plot(t,y(1,y0,t),label='a = 1')
# plt.plot(t,y(10,y0,t),label='a = 10')
plt.plot(t,y(1000,y0,t),label='Exact')
plt.legend()
plt.grid()
plt.savefig('sharp-transient-exact.pdf')

In [4]:
a = 1000.0
def rhs_sharp_transient(f, t):
return - a * f
f0 = 1.0
tend = 0.1
dt = 0.0005
# dt = 2.0/a
# dt = 0.5*dt
print('dt',dt)
# a = 1000
t,ffe = forward_euler(rhs_sharp_transient,f0, tend, dt)
# t,fbe = backward_euler(rhs_sharp_transient, f0, tend, dt)
# t, fcn = crank_nicholson(rhs_sharp_transient,f0,tend,dt)
plt.plot(t,ffe,'r-o',label='Forward Euler, dt='+str(dt))
# plt.plot(t,fbe, 'r*-',label='Backward Euler')
# plt.plot(t,fcn, 'g*-',label='Crank-Nicholson')

texact = np.linspace(0,tend,400)
plt.plot(texact,y(a,f0,texact),label='Exact')
plt.legend()
plt.grid()
plt.savefig('motivation for implicit dt='+str(dt)+'.pdf')

dt 0.0005

In [5]:
time = np.linspace(0,0.1,100)
sol = odeint(rhs_sharp_transient,f0,time)
plt.plot(time,sol, 'k-o', label='ODEint')
plt.legend()

Out[5]:
<matplotlib.legend.Legend at 0x1514125320>