{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Stiff ODE" ] }, { "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": 11, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from matplotlib import pyplot as plt" ] }, { "cell_type": "code", "execution_count": 12, "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", "$$" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "def ForwardEuler(t0,y0,T,h):\n", " N = int((T-t0)/h)\n", " y = np.zeros(N)\n", " t = np.zeros(N)\n", " t[0] = t0\n", " y[0] = y0\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": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/svg+xml": [ "\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "t0 = 0\n", "y0 = 1\n", "T = np.pi\n", "\n", "h = 1.0/52.0\n", "t,y = ForwardEuler(t0,y0,T,h)\n", "\n", "plt.plot(t,y,t,yexact(t))\n", "plt.xlim([-0.1,4])\n", "plt.title(\"Forward Euler, h =\"+str(h))\n", "plt.legend((\"Euler\",\"Exact\"))" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/svg+xml": [ "\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "h = 0.01\n", "t,y = ForwardEuler(t0,y0,T,h)\n", "\n", "plt.plot(t,y,t,yexact(t))\n", "plt.xlim([-0.1,4])\n", "plt.title(\"Forward Euler, h =\"+str(h))\n", "plt.legend((\"Euler\",\"Exact\"))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Backward Euler scheme" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "y_n = y_{n-1} + h [ -100(y_n - \\sin(t_n) ]\n", "$$\n", "or\n", "$$\n", "y_n = \\frac{ y_{n-1} + 100 h \\sin(t_n)}{1 + 100 h}\n", "$$" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "def BackwardEuler(t0,y0,T,h):\n", " N = int((T-t0)/h)\n", " y = np.zeros(N)\n", " t = np.zeros(N)\n", " t[0] = t0\n", " y[0] = y0\n", " for n in range(1,N):\n", " t[n] = t[n-1] + h\n", " y[n] = (y[n-1] + 100.0*h*np.sin(t[n]))/(1.0 + 100.0*h)\n", " return t,y" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Trapezoidal scheme" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "y_n = y_{n-1} + \\frac{h}{2}[ -100(y_{n-1}-\\sin(t_{n-1})) - 100 (y_n - \\sin(t_n))]\n", "$$\n", "or\n", "$$\n", "y_n = \\frac{(1-50h)y_{n-1} + 50h(\\sin(t_{n-1}) + \\sin(t_n))}{1 + 50 h}\n", "$$" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "def Trapezoidal(t0,y0,T,h):\n", " N = int((T-t0)/h)\n", " y = np.zeros(N)\n", " t = np.zeros(N)\n", " t[0] = t0\n", " y[0] = y0\n", " for n in range(1,N):\n", " t[n] = t[n-1] + h\n", " y[n] = ((1.0-50.0*h)*y[n-1] + 50.0*h*(np.sin(t[n-1])+np.sin(t[n])))/(1.0 + 50.0*h)\n", " return t,y" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/svg+xml": [ "\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "h = 0.1\n", "t,y = BackwardEuler(t0,y0,T,h)\n", "plt.plot(t,y)\n", "\n", "t,y = Trapezoidal(t0,y0,T,h)\n", "plt.plot(t,y)\n", "\n", "plt.plot(t,yexact(t))\n", "plt.xlim([-0.1,4])\n", "plt.title(\"Step size h =\"+str(h))\n", "plt.legend((\"Backward Euler\",\"Trapezoid\",\"Exact\"))" ] } ], "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 }