{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Instability of multistep method"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Conside the ODE\n",
"$$\n",
"y' = -y\n",
"$$\n",
"with initial condition\n",
"$$\n",
"y(0) = 1\n",
"$$\n",
"The exact solution is\n",
"$$\n",
"y(t) = \\exp(-t)\n",
"$$\n",
"Let us use a 2-step, second order method\n",
"$$\n",
"y_{n+2} - 2.01 y_{n+1} + 1.01 y_n = h[0.995 f_{n+1} - 1.005 f_n]\n",
"$$\n",
"We will use the initial conditions\n",
"$$\n",
"y_0 =1, \\qquad y_1 = \\exp(-h)\n",
"$$\n",
"Then solving for $y_{n+2}$\n",
"$$\n",
"y_{n+2} = 2.01 y_{n+1} - 1.01 y_n + h[0.995 f_{n+1} - 1.005 f_n]\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 5,
"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": [
"Exact solution"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"def yexact(t):\n",
" return np.exp(-t)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The next function implements the 2-step method."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"def f(t,y):\n",
" return -y\n",
"\n",
"def solve(t0,T,y0,h):\n",
" N = int((T-t0)/h)\n",
" y = np.zeros(N)\n",
" t = np.zeros(N)\n",
" t[0], y[0] = t0, y0\n",
" t[1], y[1] = t0+h, np.exp(-(t0+h))\n",
" for n in range(2,N):\n",
" t[n] = t[n-1] + h\n",
" y[n] = 2.01*y[n-1] - 1.01*y[n-2] \\\n",
" + h*(0.995*f(t[n-1],y[n-1]) - 1.005*f(t[n-2],y[n-2]))\n",
" return t, y"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Solve for decreasing h"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"T = 15.0\n",
"t0, y0 = 0.0, 1.0\n",
"H = [1.0/10.0,1.0/20.0,1.0/40]\n",
"\n",
"te = np.linspace(t0,T,100)\n",
"ye = yexact(te)\n",
"plt.figure(figsize=(5,3))\n",
"plt.plot(te,ye,'--')\n",
"\n",
"for h in H:\n",
" t,y = solve(t0,T,y0,h)\n",
" plt.plot(t,y)\n",
"plt.legend(('Exact','h=1/10','h=1/20','h=1/40'))\n",
"plt.xlabel('t')\n",
"plt.ylabel('y')\n",
"plt.axis([0,15,-0.1,1])\n",
"plt.grid(True);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The solution becomes worse with decreasing values of $h$. The scheme has the polynomials\n",
"$$\n",
"\\rho(w) = w^2 - 2.01 w + 1.01, \\qquad \\sigma = 0.995 w - 1.005\n",
"$$\n",
"The roots of $\\rho$ are $1$ and $1.01$. For the homogeneous problem, the 2-step scheme has the solution\n",
"$$\n",
"y_n = A + B (1.01)^n\n",
"$$\n",
"The second root is responsible for existence of growing solutions in the above scheme. As we reduce $h$ we take more steps to reach the final time, which leads to larger growth of the spurious solution."
]
}
],
"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
}