{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Trapezoidal method: Error convergence"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Conside the ODE\n",
"$$\n",
"y' = -y + 2 \\exp(-t) \\cos(2t)\n",
"$$\n",
"with initial condition\n",
"$$\n",
"y(0) = 0\n",
"$$\n",
"The exact solution is\n",
"$$\n",
"y(t) = \\exp(-t) \\sin(2t)\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": true
},
"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": 2,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def yexact(t):\n",
" return np.exp(-t)*np.sin(2.0*t)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This implements Trapezoidal method\n",
"$$\n",
"y_n = y_{n-1} + \\frac{h}{2}[ f(t_{n-1},y_{n-1}) + f(t_n,y_n)]\n",
"$$\n",
"For the present example we get\n",
"$$\n",
"y_n = y_{n-1} + \\frac{h}{2}[ -y_{n-1} + 2 \\exp(-t_{n-1}) \\cos(2t_{n-1}) - y_n + 2 \\exp(-t_n) \\cos(2t_n) ]\n",
"$$\n",
"Solving for $y_n$\n",
"$$\n",
"y_n = \\frac{1}{1 + \\frac{h}{2}} \\left\\{ (1 - \\frac{h}{2}) y_{n-1} + h [\\exp(-t_{n-1}) \\cos(2t_{n-1}) + \\exp(-t_n) \\cos(2t_n) ] \\right\\}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def trap(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",
" t[n] = t[n-1] + h\n",
" y[n] = (1.0-0.5*h)*y[n-1] + \\\n",
" h*(np.exp(-t[n-1])*np.cos(2.0*t[n-1]) + np.exp(-t[n])*np.cos(2.0*t[n]))\n",
" y[n] = y[n]/(1.0 + 0.5*h)\n",
" return t, y"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Solve for fixed h"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"t0,y0 = 0.0,0.0\n",
"T = 10.0\n",
"h = 1.0/20.0\n",
"t,y = trap(t0,T,y0,h)\n",
"te = np.linspace(t0,T,100)\n",
"ye = yexact(te)\n",
"plt.plot(t,y,te,ye,'--')\n",
"plt.legend(('Numerical','Exact'))\n",
"plt.xlabel('t')\n",
"plt.ylabel('y')\n",
"plt.title('Step size = ' + str(h));"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Solve for several h"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Study the effect of decreasing step size. The error is plotted in log scale."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"5.000000e-02 5.573877e-04 2.001287e+00\n",
"2.500000e-02 1.393158e-04 2.000322e+00\n",
"1.250000e-02 3.482702e-05 2.000080e+00\n",
"6.250000e-03 8.706634e-06 2.000020e+00\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"hh = 0.1/2.0**np.arange(5)\n",
"err=np.zeros(len(hh))\n",
"for (i,h) in enumerate(hh):\n",
" t,y = trap(t0,T,y0,h)\n",
" ye = yexact(t)\n",
" err[i] = np.abs(y-ye).max()\n",
" plt.semilogy(t,np.abs(y-ye))\n",
" \n",
"plt.legend(hh)\n",
"plt.xlabel('t')\n",
"plt.ylabel('log(error)')\n",
"\n",
"for i in range(1,len(hh)):\n",
" print('%e %e %e'%(hh[i],err[i],np.log(err[i-1]/err[i])/np.log(2.0)))\n",
"\n",
"# plot error vs h\n",
"plt.figure()\n",
"plt.loglog(hh,err,'o-')\n",
"plt.xlabel('h')\n",
"plt.ylabel('Error norm');"
]
}
],
"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
}