{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "10df9104-0ba0-491a-89ae-339ab0eff440",
   "metadata": {},
   "source": [
    "Before you turn this problem in, make sure everything runs as expected.\n",
    "First, restart the kernel (in the menubar, select Kernel → Restart) and then run all cells (in the menubar, select Cell → Run All).\n",
    "\n",
    "Make sure you fill in any place that says `BEGIN SOLUTION` and `END SOLUTION`, as well as your name and collaborators below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a3b03a68-e47f-438d-a198-0329b1060dbf",
   "metadata": {},
   "outputs": [],
   "source": [
    "NAME = \"\"\n",
    "COLLABORATORS = \"\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1801061f-7435-4e17-88bf-5ba35a2bda11",
   "metadata": {},
   "source": [
    "# 2025-11-14 Numerical Differentiation\n",
    "\n",
    "This activity will explore numerical approximation of derivatives using methods of the from\n",
    "\n",
    "$$ f' \\left( x \\right) \\approx \\frac{f \\left( x + h \\right) - f \\left( x \\right)}{h} $$\n",
    "\n",
    "for some small $h > 0$.\n",
    "This is know as the **first-order forward difference**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4d9f7188-fde5-40da-89bc-c435316333dd",
   "metadata": {},
   "outputs": [],
   "source": [
    "using LinearAlgebra\n",
    "using Plots\n",
    "using Test\n",
    "default(linewidth=3)\n",
    "\n",
    "# Approximate the derivative with a forward finite difference\n",
    "function diff_fwd(f, x, h)\n",
    "    (f.(x .+ h) - f.(x)) / h\n",
    "end\n",
    "\n",
    "# And plot\n",
    "plot([cos,\n",
    "      x -> diff_fwd(sin, x, .1),\n",
    "      x -> diff_fwd(sin, x, .5)],\n",
    "    label=[\"\\$cos ( x )\\$\" \"\\$h = 0.1\\$\" \"\\$h = 0.5\\$\"],\n",
    "    xlims=(-4, 4))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "44a4d776-795a-4e89-9803-62b607ec545f",
   "metadata": {},
   "source": [
    "In this figure, every point in the $h = 0.1$ and $h = 0.5$ lines is calculated by a finite difference estimate of the derivative.\n",
    "The analytic derivative is shown in blue.\n",
    "The $h = 0.1$ line appears to be more accurate than $h = 0.5$, which is a good sign.\n",
    "\n",
    "To be quantitative, let's examine how the error (taken to be the maximum error on a set of sample points) depends upon the parameter $h$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "28ff658c-169b-42d7-b168-75b3792f4e6e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute the error from the true derivative\n",
    "function compute_errors(diff, f, fprime; x=LinRange(-4, 4, 50))\n",
    "    hs = 10. .^ (-15:.5:0)\n",
    "    errors = [norm(diff(f, x, h) - fprime.(x)) for h in hs]\n",
    "    hs, errors\n",
    "end\n",
    "\n",
    "# And let's plot\n",
    "hs, errors = compute_errors(diff_fwd, sin, cos)\n",
    "scatter(hs, errors, label=\"forward diff\", xlabel=\"\\$h\\$\", ylabel=\"error\", xscale=:log10, yscale=:log10)\n",
    "plot!(hs, hs, label=\"\\$O(h)\\$\")\n",
    "plot!(hs, 1e-15 ./ hs, label=\"\\$O(1/h)\\$\", legend=:bottomleft)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "68fcf3f3-83bb-49ce-88cd-ea3c84b201de",
   "metadata": {},
   "source": [
    "Evidently the error behaves like $\\mathcal{O} \\left( h \\right)$ for large values of $h$ and like $\\mathcal{O} \\left( 1 / h \\right)$ for small values of $h$.\n",
    "To minimize error, we should choose a value of $h \\approx 10^{-8} \\approx \\sqrt{\\epsilon_\\text{machine}}$ in this case, but the specific choice depends on the function."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "785d20c8-d3a1-4676-bdb2-cad301515fc9",
   "metadata": {},
   "source": [
    "## Backward difference\n",
    "\n",
    "There is nothing special about adding $h$.\n",
    "We can also consider the **first-ordre backward difference**\n",
    "\n",
    "$$ f' \\left( x \\right) \\approx \\frac{f \\left( x \\right) - f \\left( x - h \\right)}{h} $$\n",
    "\n",
    "for some small $h > 0$.\n",
    "\n",
    "Complete the function below and compare with our forward difference."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2f1dfaf0-796e-456a-b386-9aac2759b734",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Approximate the derivative with a backward finite difference.\n",
    "function diff_back(f, x, h)\n",
    "    ### BEGIN SOLUTION\n",
    "\n",
    "    ### END SOLUTION\n",
    "end\n",
    "\n",
    "# And let's plot\n",
    "hs, errors = compute_errors(diff_back, sin, cos)\n",
    "scatter(hs, errors, label=\"backward diff\", xlabel=\"\\$h\\$\", ylabel=\"error\", xscale=:log10, yscale=:log10)\n",
    "plot!(hs, hs, label=\"\\$O(h)\\$\")\n",
    "plot!(hs, 1e-15 ./ hs, label=\"\\$O(1/h)\\$\", legend=:bottomleft)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "75aaf6a7-8bd4-4689-b978-c475ea4ef9c4",
   "metadata": {},
   "source": [
    "It appears that we don't see particularly different behavior between the forward and backward difference equations.\n",
    "Perhaps if we combined them?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c991e71c-8c28-44b3-96ac-80e7f849039f",
   "metadata": {},
   "source": [
    "## Centered difference formula\n",
    "\n",
    "We can do similar experiments for the \"centered\" difference formula\n",
    "\n",
    "$$ f' \\left( x \\right) \\approx \\frac{f \\left( x + h \\right) - f \\left( x - h \\right)}{2 h} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ec5605e9-017e-4d48-b4cd-3386e73ce3ed",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Finite difference with centered difference formula\n",
    "function diff_cent(f, x, h)\n",
    "    ### BEGIN SOLUTION\n",
    "\n",
    "    ### END SOLUTION\n",
    "end\n",
    "\n",
    "hs, errors = compute_errors(diff_cent, sin, cos)\n",
    "scatter(hs, errors, label=\"centered diff\",\n",
    "    xlabel=\"\\$h\\$\", ylabel=\"error\",\n",
    "    xscale=:log10, yscale=:log10, ylims=(1e-12, 10))\n",
    "plot!(hs, hs.^2, label=\"\\$O(h^2)\\$\")\n",
    "plot!(hs, 1e-15 ./ hs, label=\"\\$O(1/h)\\$\", legend=:bottomleft)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5fbb591c-8ac6-467e-aba6-e62c53cf8b6e",
   "metadata": {},
   "source": [
    "### Caption\n",
    "\n",
    "Write a caption for this figure, explaining the effect of **rounding error** and **discretization error**.\n",
    "Explain how to use the estimates $1/h$, $h$, and/or $h^2$ to produce an expression involving $\\epsilon_\\text{machine}$ to justify a good choice of $h$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "517c095a-97aa-493f-85df-0197886c6981",
   "metadata": {},
   "source": [
    "### BEGIN SOLUTION\n",
    "\n",
    "### END SOLUTION"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a7d70a37-8404-4c44-ad4f-d6cde5b88d4e",
   "metadata": {},
   "source": [
    "## Shifting the domain\n",
    "\n",
    "Now, we consider the effect of shifting from comparing accuracy on the interval $\\left( -4, 4 \\right)$ to the interval $\\left( 10^5 \\pi - 4, 10^5 \\pi + 4 \\right)$.\n",
    "Note that our test function $f \\left( x \\right) = \\sin \\left( x \\right)$ is $2 \\pi$ periodic."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ed4da580-910b-4520-99b4-4422f57d4005",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Gather our data\n",
    "hs, errors1 = compute_errors(diff_fwd,  sin, cos, x=1e5*π.+LinRange(-4, 4, 50))\n",
    "_,  errors2 = compute_errors(diff_cent, sin, cos, x=1e5*π.+LinRange(-4, 4, 50))\n",
    "\n",
    "# And plot\n",
    "scatter(hs, errors1, label=\"forward diff\",\n",
    "    xlabel=\"\\$h\\$\", ylabel=\"error\",\n",
    "    xscale=:log10, yscale=:log10, ylims=(1e-12, 10))\n",
    "scatter!(hs, errors2, label=\"centered diff\")\n",
    "plot!(hs, hs, label=\"\\$O(h)\\$\")\n",
    "plot!(hs, hs.^2, label=\"\\$O(h^2)\\$\")\n",
    "plot!(hs, 1e-15 ./ hs, label=\"\\$O(1/h)\\$\", legend=:bottomleft)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b3db2c25-a76a-49d2-9d77-b100352a58b8",
   "metadata": {},
   "source": [
    "This definitely is higher error than we saw before.\n",
    "Create a caption discussing the difference between the plots above and this most recent plot."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4686a28-3e67-447a-b998-68e6db9314b3",
   "metadata": {},
   "source": [
    "### BEGIN SOLUTION\n",
    "\n",
    "### END SOLUTION"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "75905fea-b87d-458a-bf9a-12ccd77cc3e6",
   "metadata": {},
   "source": [
    "## Estimate constants\n",
    "\n",
    "Suppose you are given two positive values $h_0$ and $h_1$ as well as two errors $e_0 = e \\left( h_0 \\right)$ and $e_1 = e \\left( h_1 \\right)$.\n",
    "Use the ratios $h_1 / h_0$ and $e_1 / e_0$ to estimate $p$ assuming our error is of the form $e \\left( h \\right) = c h^p$.\n",
    "\n",
    "Hint: take ratios of the equations and use a logarithm (`log`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a20ec8b6-5269-492c-800f-a7d0373c0960",
   "metadata": {},
   "outputs": [],
   "source": [
    "function convergence_order(h0, e0, h1, e1)\n",
    "    hratio = h1 / h0\n",
    "    eratio = e1 / e0\n",
    "    ### BEGIN SOLUTION\n",
    "\n",
    "    ### END SOLUTION\n",
    "end\n",
    "\n",
    "convergence_order(hs[end-6], errors1[end-6], hs[end-7], errors1[end-7])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b75761d2-4a14-4aab-a78e-4abac559705f",
   "metadata": {},
   "outputs": [],
   "source": [
    "@test abs(convergence_order(hs[end-5], errors1[end-5], hs[end-7], errors1[end-7]) - 1) < 1e-4\n",
    "@test abs(convergence_order(hs[end-4], errors2[end-4], hs[end-3], errors2[end-3]) - 2) < 1e-4"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "303aed43-444e-4104-9e5a-2d4417a8ec4d",
   "metadata": {},
   "source": [
    "Note: While the last assignment gave insight into some ideas behind the [**Finite Element Method**](https://en.wikipedia.org/wiki/Finite_element_method), this assignment gives insights to [**Finite Difference Methods**](https://en.wikipedia.org/wiki/Finite_difference_method)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e7a26ea-e1b6-45fb-90c0-d39e815f3365",
   "metadata": {},
   "source": [
    "## Extension\n",
    "\n",
    "We can iteratively apply these ideas to approximate higher order derivatives.\n",
    "Read the [Wikipedia article](https://en.wikipedia.org/wiki/Finite_difference#Higher-order_differences) about higher order differences and pick one of the formulations for estimating the second derivative.\n",
    "That is to say, pick one of the formulations\n",
    "\n",
    "$$ f'' \\left( x \\right) \\approx \\dots $$\n",
    "\n",
    "Implement the function below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d400acd0-5726-4be8-9d53-40cb39c023d1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Second derivative via finite differencing\n",
    "function diff2(f, x, h)\n",
    "    ### BEGIN SOLUTION\n",
    "\n",
    "    ### END SOLUTION\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7f44a153-52d0-47f3-9018-0023c787d29d",
   "metadata": {},
   "source": [
    "Let's be bold here.\n",
    "Here's a function for computing the error in a second order.\n",
    "Make a plot like the ones above showing the error in this second order finite difference.\n",
    "Add appropriate reference lines ($h$, $h^2$, $1/h$, etc) that match the error we see."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15f02d8f-eb37-4d69-80c9-a6f7d57e96d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute the error from the true second derivative\n",
    "function compute_errors_diff2(diff2, f, fprimeprime; x=LinRange(-4, 4, 50))\n",
    "    hs = 10. .^ (-15:.5:0)\n",
    "    errors = [norm(diff2(f, x, h) - fprimeprime.(x)) for h in hs]\n",
    "    hs, errors\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31cee6a9-b534-4b03-8078-bf919c579f70",
   "metadata": {},
   "outputs": [],
   "source": [
    "### BEGIN SOLUTION\n",
    "\n",
    "### END SOLUTION"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2d439be6-cca4-4f12-adc1-55bd7b56fb57",
   "metadata": {},
   "source": [
    "Caption this plot as you did above."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5468272d-7e3b-4d07-98a6-f32841470c82",
   "metadata": {},
   "source": [
    "### BEGIN SOLUTION\n",
    "\n",
    "### END SOLUTION"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.11.6",
   "language": "julia",
   "name": "julia-1.11"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.11.6"
  },
  "nbsphinx": {
   "execute": "never"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}