{ "cells": [ { "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": {}, "outputs": [], "source": [ "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": {}, "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": {}, "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": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/svg+xml": [ "\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "t0 = 0\n", "T = 10\n", "h = 1.0/20.0\n", "t,y = trap(t0,T,0,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": [ "Study the effect of decreasing step size. The error is plotted in log scale." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "hh = [1.0/2.0, 1.0/10.0, 1.0/50.0]\n", "for h in hh:\n", " t,y = trap(t0,T,0,h)\n", " ye = yexact(t)\n", " plt.semilogy(t,np.abs(y-ye))\n", " plt.legend(hh)\n", " plt.xlabel('t')\n", " plt.ylabel('log(error)')" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 2", "language": "python", "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.13" } }, "nbformat": 4, "nbformat_minor": 1 }