{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d27068c4",
   "metadata": {},
   "source": [
    "# Adams-Bashforth and Adams-Moulton methods"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1c2b9b77",
   "metadata": {},
   "source": [
    "This notebook illustrates the application of ``SymPy`` to construct Adams-Bashforth and Adams-Moulton methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "cf00bb81",
   "metadata": {},
   "outputs": [],
   "source": [
    "using SymPy"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "44e87b3c",
   "metadata": {},
   "source": [
    "## 1. A 3-point Adams-Bashforth Method"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba1eb651",
   "metadata": {},
   "source": [
    "We declare the parameter h, the weights, and a free variable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "859bd8f9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(h, c1, c2, c3, z)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "h, c1, c2, c3, z = Sym(\"h, c1, c2, c3, z\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d39e05d8",
   "metadata": {},
   "source": [
    "If ``xn`` is ``0``, then ``xnminus1`` must be ``-h``, ``xnminus2`` must be ``-2*h``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "812e3d85",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0, -h, -2*h)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(xn, xnminus1, xnminus2) = (0, -h, -2*h)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "42d9bba1",
   "metadata": {},
   "source": [
    "We can now define the rule symbolically."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "0ad8bda5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "rule (generic function with 1 method)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rule(f) = c1*f(xn) + c2*f(xnminus1) + c3*f(xnminus2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a7346256",
   "metadata": {},
   "source": [
    "The basis functions are $1$, $x$, and $x^2$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "aa4df883",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "b3 (generic function with 1 method)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b1(x) = Sym(\"1\")\n",
    "b2(x) = x\n",
    "b3(x) = x^2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2b74794d",
   "metadata": {},
   "source": [
    "The right hand sides of the equations express that the rule must integrate all basis functions correctly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "983d0447",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\frac{h^{3}}{3}$"
      ],
      "text/plain": [
       " 3\n",
       "h \n",
       "--\n",
       "3 "
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rhs1 = SymPy.integrate(b1(z), (z, xn, xn+h))\n",
    "rhs2 = SymPy.integrate(b2(z), (z, xn, xn+h))\n",
    "rhs3 = SymPy.integrate(b3(z), (z, xn, xn+h))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d913a4fd",
   "metadata": {},
   "source": [
    "The rule must integrate the basis functions correctly:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "00fd4a56",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The equations :\n",
      "1*c1 + 1*c2 + 1*c3 - 1*h\n",
      "-c2*h - 2*c3*h - h^2/2\n",
      "c2*h^2 + 4*c3*h^2 - h^3/3\n"
     ]
    }
   ],
   "source": [
    "eq1 = rule(b1) - rhs1\n",
    "eq2 = rule(b2) - rhs2\n",
    "eq3 = rule(b3) - rhs3\n",
    "println(\"The equations :\")\n",
    "println(eq1)\n",
    "println(eq2)\n",
    "println(eq3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "d930ac1b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Dict{Any, Any} with 3 entries:\n",
       "  c1 => 23*h/12\n",
       "  c2 => -4*h/3\n",
       "  c3 => 5*h/12"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sol = SymPy.solve([eq1,eq2,eq3],[c1, c2, c3])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "db7551c3",
   "metadata": {},
   "source": [
    "Thus we obtain the following method:\n",
    "\n",
    "$$\n",
    "    y_{n+1} = y_n + \\frac{h}{12} \n",
    "    \\left( \\vphantom{\\frac{h}{12}} 23 f_n - 16 f_{n-1} + 5 f_{n-2} \\right).\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6380e57e",
   "metadata": {},
   "source": [
    "## 2. A 3-point Adams-Moulton Method"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc9da9db",
   "metadata": {},
   "source": [
    "We proceed in the same fashion to derive a 3-point Adams-Moulton method.  The main difference is in the definition of the ``xnplus1`` and the ``xnminus1``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "23d1ad90",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The equations :\n",
      "1*c1 + 1*c2 + 1*c3 - 1*h\n",
      "c1*h - c3*h - h^2/2\n",
      "c1*h^2 + c3*h^2 - h^3/3\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Dict{Any, Any} with 3 entries:\n",
       "  c1 => 5*h/12\n",
       "  c2 => 2*h/3\n",
       "  c3 => -h/12"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# declare the parameter h, the weights, and a free variable\n",
    "h, c1, c2, c3, z = Sym(\"h, c1, c2, c3, z\")\n",
    "# if xn is 0, then xnplus1 must be +h, xnminus1 must be -h\n",
    "(xnplus1, xn, xnminus1) = (h, 0, -h)\n",
    "# we can now define the rule symbolically\n",
    "rule(f) = c1*f(xnplus1) + c2*f(xn) + c3*f(xnminus1)\n",
    "# the basis functions are\n",
    "b1(x) = Sym(\"1\")\n",
    "b2(x) = x\n",
    "b3(x) = x^2\n",
    "# the right hand sides of the equations\n",
    "rhs1 = SymPy.integrate(b1(z), (z, xn, xn+h))\n",
    "rhs2 = SymPy.integrate(b2(z), (z, xn, xn+h))\n",
    "rhs3 = SymPy.integrate(b3(z), (z, xn, xn+h))\n",
    "# the rule must integrate the basis functions correctly\n",
    "eq1 = rule(b1) - rhs1\n",
    "eq2 = rule(b2) - rhs2\n",
    "eq3 = rule(b3) - rhs3\n",
    "println(\"The equations :\")\n",
    "println(eq1)\n",
    "println(eq2)\n",
    "println(eq3)\n",
    "sol = SymPy.solve([eq1,eq2,eq3],[c1, c2, c3])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "47238c1b",
   "metadata": {},
   "source": [
    "Thus we obtain the following method:\n",
    "\n",
    "$$\n",
    "    y_{n+1} = y_n + \\frac{h}{12} \n",
    "    \\left( \\vphantom{\\frac{h}{12}} 5 f_{n+1} + 8 f_n - f_{n-1} \\right).\n",
    "$$"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.8.0",
   "language": "julia",
   "name": "julia-1.8"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.8.0"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}