# Motivation for Implicit Time Integration
## CH EN 2450 - Numerical Methods
**Prof. Tony Saad (<a>www.tsaad.net</a>) <br/>Department of Chemical Engineering <br/>University of Utah**
<hr/>


```python
import numpy as np
from numpy import *
%matplotlib inline
%config InlineBackend.figure_format = 'svg'
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from scipy.optimize import fsolve
```

# Solve the following initial value problem using the Forward Euler method
Solve the ODE:
\begin{equation}
\frac{\text{d}y}{\text{d}t} = -ay;\quad y(0) = y_0
\end{equation}
The analytical solution is 
\begin{equation}
y(t) = y_0 e^{-a t}
\end{equation}

We first define the RHS for this ODE


```python
a = 1000.0
def rhs_sharp_transient(f, t):
    return - a * f
```

we now plot it


```python
yexact = lambda a,y0,t : y0 * np.exp(-a*t)
t = linspace(0,0.1,300)
y0 = 1
plt.plot(t,yexact(1000,y0,t),label='Exact')
plt.legend()
plt.grid()
# plt.savefig('sharp-transient-exact.pdf')
```


    

    


# Using Forward Euler
Using the forward Euler method, solve the above ODE using various timestep sizes, starting with $\Delta t = 0.003$.


```python
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
```


```python
f0 = 1.0
tend = 0.1
dt = 0.003
# dt = 2.0/a
# dt = 0.5*dt
print('dt',dt)
# a = 1000
t,ffe = forward_euler(rhs_sharp_transient,f0, tend, dt)
plt.plot(t,ffe,'r-o',label='Forward Euler, dt='+str(dt))

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

    dt 0.003



    

    


# Using Backward Euler


```python
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: a 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
```


```python
f0 = 1.0
tend = 0.1
dt = 0.0025

t,fbe = backward_euler(rhs_sharp_transient, f0, tend, dt)
plt.plot(t,fbe, 'r*-',label='Backward Euler')

texact = np.linspace(0,tend,400)
plt.plot(texact,yexact(a,f0,texact),label='Exact')
```




    [<matplotlib.lines.Line2D at 0x151b590190>]




    

    


# Using Crank-Nicolson


```python
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
```


```python
f0 = 1.0
tend = 0.1
dt = 0.0025

t,fcn = crank_nicolson(rhs_sharp_transient, f0, tend, dt)
plt.plot(t,fcn, 'r*-',label='Backward Euler')

texact = np.linspace(0,tend,400)
plt.plot(texact,yexact(a,f0,texact),label='Exact')
```




    [<matplotlib.lines.Line2D at 0x151c892cd0>]




    

    


# Using odeint


```python
time = np.linspace(0,0.1,20)
sol = odeint(rhs_sharp_transient,f0,time)
texact = np.linspace(0,tend,400)
plt.plot(texact,y(a,f0,texact),label='Exact')
plt.plot(time,sol, 'k-o', label='ODEint')
plt.legend()
```




    <matplotlib.legend.Legend at 0x152621ed50>




    

    


4e656db85f3b54f438a1f3df742a82a45dbadb79
ipynb
Jupyter Notebook
topics/initial-value-problems/Motivation for Implicit Methods.ipynb
jomorodi/NumericalMethods
e040693001941079b2e0acc12e0c3ee5c917671c
MIT
2019-03-27T05:22:34.000Z
2021-01-27T10:49:13.000Z
topics/initial-value-problems/Motivation for Implicit Methods.ipynb
jomorodi/NumericalMethods
e040693001941079b2e0acc12e0c3ee5c917671c
MIT
topics/initial-value-problems/Motivation for Implicit Methods.ipynb
jomorodi/NumericalMethods
e040693001941079b2e0acc12e0c3ee5c917671c
MIT
2019-12-29T23:31:56.000Z
2021-12-28T19:04:10.000Z
Qwen/Qwen-72B
1. YES
2. YES