{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "50e88f84-b48e-4671-be52-c1feae47d0ca",
   "metadata": {},
   "source": [
    "(module:NumericalMethods)=\n",
    "# Module `NumericalMethods`\n",
    "\n",
    "**Last revised on August 25, 2025**."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e9ae6a31-817a-4bb9-a623-342bf3016526",
   "metadata": {},
   "source": [
    "The following code cells are the content of file `NumericalMethods.jl`, used to define the module `NumericalMethods`\n",
    "\n",
    "This file exists for two reasons:\n",
    "\n",
    "1. It can be convenient to gather the cells defining various functions for the module in a notebook like this one,\n",
    "which allows documentation, and then convert to the the \".jl\" file with the JupterLab command\n",
    "<br>`File > Save and Export Notebook As ... > Executable Script`\n",
    "This gathers the contents of the code cells, ignorig any markdown cells.\n",
    "\n",
    "2. This description of the module's definitions can be used as a section in a Jupyter Book."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d6e193db-e778-44cc-869f-e7f45cb98199",
   "metadata": {},
   "source": [
    "Usage is:\n",
    "\n",
    "    include(\"NumericalMethods.jl\")\n",
    "then\n",
    "\n",
    "    using .NumericalMethods\n",
    "\n",
    "or for a particular function, like\n",
    "\n",
    "    import .NumericalMethods: newtonmethod"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "c7faa07e-4350-41b9-a74d-25480a1868ca",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Module `NumericalMethods`\n",
    "\n",
    "# Version of 2022-11-13\n",
    "\n",
    "# The following code cells are the content of file `NumericalMethods.jl`, used to define the module `NumericalMethods`\n",
    "\n",
    "# The notebook file version exists for two reasons:\n",
    "#\n",
    "# 1. It can be convenient to gather the cells defining various functions for the module in a notebook like this one,\n",
    "#    which allows documentation, and then convert to the the \".jl\" file with the JupterLab command\n",
    "#    `File > Save and Export Notebook As ... > Executable Script`\n",
    "#\n",
    "# This gathers the contents of the code cells, ignorig any markdown cells.\n",
    "\n",
    "# 2. This description of the module's definitions can be used as a section in a Jupyter Book."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "d521b961-5da6-4346-908e-0f14c7c2c101",
   "metadata": {},
   "outputs": [
    {
     "ename": "LoadError",
     "evalue": "ParseError:\n\u001b[90m# Error @ \u001b[0;0m\u001b]8;;file:///Users/brenton/Library/CloudStorage/OneDrive-CollegeofCharleston/numerical-methods-and-analysis/NMAA-julia/docs/In[2]#1:24\u001b\\\u001b[90mIn[2]:1:24\u001b[0;0m\u001b]8;;\u001b\\\nmodule NumericalMethods\u001b[48;2;120;70;70m\u001b[0;0m\n\u001b[90m#                      └ ── \u001b[0;0m\u001b[91mpremature end of input\u001b[0;0m",
     "output_type": "error",
     "traceback": [
      "ParseError:\n\u001b[90m# Error @ \u001b[0;0m\u001b]8;;file:///Users/brenton/Library/CloudStorage/OneDrive-CollegeofCharleston/numerical-methods-and-analysis/NMAA-julia/docs/In[2]#1:24\u001b\\\u001b[90mIn[2]:1:24\u001b[0;0m\u001b]8;;\u001b\\\nmodule NumericalMethods\u001b[48;2;120;70;70m\u001b[0;0m\n\u001b[90m#                      └ ── \u001b[0;0m\u001b[91mpremature end of input\u001b[0;0m",
      "",
      "Stacktrace:",
      " [1] top-level scope",
      "   @ In[2]:1"
     ]
    }
   ],
   "source": [
    "module NumericalMethods"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1da5e83f-16be-4bca-ace5-b3d0dc838212",
   "metadata": {
    "tags": []
   },
   "source": [
    "(root-finding)=\n",
    "## Root-finding"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c4214f1-a7aa-4177-9fee-871a2773c495",
   "metadata": {},
   "source": [
    "(newtons-method)=\n",
    "### Newton's method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "94d44a69-875e-4b3e-9870-2ada4f52cb19",
   "metadata": {},
   "outputs": [],
   "source": [
    "function newton_method(f, Df, x0, errortolerance; maxiterations=20, demomode=false)\n",
    "    # Basic usage is:\n",
    "    # (rootApproximation, errorEstimate, iterations) = newton(f, Df, x0, errorTolerance)\n",
    "    # There is an optional input parameter \"demomode\" which controls whether to\n",
    "    # - println intermediate results (for \"study\" purposes), or to\n",
    "    # - work silently (for \"production\" use).\n",
    "    # The default is silence.\n",
    "\n",
    "    if demomode\n",
    "        println(\"Solving by Newton's Method.\")\n",
    "        println(\"maxiterations = $maxiterations\")\n",
    "        println(\"errortolerance = $errortolerance\")\n",
    "    end\n",
    "    x = x0\n",
    "    global errorestimate  # make it global to this function; without this it would be local to the \"for\" loop.\n",
    "    for iteration in 1:maxiterations\n",
    "        fx = f(x)\n",
    "        Dfx = Df(x)\n",
    "        # Note: a careful, robust code would check for the possibility of division by zero here,\n",
    "        # but for now I just want a simple presentation of the basic mathematical idea.\n",
    "        dx = fx/Dfx\n",
    "        x -= dx  # Aside: this is shorthand for \"x = x - dx\"\n",
    "        errorestimate = abs(dx);\n",
    "        if demomode\n",
    "            println(\"At iteration $iteration, x = $x with estimated error $errorestimate and backward error $(abs(f(x)))\")\n",
    "        end\n",
    "        if errorestimate <= errortolerance\n",
    "            if demomode\n",
    "                println(\"Done!\")\n",
    "            end\n",
    "            return (x, errorestimate, iteration)\n",
    "        end\n",
    "    end\n",
    "    # Note: if we get to here (no \"return\" above), it completed maxIterations iterations without satisfying the accuracy target,\n",
    "    # but we still return the information that we have.\n",
    "    return (x, errorestimate, maxiterations)\n",
    "end;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "43c5912e-24f3-4504-a63a-5e2348c83cc5",
   "metadata": {},
   "source": [
    "(secant-method)=\n",
    "### The secant method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "6966d496-41ee-459e-9989-b290e5ac6719",
   "metadata": {},
   "outputs": [],
   "source": [
    "function secant_method(f, a, b, errortolerance; maxiterations=20, demomode=false)\n",
    "    # Solve f(x)=0 in the interval [a, b] by the Secant Method.\n",
    "    if demomode\n",
    "        print(\"Solving by the Secant Method.\")\n",
    "    end;\n",
    "    # Some more descriptive names\n",
    "    x_older = a\n",
    "    x_more_recent = b\n",
    "    f_x_older = f(x_older)\n",
    "    f_x_more_recent = f(x_more_recent)\n",
    "    for iteration in 1:maxiterations\n",
    "        global x_new, errorestimate\n",
    "        if demomode\n",
    "            println(\"\\nIteration $(iteration):\")\n",
    "        end;\n",
    "        x_new = (x_older * f_x_more_recent - f_x_older * x_more_recent)/(f_x_more_recent - f_x_older)\n",
    "        f_x_new = f(x_new)\n",
    "        (x_older, x_more_recent) = (x_more_recent, x_new)\n",
    "        (f_x_older, f_x_more_recent) = (f_x_more_recent, f_x_new)\n",
    "        errorestimate = abs(x_older - x_more_recent)\n",
    "        if demomode\n",
    "            println(\"The latest pair of approximations are $x_older and $x_more_recent,\")\n",
    "            println(\"where the function's values are $f_x_older and $f_x_more_recent respectively.\")\n",
    "            println(\"The new approximation is $x_new with estimated error $errorestimate and backward error $(abs(f_x_new))\")\n",
    "        end;\n",
    "        if errorestimate < errortolerance\n",
    "            break\n",
    "        end;\n",
    "    end;\n",
    "    # Whether we got here due to accuracy of running out of iterations,\n",
    "    # return the information we have, including an error estimate:\n",
    "    return (x_new,  errorestimate)\n",
    "end;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "543a2ac2-93eb-4585-9b69-943e6edeec04",
   "metadata": {},
   "source": [
    "(linear-algebra)=\n",
    "## Linear Algebra and Simultaneous Equations"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bb6fbff6-cd58-4533-a3d1-6796864e1eb6",
   "metadata": {},
   "source": [
    "(row-reduction)=\n",
    "### Row Reduction\n",
    "(with no pivoting)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "daac9312-13a1-43e1-92eb-48fbd90ee060",
   "metadata": {},
   "outputs": [],
   "source": [
    "function row_reduce(A, b)\n",
    "    # To avoid modifying the matrix and vector specified as input,\n",
    "    # they are copied to new arrays, with the function copy().\n",
    "    # Warning: it does not work to say \"U = A\" and \"c = b\";\n",
    "    # this makes these names synonyms, referring to the same stored data.\n",
    "\n",
    "    U = copy(A)  # not \"U=A\", which makes U and A synonyms\n",
    "    c = copy(b)\n",
    "    n = length(b)\n",
    "    L = zeros(n, n)\n",
    "    for k in 1:n-1\n",
    "        for i in k+1:n\n",
    "            # compute all the L values for column k:\n",
    "            L[i,k] = U[i,k] / U[k,k]  # Beware the case where U[k,k] is 0\n",
    "            for j in k+1:n\n",
    "                U[i,j] -= L[i,k] * U[k,j]\n",
    "            end\n",
    "            # Put in the zeros below the main diagonal in column k of U;\n",
    "            # this is not important for calculations, since those elements of U are not used in backward substitution,\n",
    "            # but it helps for displaying results and for checking the results via residuals.\n",
    "            U[i,k] = 0.\n",
    "            \n",
    "            c[i] -= L[i,k] * c[k]\n",
    "        end\n",
    "    end\n",
    "    for i in 2:n\n",
    "        for j in 1:i-1\n",
    "            U[i,j] = 0.\n",
    "        end\n",
    "    end\n",
    "    return (U, c)\n",
    "end;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "927eb85e-a195-42a6-b84a-73f2987fd227",
   "metadata": {},
   "source": [
    "(backward_substitution)=\n",
    "### Backward substitution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "1137af9a-a821-4ba3-a749-b5e9a70ed2a8",
   "metadata": {},
   "outputs": [],
   "source": [
    "function backward_substitution(U, c; demomode=false)\n",
    "    n = length(c)\n",
    "    x = zeros(n)\n",
    "    x[end] = c[end]/U[end,end]\n",
    "    if demomode\n",
    "        println(\"x_$n = $(x[n])\")\n",
    "    end\n",
    "    for i in n-1:-1:1\n",
    "        if demomode\n",
    "            println(\"i=$i\")\n",
    "        end\n",
    "        x[i] = ( c[i] - sum(U[i,i+1:end] .* x[i+1:end]) ) / U[i,i]\n",
    "        if demomode\n",
    "            print(\"x_$i = $(x[i])\")\n",
    "        end\n",
    "    end\n",
    "    return x\n",
    "end;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4a01bc64-48e2-4c0a-8374-b8559cf0fb5b",
   "metadata": {},
   "source": [
    "(solve-linear-system)=\n",
    "### Solve a linear system (no pivoting)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "5b0041db-9f9d-4fd1-9318-51da08ec6043",
   "metadata": {},
   "outputs": [],
   "source": [
    "solvelinearsystem(A, b) = backwardsubstitution(rowreduce(A, b)...);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f20f6502-4c5c-4f8a-ae41-18c053e520ef",
   "metadata": {
    "tags": []
   },
   "source": [
    "(lu-factorization)=\n",
    "### LU factorization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "0ff3daf4-90cf-49e0-86a1-4bba4c75e8ff",
   "metadata": {},
   "outputs": [],
   "source": [
    "function lu_factorize(A; demomode=false)\n",
    "    # Compute the Doolittle LU factorization of A.\n",
    "    # Sums like $\\sum_{s=1}^{k-1} l_{k,s} u_{s,j}$ are done as matrix products;\n",
    "    # in the above case, row matrix L[k, 1:k-1] by column matrix U[1:k-1,j] gives the sum for a give j,\n",
    "    # and row matrix L[k, 1:k-1] by matrix U[1:k-1,k:n] gives the relevant row vector.\n",
    "    n = size(A)[1]  # First component of the array's size; size(A) returns \"(rows, columns)\"\n",
    "    # Initialize U as a zero matrix;\n",
    "    # correct below the main diagonal, with the other entries to be computed and filled below.\n",
    "    U = zeros(n,n)\n",
    "    # Initialize L as a zero matrix;\n",
    "    # correct above the main diagonal, with the other entries to be computed and filled in below.\n",
    "    L = zeros(n,n)\n",
    "    # Column and row 1 are special:\n",
    "    U[1,:] = A[1,:]\n",
    "    L[1,1] = 1.\n",
    "    L[2:end,1] = A[2:end,1]/U[1,1]\n",
    "    if demomode\n",
    "        println(\"After step k=1\")\n",
    "        println(\"U=\"); printmatrix(U)\n",
    "        println(\"L=\"); printmatrix(L)\n",
    "    end;\n",
    "    for k in 2:n-1\n",
    "        # Julia note: it is necessary to use indices \"[k]\" and so on to get a one-row matrix instead of a vector.\n",
    "        U[[k],k:end] = A[[k],k:end] - L[[k],1:k] * U[1:k,k:end]\n",
    "        L[k,k] = 1.\n",
    "        L[k+1:end,k] = ( A[k+1:end,k] - L[k+1:end,1:k] * U[1:k,k] )/U[k,k]\n",
    "        if demomode\n",
    "            println(\"After step k=$k\")\n",
    "            println(\"U=\"); printmatrix(U)\n",
    "            println(\"L=\"); printmatrix(L)\n",
    "        end;\n",
    "    end;\n",
    "    # The last row is also special: nothing to do for L\n",
    "    L[end,end] = 1.\n",
    "    U[end,end] = A[end,end] - sum(L[[n],1:end-1] * U[1:end-1,end])\n",
    "    if demomode\n",
    "        println(\"After step k=$n\")\n",
    "        println(\"U=\"); printmatrix(U)\n",
    "    end;\n",
    "    return [L, U]\n",
    "end;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5ff33272-988f-4207-97d8-b57c7165f401",
   "metadata": {},
   "source": [
    "(forward_substitution)=\n",
    "### Forward substitution\n",
    "(without pivoting)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "19efc78f-410d-49ba-98ad-8d09e7884134",
   "metadata": {},
   "outputs": [],
   "source": [
    "function forward_substitution(L, b)\n",
    "    # Solve L c = b for c.\n",
    "    n = length(b)\n",
    "    c = zeros(n)\n",
    "    c[1] = b[1]\n",
    "    for i in 2:n\n",
    "        c[i] = b[i] - sum(L[i:i,1:i] * c[1:i])\n",
    "    end;\n",
    "    return c\n",
    "end;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a8bb584-29bf-40de-861a-ecdf58a41d61",
   "metadata": {},
   "source": [
    "(plu-factorization)=\n",
    "### PLU factorization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "f6a5a07e-972a-467b-9fb1-eedd66711dc6",
   "metadata": {},
   "outputs": [],
   "source": [
    "function plu(A; demomode=false)\n",
    "    # Compute the Doolittle PA=LU factorization of A —\n",
    "    # but with the permutation recorded as permutation vector, not as the permutation matrix P.\n",
    "    # Sums like $\\sum_{s=1}^{k-1} l_{k,s} u_{s,j}$ are done as matrix products;\n",
    "    # in the above case, row matrix L[k, 1:k-1] by column matrix U[1:k-1,j] gives the sum for a give j,\n",
    "    # and row matrix L[k, 1:k-1] by matrix U[1:k-1,k:n] gives the relevant row vector.\n",
    "\n",
    "    n = size(A)[1]  # gives the number of rows in the 2D array.\n",
    "    π = zeros(Int, n)\n",
    "    # Julia can use Greek letters (and in fact, UNICODE):\n",
    "    # to insert character π, type \\pi, hit tab, and select \"π\" from the menu.\n",
    "    # Or just call it \"perm\" or such.\n",
    "    π = collect(1:n)\n",
    "    # Julia language note: function \"collect\" converts the abstract entity \"1:n\" into an array of numbers.\n",
    "    \n",
    "    # Initialize U as the zero matrix;\n",
    "    # correct below the main diagonal, with the other entries to be computed below.\n",
    "    U = zeros(n,n)\n",
    "\n",
    "    # Initialize L as zeros;\n",
    "    # correct above the main diagonal, with the other entries to be computed below,\n",
    "    # including the ones on the diagonal.\n",
    "    L = zeros(n,n)\n",
    "\n",
    "    for k in 1:n-1\n",
    "        if demomode; println(\"k=$k\"); end\n",
    "        # Find the pivot element in column k:\n",
    "        pivotrow = k\n",
    "        abs_u_ik_max = abs(A[π[k],k])\n",
    "        for row in k+1:n\n",
    "            abs_u_ik = abs(A[π[row],k])\n",
    "            if abs_u_ik > abs_u_ik_max\n",
    "                pivotrow = row\n",
    "                abs_u_ik_max = abs_u_ik\n",
    "            end\n",
    "        end\n",
    "        if pivotrow > k # swap rows, virtually\n",
    "            if demomode; println(\"Swap row $k with row $pivotrow\"); end\n",
    "            (π[k], π[pivotrow]) = (π[pivotrow], π[k])\n",
    "        else\n",
    "            if demomode; println(\"No row swap needed.\"); end\n",
    "        end\n",
    "        U[k,k:end] = A[[π[k]],k:end] - L[[π[k]],1:k] * U[1:k,k:end]\n",
    "        L[π[k],k] = 1.\n",
    "        for row in k+1:n\n",
    "            L[π[row],k] = ( A[π[row],k] - L[π[row],1:k] ⋅ U[1:k,k] ) / U[k,k]\n",
    "            # Julia note: To enter the centered dot notation for the dot product, type \"\\cdot\" and then hit the tab key.\n",
    "        end\n",
    "        if demomode\n",
    "            println(\"permuted A is:\")\n",
    "            for row in 1:n\n",
    "                println(A[π[row],:])\n",
    "            end\n",
    "            println(\"Intermediate L is\"); printmatrix(L)\n",
    "            println(\"Intermediate U is\"); printmatrix(U)\n",
    "        end\n",
    "    end\n",
    "    # The last row (index \"end\") is special: nothing to do for L except put in the 1 on the \"permuted main diagonal\"\n",
    "    L[π[end],end] = 1.\n",
    "    U[end,end] = A[π[end],end] - L[π[end],1:end-1] ⋅ U[1:end-1,end]\n",
    "    if demomode\n",
    "        println(\"After the final step, k=$(n-1)\")\n",
    "        println(\"L is\"); printmatrix(L)\n",
    "        println(\"U is\"); printmatrix(U)\n",
    "    end\n",
    "    return (L, U, π)\n",
    "end;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b6387a0a-1cb0-4039-83e8-3945c61faf3c",
   "metadata": {},
   "source": [
    "(forward-substitution-pivoting)=\n",
    "### Forward substitution with pivoting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "314ee2db-c40c-4e85-8c82-eaa448c448b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "function forward_substitution(L, b, π)\n",
    "    # Version 2: with permutation of rows\n",
    "    # Solve L c = b for c, with permutation of the rows of L and of b.\n",
    "    n = length(b)\n",
    "    c = zeros(n)\n",
    "    c[1] = b[π[1]]\n",
    "    for i in 2:n\n",
    "        c[i] = b[π[i]] - L[π[i], 1:i] ⋅ c[1:i]\n",
    "    end\n",
    "    return c\n",
    "end;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "59824813-8c68-49d3-90cd-ecae3f81b775",
   "metadata": {
    "tags": []
   },
   "source": [
    "(collocation)=\n",
    "## Collocation and Data Fitting"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "023dac37-c2b3-4586-a477-e1b51bffac49",
   "metadata": {},
   "source": [
    "(polynomial-collocation)=\n",
    "### Polynomial collocation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "e99b25b5-5f97-4957-82fe-b11a58fd6e5c",
   "metadata": {},
   "outputs": [],
   "source": [
    "function poly_fit(x, y)\n",
    "    # Version 1: exact collocation.\n",
    "    # Compute the coeffients c_i of the polynomial of lowest degree that collocates the points (x[i], y[i]).\n",
    "    # These are returned in an array c of the same length as x and y, even if the degree is less than the normal length(x)-1,\n",
    "    # in which case the array has some trailing zeroes.\n",
    "    # The polynomial is thus p(x) = c[1] + c[2]x + ... c[d+1] x^d where d =length(x)-1, the nominal degree.\n",
    "    n_nodes = length(x)\n",
    "    degree = n_nodes - 1\n",
    "    V = zeros(n_nodes,n_nodes)\n",
    "    for i in 0:degree\n",
    "        for j in 0:degree\n",
    "             V[i+1,j+1] = x[i+1]^j  # Shift the array indices up by one, since Julia counts from 1, not 0.\n",
    "        end\n",
    "    end\n",
    "    c = solvelinearsystem(V, y)\n",
    "    return c\n",
    "end;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "df9cd722-8328-489d-a96f-dd0d0d5122df",
   "metadata": {},
   "source": [
    "(polyfit-least-squares)=\n",
    "### Least squares polynomial approximation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "0f12922f-0667-4802-81d7-b8102f610ab5",
   "metadata": {},
   "outputs": [],
   "source": [
    "function poly_fit(x, y, n)\n",
    "    # Version 2: least squares fitting.\n",
    "    # Compute the coeffients c_i of the polynomial of degree n that give the best least squares fit to data (x[i], y[i]).\n",
    "    N = length(x)\n",
    "    m = zeros(2n+1)\n",
    "    for k in 0:2n\n",
    "        m[k+1] = sum(x.^k)  # Here and below, shift the indices up by one, since Julia counts from 1, not 0.\n",
    "    end\n",
    "    M = zeros(n+1,n+1)\n",
    "    for i in 0:n\n",
    "        for j in 0:n \n",
    "             M[i+1, j+1] = m[i+j+1]\n",
    "        end\n",
    "    end\n",
    "    p = zeros(n+1)\n",
    "    for k in 0:n\n",
    "        p[k+1] = sum(x.^k .* y)\n",
    "    end\n",
    "    c = solvelinearsystem(M, p)\n",
    "    return c\n",
    "end;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38826586-e528-4779-a1e9-7aa6f6bdba6b",
   "metadata": {},
   "source": [
    "(poly_val)=\n",
    "### Evaluate a polynomial"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "be0d5999-b86c-40bb-9434-5a44ceac901b",
   "metadata": {},
   "outputs": [],
   "source": [
    "function poly_val(x; coeffs)  # coeffs has to be a keyword argument in order that only x gets vectorized\n",
    "    # Evaluate the polynomial with coefficients in c (as given by polyfit, for example).\n",
    "    # If x is an array, the usage becomes y = polyval.(c, x)\n",
    "    # for each element of array x.\n",
    "    y = coeffs[1]\n",
    "    for i in 2:length(coeffs)\n",
    "        y += coeffs[i]*x^(i-1)\n",
    "    end\n",
    "    return y\n",
    "end;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c9a57292-9a94-481a-97bf-3c2db7c649a4",
   "metadata": {},
   "source": [
    "(derivatives-definite-integrals)=\n",
    "## Derivatives and Definite Integrals"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d14b2c23-aa63-4dcb-b58d-f11df7e7081e",
   "metadata": {},
   "source": [
    "(minimization)=\n",
    "## Minimization"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "026a78ce-d251-4f64-9359-59ae5384d89e",
   "metadata": {
    "tags": []
   },
   "source": [
    "(differential-equations)=\n",
    "## Differential Equations"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "325030c1-715a-4012-8ecd-66142d9ddfdd",
   "metadata": {},
   "source": [
    "(euler_method)=\n",
    "### Euler's method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "9e2999b3-720e-4d34-a7fd-fd1fce06696b",
   "metadata": {},
   "outputs": [],
   "source": [
    "function euler_method(f, a, b, u_0; n=100)\n",
    "    # Solve du/dt = f(t, u) for t in [a, b], with initial value u(a) = u_0\n",
    "    h = (b-a)/n\n",
    "    t = range(a, b, n+1)  # Note: \"n\" counts steps, so there are n+1 values for t.\n",
    "    u = zeros(n+1)\n",
    "    u[1] = u_0\n",
    "    for i in 1:n\n",
    "        u[i+1] = u[i] + f(t[i], u[i])*h\n",
    "    end\n",
    "    return (t, u)\n",
    "end;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "953b6a4b-1dcd-49cd-8348-b490750453b4",
   "metadata": {},
   "source": [
    "(euler_method_error_control)="
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "1b75db9d-1ddc-4b60-be87-6160bd5c93c5",
   "metadata": {},
   "outputs": [],
   "source": [
    "function euler_method_error_control(f, a, b, u_0; errortolerance=1e-3, h_min=1e-6, h_max=0.1, steps_max=1000, demomode=false)\n",
    "    steps = 0\n",
    "    t_i = a\n",
    "    U_i = u_0\n",
    "    t = [t_i]\n",
    "    U = [U_i]\n",
    "    h = h_max  # Start optimistically!\n",
    "    while t_i < b && steps < steps_max\n",
    "        K_1 = h*f(t_i, U_i)\n",
    "        K_2 = h*f(t_i + h/2, U_i + K_1/2)\n",
    "        errorestimate = abs(K_1 - K_2)\n",
    "        s = 0.9 * sqrt(errortolerance/errorestimate)\n",
    "        if errorestimate <= errortolerance  # Success!\n",
    "            t_i += h\n",
    "            U_i += K_1\n",
    "            append!(t, t_i)\n",
    "            append!(U, U_i)\n",
    "            # Adjust step size up, but not too big\n",
    "            h = min(s*h, h_max)\n",
    "        else  # Innacurate; reduce step size and try again\n",
    "            h = max(s*h, h_min)\n",
    "            if demomode\n",
    "                println(\"t_i=$t_i: Decreasing step size to $(about(h)) and trying again.\")\n",
    "            end\n",
    "        end\n",
    "        # A refinement not mentioned above; the next step should not overshoot t=b:\n",
    "        if t_i + h > b\n",
    "            h = b - t_i\n",
    "        end\n",
    "        steps += 1\n",
    "    end\n",
    "    return (t, U)\n",
    "    # Note: if the step count ran out, this does not reach t=b, but at least it is correct as far as it goes\n",
    "end;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d100426d-5ead-4ba4-bc15-2877f7e6ab95",
   "metadata": {},
   "source": [
    "(euler_method_system)="
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "bf2cbb40-b320-42f0-bed3-8b56cd68336c",
   "metadata": {},
   "outputs": [],
   "source": [
    "function euler_method_system(f, a, b, u_0, n)\n",
    "    # TO DO: one could use multiple dispatch to keep the name \"eulermethod\".\n",
    "    # The conversion for the system version is mainly \"U[i] -> U[i,:]\"\n",
    "    \n",
    "    h = (b-a)/n\n",
    "    t = range(a, b, n+1)\n",
    "    \n",
    "    # The following three lines and the one in the for loop below change for the system version\n",
    "    n_unknowns = length(u_0)\n",
    "    U = zeros(n+1, n_unknowns)\n",
    "    U[1,:] = u_0  # Only for system version\n",
    "\n",
    "    for i in 1:n\n",
    "        U[i+1,:] = U[i,:] + f(t[i], U[i,:])*h  # For the system version\n",
    "    end\n",
    "    return (t, U)\n",
    "end;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "82bdff13-83cf-47b1-b443-3718259b233a",
   "metadata": {},
   "source": [
    "(explicit_trapezoid)=\n",
    "### The explicit trapezoid method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "202e1871-f3aa-494f-946c-49cc8f24e206",
   "metadata": {},
   "outputs": [],
   "source": [
    "function explicit_trapezoid(f, a, b, u_0; n=100, demomode=false)\n",
    "    # Use the Explict Trapezoid Method (a.k.a Improved Euler) to solve\n",
    "    #     du/dt = f(t, u)\n",
    "    # for t in [a, b], with initial value u(a) = u_0\n",
    "    \n",
    "    h = (b-a)/n\n",
    "    t = range(a, b, n+1)  # Note: \"n\" counts steps, so there are n+1 values for t.\n",
    "    u = zeros(n+1)\n",
    "    u[1] = u_0\n",
    "    for i in 1:n\n",
    "        K_1 = f(t[i], u[i])*h\n",
    "        K_2 = f(t[i]+h, u[i]+K_1)*h\n",
    "        u[i+1] = u[i] + (K_1 + K_2)/2.0\n",
    "    end;\n",
    "    return (t, u)\n",
    "end;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b22a5ee1-36e0-4ef4-b666-71bd1f3df025",
   "metadata": {},
   "source": [
    "(explicit_trapezoid_system)="
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "58072bf9-9121-4563-9261-2fa8b4756280",
   "metadata": {},
   "outputs": [],
   "source": [
    "function explicit_trapezoid_system(f, a, b, u_0, n)\n",
    "    # Use the Explict Trapezoid Method (a.k.a Improved Euler) to solve the system\n",
    "    #    du/dt = f(t, u) for t in [a, b], with initial value u(a) = u_0 \n",
    "    # The conversion for the system version is mainly \"u[i] -> u[i,:]\"\n",
    "\n",
    "    h = (b-a)/n\n",
    "    t = range(a, b, n+1)\n",
    "    n_unknowns = length(u_0)\n",
    "    u = zeros(n+1, n_unknowns)\n",
    "    u[1,:] = u_0\n",
    "    for i in 1:n\n",
    "        K_1 = f(t[i], u[i,:])*h\n",
    "        K_2 = f(t[i]+h, u[i,:]+K_1)*h\n",
    "        u[i+1,:] = u[i,:] + (K_1 + K_2)/2.0\n",
    "    end\n",
    "    return (t, u)\n",
    "end;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "84dcbc71-e804-48bc-8a98-037f7ca9aaac",
   "metadata": {},
   "source": [
    "(explicit_midpoint)=\n",
    "### The explicit midpoint method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "e6535782-5db3-4b43-8197-1f351d13ec7d",
   "metadata": {},
   "outputs": [],
   "source": [
    "function explicit_midpoint(f, a, b, u_0; n=100, demomode=false)\n",
    "    # Use the Explicit Midpoint Method (a.k.a Modified Euler) to solve\n",
    "    #    du/dt = f(t, u) for t in [a, b], with initial value u(a) = u_0\n",
    "\n",
    "    h = (b-a)/n\n",
    "    t = range(a, b, n+1)  # Note: \"n\" counts steps, so there are n+1 values for t.\n",
    "    u = zeros(length(t))\n",
    "    u[1] = u_0\n",
    "    for i in 1:n\n",
    "        K_1 = f(t[i], u[i])*h\n",
    "        K_2 = f(t[i]+h/2, u[i]+K_1/2)*h\n",
    "        u[i+1] = u[i] + K_2\n",
    "    end;\n",
    "    return (t, u)\n",
    "    end;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bde8200d-bceb-4f51-ad4c-65613af90b16",
   "metadata": {},
   "source": [
    "(explicit_midpoint_system)="
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "085fa8a1-9e98-4d7d-91f6-da38a24450a5",
   "metadata": {},
   "outputs": [],
   "source": [
    "function explicit_midpoint_system(f, a, b, u_0, n)\n",
    "    # Use the Explict Midpoint Method (a.k.a Modified Euler) to solve the system\n",
    "    #    du/dt = f(t, u) for t in [a, b], with initial value u(a) = u_0 \n",
    "    # The conversion for the system version is mainly \"u[i] -> u[i,:]\"\n",
    "\n",
    "    h = (b-a)/n\n",
    "    t = range(a, b, n+1)\n",
    "    n_unknowns = length(u_0)\n",
    "    u = zeros(n+1, n_unknowns)\n",
    "    u[1,:] = u_0\n",
    "    for i in 1:n\n",
    "        K_1 = f(t[i], u[i,:])*h\n",
    "        K_2 = f(t[i]+h/2, u[i,:]+K_1/2)*h\n",
    "        u[i+1,:] = u[i,:] + K_2\n",
    "    end\n",
    "    return (t, u)\n",
    "end;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3ef3c1df-f9eb-4df3-9092-4838547165a7",
   "metadata": {},
   "source": [
    "(runge_kutta)=\n",
    "### The Runge-Kutta method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "90f546c5-ac7d-4939-886b-387bfc810c7e",
   "metadata": {},
   "outputs": [],
   "source": [
    "function runge_kutta(f, a, b, u_0; n=100, demomode=false)\n",
    "    # Use the (classical) Runge-Kutta Method to solve\n",
    "    #    du/dt = f(t, u) for t in [a, b], with initial value u(a) = u_0\n",
    "    h = (b-a)/n\n",
    "    t = range(a, b, n+1)  # Note: \"n\" counts steps, so there are n+1 values for t.\n",
    "    u = zeros(length(t))\n",
    "    u[1] = u_0\n",
    "    for i in 1:n\n",
    "        K_1 = f(t[i], u[i])*h\n",
    "        K_2 = f(t[i]+h/2, u[i]+K_1/2)*h\n",
    "        K_3 = f(t[i]+h/2, u[i]+K_2/2)*h\n",
    "        K_4 = f(t[i]+h, u[i]+K_3)*h\n",
    "        u[i+1] = u[i] + (K_1 + 2*K_2 + 2*K_3 + K_4)/6\n",
    "    end;\n",
    "    return (t, u)\n",
    "end;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a0de07e2-963c-4417-a57a-4f88a8cbd8ba",
   "metadata": {},
   "source": [
    "(runge_kutta_system)="
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "9af6e27c-8ff0-442c-bcfe-ac3af6804028",
   "metadata": {},
   "outputs": [],
   "source": [
    "function runge_kutta_system(f, a, b, u_0, n)\n",
    "    # Use the (classical) Runge-Kutta Method to solve\n",
    "    #    du/dt = f(t, u) for t in [a, b], with initial value u(a) = u_0\n",
    "    # The conversion for the system version is mainly \"u[i] -> u[i,:]\"\n",
    "\n",
    "    h = (b-a)/n\n",
    "    t = range(a, b, n+1)\n",
    "    n_unknowns = length(u_0)\n",
    "    u = zeros(n+1, n_unknowns)\n",
    "    u[1,:] = u_0\n",
    "    for i in 1:n\n",
    "        K_1 = f(t[i], u[i,:])*h\n",
    "        K_2 = f(t[i]+h/2, u[i,:]+K_1/2)*h\n",
    "        K_3 = f(t[i]+h/2, u[i,:]+K_2/2)*h\n",
    "        K_4 = f(t[i]+h, u[i,:]+K_3)*h\n",
    "        u[i+1,:] = u[i,:] + (K_1 + 2*K_2 + 2*K_3 + K_4)/6\n",
    "    end\n",
    "    return (t, u)\n",
    "end;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "50dd975f-a658-4eba-8839-1ad61f210282",
   "metadata": {},
   "source": [
    "## Some auxilliary functions\n",
    "\n",
    "For examples, presentation of results, etc."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "58eb6eea-2028-42ce-9bb9-f6cf3615cb52",
   "metadata": {},
   "source": [
    "(print_matrix)=\n",
    "### Helper function `print_matrix`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "45b39442-79ba-469e-bd11-65405d14203c",
   "metadata": {},
   "outputs": [],
   "source": [
    "function print_matrix(A)\n",
    "    # A helper function to \"pretty print\" matrices\n",
    "    (rows, cols) = size(A)\n",
    "    print(\"[ \")\n",
    "    for row in 1:rows\n",
    "        if row > 1\n",
    "            print(\"  \")\n",
    "        end\n",
    "        for col in 1:cols\n",
    "            print(A[row,col], \" \")\n",
    "        end\n",
    "        if row < rows;\n",
    "            println()\n",
    "        else\n",
    "            println(\"]\")\n",
    "        end\n",
    "    end\n",
    "end;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "3195bcda-6eb6-4525-ad45-7583f9d30208",
   "metadata": {},
   "outputs": [],
   "source": [
    "# A helper function for rounding some output to three significant digits\n",
    "approx3(x) = round(x, sigdigits=3);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "09ae2872-b4d5-4875-a7e8-d5429d48edaa",
   "metadata": {},
   "outputs": [],
   "source": [
    "# A helper function for rounding some output to four significant digits\n",
    "approx4(x) = round(x, sigdigits=4);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "69b759f6-4a1c-46f8-ac34-896fe6041ceb",
   "metadata": {},
   "source": [
    "## For some examples in Chapter {doc}`ODE-IVPs`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "2de9649f-99b6-4098-b8db-810c85e3b6c9",
   "metadata": {},
   "outputs": [],
   "source": [
    "f_mass_spring(t, u) = [ u[2], -(K/M)*u[1] - (D/M)*u[2] ];\n",
    "\n",
    "function y_mass_spring(t; t_0, u_0, K, M, D)\n",
    "    (y_0, v_0) = u_0\n",
    "    discriminant = D^2 - 4K*M\n",
    "    if discriminant < 0  # underdamped\n",
    "        omega = sqrt(4K*M - D^2)/(2M)\n",
    "        A = y_0\n",
    "        B = (v_0 + y_0*D/(2M))/omega\n",
    "        return exp(-D/(2M)*(t-t_0)) * ( A*cos(omega*(t-t_0)) + B*sin(omega*(t-t_0)))\n",
    "    elseif discriminant > 0  # overdamped\n",
    "        Delta = sqrt(discriminant)\n",
    "        lambda_plus = (-D + Delta)/(2M)\n",
    "        lambda_minus = (-D - Delta)/(2M)\n",
    "        A = M*(v_0 - lambda_minus * y_0)/Delta\n",
    "        B = y_0 - A\n",
    "        return A*exp(lambda_plus*(t-t_0)) + B*exp(lambda_minus*(t-t_0))\n",
    "    else\n",
    "        lambda = -D/(2M)\n",
    "        A = y_0\n",
    "        B = v_0 - A * lambda\n",
    "        return (A + B*t)*exp(lambda*(t-t_0))\n",
    "    end\n",
    "end;\n",
    "\n",
    "function damping(K, M, D)\n",
    "    if D == 0\n",
    "        println(\"Undamped\")\n",
    "    else\n",
    "        discriminant = D^2 - 4K*M\n",
    "        if discriminant < 0\n",
    "            println(\"Underdamped\")\n",
    "        elseif discriminant > 0\n",
    "            println(\"Overdamped\")\n",
    "        else\n",
    "            println(\"Critically damped\")\n",
    "        end\n",
    "    end\n",
    "end;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "ec71ccdb-f91c-484b-9229-e03d05e581a9",
   "metadata": {},
   "outputs": [
    {
     "ename": "LoadError",
     "evalue": "ParseError:\n\u001b[90m# Error @ \u001b[0;0m\u001b]8;;file:///Users/brenton/Library/CloudStorage/OneDrive-CollegeofCharleston/numerical-methods-and-analysis/NMAA-julia/docs/In[28]#1:1\u001b\\\u001b[90mIn[28]:1:1\u001b[0;0m\u001b]8;;\u001b\\\n\u001b[48;2;120;70;70mend\u001b[0;0m;\n\u001b[90m└─┘ ── \u001b[0;0m\u001b[91minvalid identifier\u001b[0;0m",
     "output_type": "error",
     "traceback": [
      "ParseError:\n\u001b[90m# Error @ \u001b[0;0m\u001b]8;;file:///Users/brenton/Library/CloudStorage/OneDrive-CollegeofCharleston/numerical-methods-and-analysis/NMAA-julia/docs/In[28]#1:1\u001b\\\u001b[90mIn[28]:1:1\u001b[0;0m\u001b]8;;\u001b\\\n\u001b[48;2;120;70;70mend\u001b[0;0m;\n\u001b[90m└─┘ ── \u001b[0;0m\u001b[91minvalid identifier\u001b[0;0m",
      "",
      "Stacktrace:",
      " [1] top-level scope",
      "   @ In[28]:1"
     ]
    }
   ],
   "source": [
    "end;"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.11.4",
   "language": "julia",
   "name": "julia-1.11"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}