{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "1426c8d8-5fab-47cb-b9f2-ab1157dc22a9",
   "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": "d64f8ca4-60de-4612-b615-74400a353916",
   "metadata": {},
   "outputs": [],
   "source": [
    "NAME = \"\"\n",
    "COLLABORATORS = \"\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "478481c8-260f-4edd-9a4d-eee78fb14b4e",
   "metadata": {},
   "source": [
    "# 2025-09-05 Taylor Series\n",
    "\n",
    "Suppose we wish to sum the series\n",
    "\n",
    "$$ \\sum_{n = 1}^{\\infty} \\frac{1}{n^2} = 1 + \\frac{1}{2^2} + \\frac{1}{3^2} + \\frac{1}{4^2} + \\cdots $$\n",
    "\n",
    "We cannot sum an infinite number of terms, but we can sum the first $k$ terms and see what happens as $k$ increases.\n",
    "In other words, we will compute\n",
    "\n",
    "$$ y_k = \\sum_{n = 1}^{k} \\frac{1}{n^2} $$\n",
    "\n",
    "and estimate what happens in the limit as $k \\rightarrow \\infty$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "90446b75-ce84-499f-860f-51658aff9683",
   "metadata": {},
   "outputs": [],
   "source": [
    "using Plots\n",
    "\n",
    "function series_n2(k)\n",
    "    sum = 0\n",
    "    for n in 1:k\n",
    "        sum += 1 / n^2\n",
    "    end\n",
    "    sum\n",
    "end\n",
    "\n",
    "series_n2(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d2e8137-8697-4d19-bbac-07e8720d7c6e",
   "metadata": {},
   "source": [
    "This [series converges](https://en.wikipedia.org/wiki/Pi#Irrationality_and_transcendence) to $\\pi^2 / 6$.\n",
    "(Proof is beyond the scope of this course.)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f35e6f75-a218-428d-9359-0edf355f498b",
   "metadata": {},
   "source": [
    "## How fast does the series converge?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "65d99360-af0a-4aab-a07c-e1dc53f2449a",
   "metadata": {},
   "outputs": [],
   "source": [
    "known_limit = π^2 / 6\n",
    "ks = 2 .^ (0:10) # [1, 2, 4, ..., 1024]\n",
    "plot(ks, [series_n2, k -> known_limit], label=[\"partial sum\" \"limit\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "04b96182-fbf4-44b1-8206-d7c7f5dcdd9e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Lets swap to a log scale\n",
    "plot(ks, [series_n2, k -> known_limit], xscale=:log10, label=[\"partial sum\" \"limit\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4bcb478c-f53c-4301-b54a-bb4ff96e4a9d",
   "metadata": {},
   "source": [
    "## Error\n",
    "\n",
    "If we have an approximation $y_k$ of an exact value $y_*$, then we can compute the **absolute error** as\n",
    "\n",
    "$$ y_k - y_* $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7e35cb9e-266a-47ab-ae05-7c7a31f2037a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot error\n",
    "plot(ks,  k -> series_n2(k) - known_limit, xscale=:log10, label=\"absolute error\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13daaab9-baaf-4278-8a6d-93e0a964af23",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Hmm, let's compare to a reference line\n",
    "plot(ks,\n",
    "    [k -> series_n2(k) - known_limit, k -> 1/k],\n",
    "    xscale=:log10,\n",
    "    label=[\"absolute error\" \"1/k\"]\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ec9cc4ba-190e-4e1a-b69c-6c819d9c1578",
   "metadata": {},
   "source": [
    "## Conclusion?\n",
    "\n",
    "This plot suggests that the error is proportional to $1 / k$, or\n",
    "\n",
    "$$ \\lvert \\sum_{n = 1}^{k} \\frac{1}{n^2} - \\pi^2 / 6 \\rvert \\sim \\frac{1}{k} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "678e37f4-8557-473b-928b-743d7a2ba4f1",
   "metadata": {},
   "source": [
    "### Notation\n",
    "\n",
    "The notation $f \\sim g$ means that\n",
    "\n",
    "$$ a g \\left( x \\right) \\leq f \\left( x \\right) \\leq b g \\left( x \\right) $$\n",
    "\n",
    "for two positive numbers $a$ and $b$ in some asymptotic regime ($x$ is sufficiently large or small).\n",
    "\n",
    "In this case, the 'asymptotic regime' is large $k$ (the variable '$x$' for this series), but sometimes we will also investigate small limits."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d05769d1-3925-4b19-8e7b-b4cd4c3948d4",
   "metadata": {},
   "source": [
    "## Case Study - exponential function\n",
    "\n",
    "Recall that [one definition](https://en.wikipedia.org/wiki/Exponential_function#Formal_definition) of the exponential function is\n",
    "\n",
    "$$ e^x = \\sum_{n = 0}^\\infty \\frac{x^n}{n!} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4acbe5b9-5fd9-4b93-9295-f6b43adf3d5b",
   "metadata": {},
   "source": [
    "### Modify this function\n",
    "\n",
    "Modify this function so at most `k` terms are summed.\n",
    "\n",
    "(Hint: `if` and `break` may help, or a modification to the `while` condition)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2ae67352-ab25-46a3-afe5-95e3b2e8f016",
   "metadata": {},
   "outputs": [],
   "source": [
    "function myexp(x, k)\n",
    "    sum = 0\n",
    "    term = 1\n",
    "    n = 1\n",
    "    # modify so at most k terms are summed\n",
    "    while sum + term != sum\n",
    "        sum += term\n",
    "        term *= x / n\n",
    "        ### BEGIN SOLUTION\n",
    "\n",
    "        ### END SOLUTION\n",
    "        n += 1\n",
    "    end\n",
    "    sum\n",
    "end\n",
    "\n",
    "myexp(1, 10)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2e22f88f-7437-47b0-929b-b2fb79a40625",
   "metadata": {},
   "source": [
    "### Investigation\n",
    "\n",
    "Why did the original loop terminate (before any changes you made)?\n",
    "Can you write an equivalent expression in terms of `sum`, `term`, and $\\epsilon_\\text{machine}$?\n",
    "(Discuss on Zulip)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9250162a-7ec7-43d8-9870-73cb1ca818b7",
   "metadata": {},
   "outputs": [],
   "source": [
    "@assert myexp(1, 100) ≈ exp(1)\n",
    "@assert myexp(1, 3) == 2.5\n",
    "@assert abs(myexp(1, 4)  - 2.6666666666666665) < eps()\n",
    "@assert abs(myexp(1, 10) - 2.7182815255731922) < eps()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5bf86df6-c480-493e-a33a-f326df223de6",
   "metadata": {},
   "source": [
    "## Relative Error\n",
    "\n",
    "If we have a function $\\tilde{f} \\left( x \\right)$ that approximates $f \\left( x \\right)$, then we can compute the **relative error** as\n",
    "\n",
    "$$ \\frac{\\tilde{f} \\left( x \\right) - f \\left( x \\right)}{\\lvert f \\left( x \\right) \\rvert} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0c5c21df-618f-4306-897d-67ef429cdbb3",
   "metadata": {},
   "outputs": [],
   "source": [
    "relative_error(x, k) = abs(myexp(x, k) - exp(x)) / exp(x)\n",
    "relative_error(1, 20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "39b574d1-69d2-4e39-af32-9bd8088d54f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot relative error as k increases\n",
    "plot(ks,\n",
    "    k -> relative_error(1, k),\n",
    "    xscale=:log10,\n",
    "    yscale=:log10,\n",
    "    xlabel=\"k\",\n",
    "    ylabel=\"relative error\",\n",
    "    label=:none\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bceb25a2-bce7-4274-bcc9-d3c0ddd537c4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot relative error as k increases\n",
    "plot(ks,\n",
    "    k -> relative_error(10, k),\n",
    "    xscale=:log10,\n",
    "    yscale=:log10,\n",
    "    xlabel=\"k\",\n",
    "    ylabel=\"relative error\",\n",
    "    label=:none\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2cfdd76b-1ffa-4e4c-9697-c0d5f344877a",
   "metadata": {},
   "source": [
    "## More Plotting\n",
    "\n",
    "Plot the relative error in `myexp(-20, x)` as `k` increases."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cb1d810e-d695-4354-b4f0-9c199db64f1b",
   "metadata": {},
   "outputs": [],
   "source": [
    "### BEGIN SOLUTION\n",
    "\n",
    "### END SOLUTION"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5b0dd39b-2096-4dc5-a9d0-630e8c168b7d",
   "metadata": {},
   "source": [
    "## Final Accuracy\n",
    "\n",
    "The intermediate accuracy can be quite bad.\n",
    "How good is the final result?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "92e0ea9a-6797-4b37-8540-9f849302a64d",
   "metadata": {},
   "outputs": [],
   "source": [
    "bottom = 1e-20 # To avoid having the plotter attempt log(0)\n",
    "plot(x -> relative_error(x, 1000) + bottom,\n",
    "    xlims=(-20, 20),\n",
    "    yscale=:log10,\n",
    "    xlabel=\"x\",\n",
    "    ylabel=\"relative error\",\n",
    "    label=:none\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a70a87c3-ab77-4550-b510-7bbf5f539359",
   "metadata": {},
   "source": [
    "## Exploration\n",
    "\n",
    "Select 1 or more of these questions to explore.\n",
    "\n",
    "1) Plot each of the terms in the  series expansion for computing $e^{-20}$.\n",
    "   * Where are rounding errors committed? Think about the size of those terms and the size of the final result ($e^{-20}$ is a small number).\n",
    "   * How would you explain the large relative errors using the concept of conditioning?\n",
    "   * How is this different different from computing $e^x$ with a positive argument?\n",
    "\n",
    "2) Use the identity $e^{-x} = 1 / e^x$ to write a new function that calls `myexp()` and is accurate for both positive and negative arguments `x`. Demonstrate this improvement with a plot and explain.\n",
    "\n",
    "3) Use the identity $e^x = \\left( e^{x / 2} \\right) ^2$ to reduce the number of terms that need to be summed. Recall that squaring is easy -- it is just multiplication. Demonstrate this improvement with a plot and explain."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2b288eb3-c69b-4401-9dc1-a6d02f67e0e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Exploration 1\n",
    "\n",
    "### BEGIN SOLUTION\n",
    "\n",
    "### END SOLUTION"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "da41319c-c6f7-4c54-87c7-bf098a7b6235",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Exploration 2\n",
    "\n",
    "### BEGIN SOLUTION\n",
    "\n",
    "### END SOLUTION"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4e29f94b-c06f-4608-8707-da6b2ba2c5af",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Exploration 3\n",
    "\n",
    "### 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
}