{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Motivation for Implicit Time Integration\n",
"## CH EN 2450 - Numerical Methods\n",
"**Prof. Tony Saad (www.tsaad.net)
Department of Chemical Engineering
University of Utah**\n",
"
"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"from numpy import *\n",
"%matplotlib inline\n",
"%config InlineBackend.figure_format = 'svg'\n",
"import matplotlib.pyplot as plt\n",
"from scipy.integrate import odeint\n",
"from scipy.optimize import fsolve"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Solve the following initial value problem using the Forward Euler method\n",
"Solve the ODE:\n",
"\\begin{equation}\n",
"\\frac{\\text{d}y}{\\text{d}t} = -ay;\\quad y(0) = y_0\n",
"\\end{equation}\n",
"The analytical solution is \n",
"\\begin{equation}\n",
"y(t) = y_0 e^{-a t}\n",
"\\end{equation}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We first define the RHS for this ODE"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"a = 1000.0\n",
"def rhs_sharp_transient(f, t):\n",
" return - a * f"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"we now plot it"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"yexact = lambda a,y0,t : y0 * np.exp(-a*t)\n",
"t = linspace(0,0.1,100)\n",
"y0 = 1\n",
"plt.plot(t,yexact(1000,y0,t),label='Exact')\n",
"plt.legend()\n",
"plt.grid()\n",
"# plt.savefig('sharp-transient-exact.pdf')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Using Forward Euler\n",
"Using the forward Euler method, solve the above ODE using various timestep sizes, starting with $\\Delta t = 0.003$."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"def forward_euler(rhs, f0, tend, dt):\n",
" ''' Computes the forward_euler method '''\n",
" nsteps = int(tend/dt)\n",
" f = np.zeros(nsteps)\n",
" f[0] = f0\n",
" time = np.linspace(0,tend,nsteps)\n",
" for n in np.arange(nsteps-1):\n",
" f[n+1] = f[n] + dt * rhs(f[n], time[n])\n",
" return time, f"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"dt 0.0005\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"f0 = 1.0\n",
"tend = 0.1\n",
"dt = 0.0005\n",
"# dt = 2.0/a\n",
"# dt = 0.5*dt\n",
"print('dt',dt)\n",
"# a = 1000\n",
"t,ffe = forward_euler(rhs_sharp_transient,f0, tend, dt)\n",
"plt.plot(t,ffe,'r-o',label='Forward Euler, dt='+str(dt))\n",
"\n",
"texact = np.linspace(0,tend,400)\n",
"plt.plot(texact,yexact(a,f0,texact),label='Exact')\n",
"plt.legend()\n",
"plt.grid()\n",
"plt.savefig('motivation for implicit dt='+str(dt)+'.pdf')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Using Backward Euler"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"def be_residual(fnp1, rhs, fn, dt, tnp1):\n",
" '''\n",
" Nonlinear residual function for the backward Euler implicit time integrator\n",
" ''' \n",
" return fnp1 - fn - dt * rhs(fnp1, tnp1)\n",
"\n",
"def backward_euler(rhs, f0, tend, dt):\n",
" ''' \n",
" Computes the backward euler method \n",
" :param rhs: a rhs function\n",
" '''\n",
" nsteps = int(tend/dt)\n",
" f = np.zeros(nsteps)\n",
" f[0] = f0\n",
" time = np.linspace(0,tend,nsteps)\n",
" for n in np.arange(nsteps-1):\n",
" fn = f[n]\n",
" tnp1 = time[n+1]\n",
" fnew = fsolve(be_residual, fn, (rhs, fn, dt, tnp1))\n",
" f[n+1] = fnew\n",
" return time, f"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[]"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"f0 = 1.0\n",
"tend = 0.1\n",
"dt = 0.01\n",
"\n",
"t,fbe = backward_euler(rhs_sharp_transient, f0, tend, dt)\n",
"plt.plot(t,fbe, 'r*-',label='Backward Euler')\n",
"\n",
"texact = np.linspace(0,tend,400)\n",
"plt.plot(texact,yexact(a,f0,texact),label='Exact')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Using Crank-Nicolson"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [],
"source": [
"def cn_residual(fnp1, rhs, fn, dt, tnp1, tn):\n",
" '''\n",
" Nonlinear residual function for the Crank-Nicolson implicit time integrator\n",
" '''\n",
" return fnp1 - fn - 0.5 * dt * ( rhs(fnp1, tnp1) + rhs(fn, tn) )\n",
"\n",
"def crank_nicolson(rhs,f0,tend,dt):\n",
" nsteps = int(tend/dt)\n",
" f = np.zeros(nsteps)\n",
" f[0] = f0\n",
" time = np.linspace(0,tend,nsteps)\n",
" for n in np.arange(nsteps-1):\n",
" fn = f[n]\n",
" tnp1 = time[n+1]\n",
" tn = time[n]\n",
" fnew = fsolve(cn_residual, fn, (rhs, fn, dt, tnp1, tn))\n",
" f[n+1] = fnew\n",
" return time, f"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[]"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"f0 = 1.0\n",
"tend = 0.1\n",
"dt = 0.0025\n",
"\n",
"t,fcn = crank_nicolson(rhs_sharp_transient, f0, tend, dt)\n",
"plt.plot(t,fcn, 'r*-',label='Backward Euler')\n",
"\n",
"texact = np.linspace(0,tend,400)\n",
"plt.plot(texact,yexact(a,f0,texact),label='Exact')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Using odeint"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"time = np.linspace(0,0.1,10)\n",
"sol = odeint(rhs_sharp_transient,f0,time)\n",
"texact = np.linspace(0,tend,400)\n",
"plt.plot(texact,yexact(a,f0,texact),label='Exact')\n",
"plt.plot(time,sol, 'k-o', label='ODEint')\n",
"plt.legend()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"anaconda-cloud": {},
"hide_input": false,
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.3"
},
"toc": {
"colors": {
"hover_highlight": "#DAA520",
"navigate_num": "#000000",
"navigate_text": "#333333",
"running_highlight": "#FF0000",
"selected_highlight": "#FFD700",
"sidebar_border": "#EEEEEE",
"wrapper_background": "#FFFFFF"
},
"moveMenuLeft": true,
"nav_menu": {
"height": "30px",
"width": "252px"
},
"navigate_menu": true,
"number_sections": true,
"sideBar": true,
"threshold": 4,
"toc_cell": false,
"toc_section_display": "block",
"toc_window_display": false,
"widenNotebook": false
}
},
"nbformat": 4,
"nbformat_minor": 2
}