{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Stiff ODE: Forward Euler with variable step"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Consider the stiff ODE\n",
"$$\n",
"y' = -100(y - \\sin(t)), \\qquad t \\ge 0\n",
"$$\n",
"$$\n",
"y(0) = 1\n",
"$$\n",
"The exact solution is\n",
"$$\n",
"y(t) = \\frac{10101}{10001}e^{-100t} - \\frac{100}{10001}(\\cos t - 100 \\sin t)\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 42,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"%config InlineBackend.figure_format = 'svg'\n",
"import numpy as np\n",
"from matplotlib import pyplot as plt"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The next two functions implement the right hand side of the ode and the exact solution."
]
},
{
"cell_type": "code",
"execution_count": 43,
"metadata": {},
"outputs": [],
"source": [
"def f(t,y):\n",
" return -100.0*(y - np.sin(t))\n",
"\n",
"def yexact(t):\n",
" return 10101.0*np.exp(-100*t)/10001.0 - 100.0*(np.cos(t) - 100.0*np.sin(t))/10001.0"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Forward Euler"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"y_n = y_{n-1} - 100 h [ y_{n-1} - \\sin(t_{n-1})]\n",
"$$\n",
"We start with a step size of $h=0.01$ until $t=0.1$ after which we increase it to $h=0.021$ which is slightly above the stability limit of 0.02."
]
},
{
"cell_type": "code",
"execution_count": 44,
"metadata": {},
"outputs": [],
"source": [
"def ForwardEuler(t0,y0,T,hlo,hhi):\n",
" t, y = [], []\n",
" t.append(t0)\n",
" y.append(y0)\n",
" h = hlo\n",
" n = 1\n",
" while t[n-1]+h <= T:\n",
" if t[n-1] > 0.1:\n",
" h = hhi\n",
" y.append(y[n-1] + h*f(t[n-1],y[n-1]))\n",
" t.append(t[n-1] + h)\n",
" n += 1\n",
" t, y = np.array(t),np.array(y)\n",
" plt.plot(t,y,t,yexact(t))\n",
" plt.title(\"Forward Euler\")\n",
" plt.grid(True)\n",
" plt.legend((\"Euler\",\"Exact\"));\n",
" return t, y"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Use a step size slightly below the stability limit of 0.02."
]
},
{
"cell_type": "code",
"execution_count": 45,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"t0, y0 = 0.0, 1.0\n",
"T = 2*np.pi\n",
"t,y = ForwardEuler(t0,y0,T,0.019,0.019)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The solution is oscillatory initially. Use a smaller step size of 0.01"
]
},
{
"cell_type": "code",
"execution_count": 46,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"t,y = ForwardEuler(t0,y0,T,0.01,0.01)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The exponential decay happens quite fast, so we may think of using a small step initially and then increase it once the decay is completed."
]
},
{
"cell_type": "code",
"execution_count": 47,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"t,y = ForwardEuler(t0,y0,T,0.01,0.021)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The exponential mode might have decayed but we have to still use a small step size for the entire simulation, since this mode can grow back and lead to instability."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"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.6.6"
}
},
"nbformat": 4,
"nbformat_minor": 1
}