{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Instability of forward Euler scheme"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Consider the ODE\n",
"$$\n",
"y' = -5 t y^2 + \\frac{5}{t} - \\frac{1}{t^2}, \\qquad t \\ge 1\n",
"$$\n",
"with initial condition\n",
"$$\n",
"y(1) = 1\n",
"$$\n",
"The exact solution is\n",
"$$\n",
"y(t) = \\frac{1}{t}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 1,
"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": [
"Right hand side function"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"def f(t,y):\n",
" return -5.0*t*y**2 + 5.0/t - 1.0/t**2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Exact solution"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"def yexact(t):\n",
" return 1.0/t"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This implements Euler method\n",
"$$\n",
"y_n = y_{n-1} + h f(t_{n-1},y_{n-1})\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"def euler(t0,T,y0,h):\n",
" N = int((T-t0)/h)\n",
" y = np.zeros(N)\n",
" t = np.zeros(N)\n",
" y[0] = y0\n",
" t[0] = t0\n",
" for n in range(1,N):\n",
" y[n] = y[n-1] + h*f(t[n-1],y[n-1])\n",
" t[n] = t[n-1] + h\n",
" return t, y"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Number of iterations= 126\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
"