{ "cells": [ { "cell_type": "markdown", "id": "299ff832-51c2-4f70-be5c-00c68f1f3d70", "metadata": {}, "source": [ "## Floating point numbers and numerical calculations" ] }, { "cell_type": "markdown", "id": "7928c2ef-dfd9-4c1f-a7ac-34a2f6b35277", "metadata": {}, "source": [ "### Instability of algorithms" ] }, { "cell_type": "markdown", "id": "6a87f8b0-a5e0-4da1-b9a5-5142a2270170", "metadata": {}, "source": [ "Let's consider the integrals \n", "\n", "$$I(n) = \\int_0^1 x^n e^{-x} \\mathrm{d}x , \\quad n = 1, 2, \\ldots$$\n", "\n", "Since the integrand is positive for $0 \\le x \\le 1$,\n", "\n", "$$I(n) > 0 .$$\n", "\n", "Moreover, since $x^{n+1} < x^n$ for $0 < x < 1$,\n", "\n", "$$I(n) > I(n+1).$$\n", "\n", "We'll use those properties of the integral to check the validity of our numerical calculations." ] }, { "cell_type": "markdown", "id": "9e6c39a5-ece2-42c5-800b-54afad0af5a5", "metadata": {}, "source": [ "Integrating by parts, we can establish the following recurrence relation\n", "\n", "$$I(n+1) = n I(n) - \\frac{1}{e}$$\n", "\n", "Also,\n", "\n", "$$I(1) = 1 - \\frac{2}{e} .$$\n", "\n", "We can use those relations to calculate $I(n)$ for any $n$." ] }, { "cell_type": "code", "execution_count": null, "id": "5464495e-09c9-4e21-b519-2e239c14582e", "metadata": {}, "outputs": [], "source": [ "\n", "using QuadGK\n", "using PyPlot" ] }, { "cell_type": "code", "execution_count": null, "id": "7622ae80-57a8-437a-8e12-04d548272344", "metadata": {}, "outputs": [], "source": [ "\n", "integrand(x, n) = x^n * exp(-x)" ] }, { "cell_type": "code", "execution_count": null, "id": "91cd8667-9e27-4a45-adc3-06206fd34444", "metadata": {}, "outputs": [], "source": [ "\n", "quadgk_count(x -> integrand(x, 1), 0.0, 1.0)" ] }, { "cell_type": "code", "execution_count": null, "id": "2d63e719-13ef-4e80-8e1b-c91257d441da", "metadata": {}, "outputs": [], "source": [ "\n", "np = 22\n", "res1 = zeros(np)\n", "res2 = zeros(np);" ] }, { "cell_type": "code", "execution_count": null, "id": "8f824241-aa2c-4d79-9000-c9b9f453cfc0", "metadata": {}, "outputs": [], "source": [ "\n", "@time for n = 1:np\n", " res1[n] = quadgk(x -> integrand(x, n), 0.0, 1.0)[1]\n", "end" ] }, { "cell_type": "code", "execution_count": null, "id": "ffe9a359-c877-43db-a74f-56b6499fe3a7", "metadata": {}, "outputs": [], "source": [ "\n", "res2[1] = 1.0 - 2*exp(-1.0)\n", "@time for n = 2:np\n", " res2[n] = n * res2[n-1] - exp(-1.0)\n", "end" ] }, { "cell_type": "code", "execution_count": null, "id": "a41d1415-0228-4bab-876d-68e5da2d60e2", "metadata": {}, "outputs": [], "source": [ "\n", "plot(1:np, res1, marker=\"o\", linestyle=\"none\", label=\"QuadGK\")\n", "plot(1:np, res2, marker=\"x\", linestyle=\"none\", label=\"Recurrence\")\n", "ylim(-0.25, 0.3)\n", "legend()\n", "grid(true)" ] }, { "cell_type": "code", "execution_count": null, "id": "fa59d0ad-645b-42c3-9ef2-e75bb7dfc285", "metadata": {}, "outputs": [], "source": [ "\n", "semilogy(1:np, abs.(res1 .- res2), marker=\"o\")\n", "grid(true)\n", "xlabel(\"n\")\n", "ylabel(\"absolute error\");" ] }, { "cell_type": "markdown", "id": "89180f28-d991-4e87-bf60-d31a944757bf", "metadata": {}, "source": [ "### Dertivatives using finite differences" ] }, { "cell_type": "markdown", "id": "a7649936-de41-4937-8200-8b1c01d5b17a", "metadata": {}, "source": [ "The test function and its analytic derivative" ] }, { "cell_type": "code", "execution_count": null, "id": "36dff565-5be9-4b2e-8c89-017a2acf09fb", "metadata": {}, "outputs": [], "source": [ "\n", "f(x) = exp(x) + sin(x)\n", "fp(x) = exp(x) + cos(x);" ] }, { "cell_type": "markdown", "id": "0eadbac4-c86d-44b9-b771-41c1e8c9c9bb", "metadata": {}, "source": [ "Approximations to the first derivatives" ] }, { "cell_type": "code", "execution_count": null, "id": "378fbcb0-c464-4768-b2d5-68cd308d2679", "metadata": {}, "outputs": [], "source": [ "forward(f, x, h) = (f(x + h) - f(x)) / h\n", "central2(f, x, h) = (f(x + h) - f(x - h)) / (2 * h)\n", "central4(f, x, h) = (-f(x + 2 * h) + 8 * f(x + h) - 8 * f(x - h) + f(x - 2 * h)) / (12 * h);" ] }, { "cell_type": "code", "execution_count": null, "id": "222c4eb8-aa3e-4c2e-9784-9ce1fe52e39c", "metadata": {}, "outputs": [], "source": [ "\n", "np = 53\n", "dx = 2.0 .^ ((-np):0)\n", "err_f = similar(dx)\n", "err_c2 = similar(dx)\n", "err_c4 = similar(dx);" ] }, { "cell_type": "code", "execution_count": null, "id": "1972332d-5c6f-421e-bf11-a1b4aebfde59", "metadata": {}, "outputs": [], "source": [ "\n", "x0 = 1.0\n", "exact_fpx0 = fp(x0) # exact value of the derivative at x0\n", "\n", "for i in 1:(np + 1)\n", " err_f[i] = abs(forward(f, x0, dx[i]) - exact_fpx0)\n", " err_c2[i] = abs(central2(f, x0, dx[i]) - exact_fpx0)\n", " err_c4[i] = abs(central4(f, x0, dx[i]) - exact_fpx0)\n", "end" ] }, { "cell_type": "code", "execution_count": null, "id": "e3fb30cb-00b7-47ee-be5a-1a41c9153430", "metadata": {}, "outputs": [], "source": [ "\n", "loglog(dx, err_f; marker=\".\", label=\"Forward [x, x+h]\")\n", "loglog(dx, err_c2; marker=\".\", label=\"Central [x-h, x+h]\")\n", "loglog(dx, err_c4; marker=\".\", label=\"Central [x-2h, x-h, x+h, x+2h]\")\n", "legend()\n", "grid(true)\n", "xlabel(\"h\")\n", "ylabel(\"absolute error\")\n", "ylim(1e-14, 1.0);" ] }, { "cell_type": "markdown", "id": "7909a137-cef8-4079-ab05-cab144889723", "metadata": {}, "source": [ "### (More) catastrophic cancellations" ] }, { "cell_type": "markdown", "id": "f88bc0db-bf6e-49f4-b2d1-2f79b5d94057", "metadata": {}, "source": [ "Consider the function\n", "\n", "$$\\phi(x) = \\frac{a^{\\frac{3}{2}}}{ \\sin(x)} \\, \\left(\\frac{1}{\\sqrt{a - x}} - \\frac{1}{\\sqrt{a + x}} \\right) ,$$\n", "\n", "where $a$ is a large constant.\n", "\n", "$\\phi$ and its derivative are continuous functions\n", "\n", "$$\\lim_{x \\to 0} = 1 .$$" ] }, { "cell_type": "code", "execution_count": null, "id": "3d61953c-8215-4265-89b2-b3327137dee1", "metadata": {}, "outputs": [], "source": [ "\n", "const a = 1e7\n", "phi(x) = (1/sqrt(a - x) - 1/sqrt(a + x))*a^(3/2) / sin(x)" ] }, { "cell_type": "markdown", "id": "420fdbe1-e8f1-476c-939e-f9ffe440a29a", "metadata": {}, "source": [ "However, for small $x$ we get:" ] }, { "cell_type": "code", "execution_count": null, "id": "3787b5ac-a13c-44d2-be12-7632bfc26005", "metadata": {}, "outputs": [], "source": [ "\n", "\n", "xx = range(1e-10, 1e-8, 30)\n", "yy = phi.(xx);" ] }, { "cell_type": "code", "execution_count": null, "id": "f9882b26-77e5-489a-8938-3f1cd2895970", "metadata": {}, "outputs": [], "source": [ "\n", "plot(xx, yy, marker=\".\", label=\"original expression\")\n", "grid(true)\n", "xlabel(L\"x\")\n", "ylabel(L\"\\phi(x)\")\n", "title(\"An example of catastrophic cancellations\");" ] }, { "cell_type": "markdown", "id": "8150c9ea-1968-4b37-b837-f213f265d750", "metadata": {}, "source": [ "On the other hand, if we rewrite $\\phi(x)$ as follows:\n", "\n", "$$\\phi(x) = \\frac{a^{\\frac{3}{2}}}{ \\sin(x)} \\, \\left(\\frac{1}{\\sqrt{a - x}} - \\frac{1}{\\sqrt{a + x}} \\right) =$$\n", "$$\\frac{a^{\\frac{3}{2}}}{ \\sin(x)} \\frac{\\sqrt{a + x} - \\sqrt{a - x}}{\\sqrt{a^2 - x^2}} = $$\n", "$$\\frac{a^{\\frac{3}{2}}}{ \\sin(x)} \\frac{\\left(\\sqrt{a + x} - \\sqrt{a - x}\\right)\n", "\\left(\\sqrt{a + x} + \\sqrt{a - x}\\right)}{\\sqrt{a^2 - x^2} \\left(\\sqrt{a + x} + \\sqrt{a - x}\\right)} = $$\n", "$$\\frac{a^{\\frac{3}{2}}}{ \\sin(x)} \\frac{ 2 \\, x}{\\sqrt{a^2 - x^2} \\left(\\sqrt{a + x} + \\sqrt{a - x}\\right)} = \\psi(x) .$$\n", "\n", "$\\phi(x)$ and $\\psi(x)$ is the same function writtent in two different forms." ] }, { "cell_type": "code", "execution_count": null, "id": "3140ee3b-500b-4c51-962c-651435ec5e34", "metadata": {}, "outputs": [], "source": [ "\n", "psi(x) = 2 * a^(3/2) * x / (sqrt(a^2 - x^2) * (sqrt(a + x) + sqrt(a - x)) * sin(x))" ] }, { "cell_type": "code", "execution_count": null, "id": "ade76676-5fd6-40e5-80d5-7cf176cea00b", "metadata": {}, "outputs": [], "source": [ "\n", "zz = psi.(xx);" ] }, { "cell_type": "code", "execution_count": null, "id": "650d097d-8c87-4c4d-a67c-248f6b58909a", "metadata": {}, "outputs": [], "source": [ "\n", "plot(xx, yy, marker=\".\", label=L\"$\\phi(x)$\")\n", "plot(xx, zz, marker=\".\", label=L\"$\\psi(x)$\")\n", "grid(true)\n", "xlabel(\"x\")\n", "legend()\n", "title(\"Catastrophic cancellations is mitigated\");" ] }, { "cell_type": "markdown", "id": "c4a03d69-b310-449a-aade-7b6a58bc6ad9", "metadata": {}, "source": [ "### Integer overflow" ] }, { "cell_type": "code", "execution_count": null, "id": "f0bc76df-add8-4ef5-91cd-71aab425d1bf", "metadata": {}, "outputs": [], "source": [ "10^20" ] }, { "cell_type": "code", "execution_count": null, "id": "923e660e-ce7e-4dc9-b07b-54af541c7a87", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Julia 1.10.6", "language": "julia", "name": "julia-1.10" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.10.6" } }, "nbformat": 4, "nbformat_minor": 5 }