{ "cells": [ { "cell_type": "markdown", "id": "9086ab74-ff6f-426e-8930-7ae1b606c155", "metadata": {}, "source": [ "# Toolbox for Asymptotics with Explicit Error Bounds" ] }, { "cell_type": "code", "execution_count": 1, "id": "14372689-d5ad-426f-b959-d051ca3b2e53", "metadata": {}, "outputs": [], "source": [ "import dependent_bterms as dbt" ] }, { "cell_type": "code", "execution_count": 2, "id": "1e3e651e-9760-482c-bfc5-a3154acd88e5", "metadata": {}, "outputs": [], "source": [ "%display typeset" ] }, { "cell_type": "markdown", "id": "7b3ad651-ed83-491e-bf12-d4277fbcd990", "metadata": {}, "source": [ "## Creating a special `AsymptoticRing`\n", "\n", "See [the documentation of SageMath's module on asymptotic expansions](https://doc.sagemath.org/html/en/reference/asymptotic/sage/rings/asymptotic/asymptotic_ring.html) for a primer on the `AsymptoticRing`.\n", "\n", "Assume that we want to carry out computations in a setting where we have a variable $n\\to\\infty$, as well as a dependent variable $k = k(n)$ for which we know $n^{\\alpha} \\leq k \\leq n^{\\beta}$ for some $0\\leq\\alpha\\leq\\beta$. Concretely, let us consider $1 = n^0 \\leq k\\leq n^{4/7}$.\n", "\n", "**Note:** This implementation currently only supports monomial asymptotic growth of the independent variable. Trying to use it with more intricate growth groups other than $n^{\\mathbb{Q}}$ will likely result in unwanted behavior." ] }, { "cell_type": "code", "execution_count": 3, "id": "ad853bea-c7ad-43fd-9ed8-0c6a62782894", "metadata": {}, "outputs": [], "source": [ "AR, n, k = dbt.AsymptoticRingWithDependentVariable(\n", " \"n^QQ\", \"k\", 0, 4 / 7, bterm_round_to=3, default_prec=5\n", ")" ] }, { "cell_type": "markdown", "id": "8427d2ac-fa8d-4482-baaf-a9d3cfc08dea", "metadata": {}, "source": [ "Technically, `k` is a plain symbolic variable -- the instantiated Asymptotic Ring uses the Symbolic Ring as its coefficient ring:" ] }, { "cell_type": "code", "execution_count": 4, "id": "e0bf9bf7-bf96-44ef-beda-27e33533bcde", "metadata": {}, "outputs": [ { "data": { "text/html": [ "\$\\displaystyle \\verb|Asymptotic|\\verb| |\\verb|Ring|\\verb| |\\verb|<n^QQ>|\\verb| |\\verb|over|\\verb| |\\verb|Symbolic|\\verb| |\\verb|Ring|\$" ], "text/latex": [ "$\\displaystyle \\verb|Asymptotic|\\verb| |\\verb|Ring|\\verb| |\\verb||\\verb| |\\verb|over|\\verb| |\\verb|Symbolic|\\verb| |\\verb|Ring|$" ], "text/plain": [ "Asymptotic Ring over Symbolic Ring" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "AR" ] }, { "cell_type": "code", "execution_count": 5, "id": "ead465bc-41f8-4188-9b97-5707a51525cc", "metadata": {}, "outputs": [ { "data": { "text/html": [ "\$\\displaystyle \\mathrm{True}\$" ], "text/latex": [ "$\\displaystyle \\mathrm{True}$" ], "text/plain": [ "True" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "k in SR" ] }, { "cell_type": "markdown", "id": "bde8805f-79c6-47ce-8106-fc458ae486bd", "metadata": {}, "source": [ "## Arithmetic with the dependent variable\n", "\n", "Arithmetic with the expansions from our special Asymptotic Rign works *as usual*:" ] }, { "cell_type": "code", "execution_count": 6, "id": "82c7154d-c435-42b2-a6ec-2502c210e4cd", "metadata": {}, "outputs": [ { "data": { "text/html": [ "\$\\displaystyle 3 n + 127 + 42 n^{-1} + 12 n^{-\\frac{4}{3}} + 4 n^{-\\frac{7}{3}}\$" ], "text/latex": [ "$\\displaystyle 3 n + 127 + 42 n^{-1} + 12 n^{-\\frac{4}{3}} + 4 n^{-\\frac{7}{3}}$" ], "text/plain": [ "3*n + 127 + 42*n^(-1) + 12*n^(-4/3) + 4*n^(-7/3)" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "(1 + 3 * n) * (4 * n ^ (-7 / 3) + 42 / n + 1)" ] }, { "cell_type": "code", "execution_count": 7, "id": "2cddbe0d-954d-440e-88ae-a3106f7ec6e9", "metadata": {}, "outputs": [ { "data": { "text/html": [ "\$\\displaystyle 1 + n^{-1} + n^{-2} + 2 n^{-3} + 2 n^{-4} + 3 n^{-5} + 4 n^{-6} + 5 n^{-7} + 6 n^{-8} + 8 n^{-9} + O\\!\\left(n^{-10}\\right)\$" ], "text/latex": [ "$\\displaystyle 1 + n^{-1} + n^{-2} + 2 n^{-3} + 2 n^{-4} + 3 n^{-5} + 4 n^{-6} + 5 n^{-7} + 6 n^{-8} + 8 n^{-9} + O\\!\\left(n^{-10}\\right)$" ], "text/plain": [ "1 + n^(-1) + n^(-2) + 2*n^(-3) + 2*n^(-4) + 3*n^(-5) + 4*n^(-6) + 5*n^(-7) + 6*n^(-8) + 8*n^(-9) + O(n^(-10))" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "prod((1 + n ^ (-j)) for j in srange(1, 10)) * (1 + O(n ^ (-10)))" ] }, { "cell_type": "code", "execution_count": 8, "id": "67c2dfce-34cc-44c8-aa05-f5913f478bbd", "metadata": {}, "outputs": [ { "data": { "text/html": [ "\$\\displaystyle 1 + n^{-1} + n^{-2} + n^{-3} + n^{-4} + O\\!\\left(n^{-5}\\right)\$" ], "text/latex": [ "$\\displaystyle 1 + n^{-1} + n^{-2} + n^{-3} + n^{-4} + O\\!\\left(n^{-5}\\right)$" ], "text/plain": [ "1 + n^(-1) + n^(-2) + n^(-3) + n^(-4) + O(n^(-5))" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "n / (n - 1)" ] }, { "cell_type": "markdown", "id": "a9b6e0ea-6520-4a50-87ca-ae21d64c9fe4", "metadata": {}, "source": [ "We can also simply use the dependent variable." ] }, { "cell_type": "code", "execution_count": 9, "id": "cc95d7a3-6bbd-4332-872b-ea1f51f8e933", "metadata": {}, "outputs": [ { "data": { "text/html": [ "\$\\displaystyle k^{3} n + k n^{2} + O\\!\\left(n^{\\frac{3}{2}}\\right)\$" ], "text/latex": [ "$\\displaystyle k^{3} n + k n^{2} + O\\!\\left(n^{\\frac{3}{2}}\\right)$" ], "text/plain": [ "k^3*n + k*n^2 + O(n^(3/2))" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "expr = k * n ^ 2 + O(n ^ (3 / 2)) + k ^ 3 * n\n", "expr" ] }, { "cell_type": "markdown", "id": "c114cac4-1e0e-4e32-985a-e10a2494361c", "metadata": {}, "source": [ "Notice that the term order is now constructed with respect to the highest potential growth of a summand \n", "(i.e., if $k = k(n)$ were to assume its upper bound). We can ask the summands to reveal their respective growth bounds; the following tuples describe the lower and the upper bound of a summand, respectively:" ] }, { "cell_type": "code", "execution_count": 10, "id": "f6b90d67-2b91-4d02-954f-028ccbfeffd3", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "O(n^(3/2)) -> (n^(3/2), n^(3/2))\n", "k*n^2 -> (n^2, n^(18/7))\n", "k^3*n -> (n, n^(19/7))\n" ] } ], "source": [ "for summand in expr.summands.elements_topological():\n", " print(f\"{summand} -> {summand.dependent_growth_range()}\")" ] }, { "cell_type": "markdown", "id": "882fef0b-0959-49d0-8df3-a6837ab23afa", "metadata": {}, "source": [ "Automatic expansions work too:" ] }, { "cell_type": "code", "execution_count": 11, "id": "0b5491c5-5580-462d-b694-ec025abb7077", "metadata": {}, "outputs": [ { "data": { "text/html": [ "\$\\displaystyle 1 + \\left(k + 1\\right) n^{-1} + \\left(\\frac{1}{2} \\, {\\left(k + 1\\right)}^{2}\\right) n^{-2} + \\left(\\frac{1}{6} \\, {\\left(k + 1\\right)}^{3}\\right) n^{-3} + \\left(\\frac{1}{24} \\, {\\left(k + 1\\right)}^{4}\\right) n^{-4} + O\\!\\left(n^{-\\frac{15}{7}}\\right)\$" ], "text/latex": [ "$\\displaystyle 1 + \\left(k + 1\\right) n^{-1} + \\left(\\frac{1}{2} \\, {\\left(k + 1\\right)}^{2}\\right) n^{-2} + \\left(\\frac{1}{6} \\, {\\left(k + 1\\right)}^{3}\\right) n^{-3} + \\left(\\frac{1}{24} \\, {\\left(k + 1\\right)}^{4}\\right) n^{-4} + O\\!\\left(n^{-\\frac{15}{7}}\\right)$" ], "text/plain": [ "1 + (k + 1)*n^(-1) + (1/2*(k + 1)^2)*n^(-2) + (1/6*(k + 1)^3)*n^(-3) + (1/24*(k + 1)^4)*n^(-4) + O(n^(-15/7))" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "auto_expansion = exp((1 + k) / n)\n", "auto_expansion" ] }, { "cell_type": "markdown", "id": "63068238-4b68-4833-b1f9-79eca6665a55", "metadata": {}, "source": [ "Note that the error term, $O(n^{-15/7})$, would be able to absorb some parts of the exact terms, after expanding the powers of $(k+1)$. This is not done automatically, but can easily be triggered via a utility function from our toolbox:" ] }, { "cell_type": "code", "execution_count": 12, "id": "2f39f0e4-67e8-44b7-8e97-1335023e9330", "metadata": {}, "outputs": [ { "data": { "text/html": [ "\$\\displaystyle 1 + \\left(k + 1\\right) n^{-1} + \\left(\\frac{1}{2} \\, k^{2} + k + \\frac{1}{2}\\right) n^{-2} + \\left(\\frac{1}{6} \\, k^{3} + \\frac{1}{2} \\, k^{2}\\right) n^{-3} + \\frac{1}{24} \\, k^{4} n^{-4} + O\\!\\left(n^{-\\frac{15}{7}}\\right)\$" ], "text/latex": [ "$\\displaystyle 1 + \\left(k + 1\\right) n^{-1} + \\left(\\frac{1}{2} \\, k^{2} + k + \\frac{1}{2}\\right) n^{-2} + \\left(\\frac{1}{6} \\, k^{3} + \\frac{1}{2} \\, k^{2}\\right) n^{-3} + \\frac{1}{24} \\, k^{4} n^{-4} + O\\!\\left(n^{-\\frac{15}{7}}\\right)$" ], "text/plain": [ "1 + (k + 1)*n^(-1) + (1/2*k^2 + k + 1/2)*n^(-2) + (1/6*k^3 + 1/2*k^2)*n^(-3) + 1/24*k^4*n^(-4) + O(n^(-15/7))" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dbt.simplify_expansion(auto_expansion)" ] }, { "cell_type": "markdown", "id": "c75519af-baab-4af7-814d-13f229bb1d55", "metadata": {}, "source": [ "Now only the certified lower-order terms remain." ] }, { "cell_type": "markdown", "id": "91c36883-7b4f-41df-9597-d705455cbd27", "metadata": {}, "source": [ "## Expansions with explicit error bounds: B-terms\n", "\n", "B-terms are, in a nutshell, O-terms with explicitly specified bound constant and validity point. The term\n", "$$ B_{n\\geq 10}(42 n^{3}) $$\n", "represents an error term, valid from $n \\geq 10$, that is bounded by $42 n^3$. Basic support for arithmetic (over univariate, monomial growth groups) is shipped with a standard installation of SageMath. Our specialized asymptotic ring extends this functionality to expansions involving the dependent variable $k = k(n)$." ] }, { "cell_type": "code", "execution_count": 13, "id": "dc0a490b-404f-4d5a-943b-b5b46f88cd08", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "dependent_bterms/structures.py:362: FutureWarning: This class/method/function is marked as experimental. It, its functionality or its interface might change without a formal deprecation.\n", "See https://github.com/sagemath/sage/issues/31922 for details.\n", " super().__init__(\n" ] }, { "data": { "text/html": [ "\$\\displaystyle 7 n + B_{n \\ge 10}\\left(\\frac{53}{10} n^{-1}\\right)\$" ], "text/latex": [ "$\\displaystyle 7 n + B_{n \\ge 10}\\left(\\frac{53}{10} n^{-1}\\right)$" ], "text/plain": [ "7*n + B(53/10*n^(-1), n >= 10)" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "7 * n + AR.B(5 / n, valid_from=10) + 3 / n ^ 2" ] }, { "cell_type": "markdown", "id": "5a1d8fb8-c383-4046-8c29-c87b5a78eab6", "metadata": {}, "source": [ "The constant $53/10$ is obtained from automatically estimating\n", "$$ \\frac{3}{n^2} = \\frac{3}{n\\cdot n} \\leq \\frac{3}{10 n}, $$\n", "and then letting the existing B-term *absorb* this bound." ] }, { "cell_type": "markdown", "id": "f8a7d10f-ac6a-4022-a89b-e86eca9e1a26", "metadata": {}, "source": [ "To avoid the accumulation of complicated symbolic expressions appearing when carrying out this automatic estimate, we have specified (via `bterm_round_to`) above that B-term constants should be rounded:" ] }, { "cell_type": "code", "execution_count": 14, "id": "d1dab575-f523-46cc-be5f-3ac20cf50030", "metadata": {}, "outputs": [ { "data": { "text/html": [ "\$\\displaystyle B_{n \\ge 10}\\left(\\frac{293}{200} n^{-1}\\right)\$" ], "text/latex": [ "$\\displaystyle B_{n \\ge 10}\\left(\\frac{293}{200} n^{-1}\\right)$" ], "text/plain": [ "B(293/200*n^(-1), n >= 10)" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "AR.B(1 / n, valid_from=10) + n ^ (-4 / 3)" ] }, { "cell_type": "markdown", "id": "a1a0a937-dce0-4f6d-ac72-a450d05235b3", "metadata": {}, "source": [ "where the constant is $\\big\\lceil (1 + \\frac{1}{10^{1/3}}) \\cdot 10^3\\big\\rceil \\cdot 10^{-3} = \\frac{293}{200}$. This mechanism significantly improves performance for large expansions." ] }, { "cell_type": "markdown", "id": "db3c5b53-aaaa-4267-9888-2646fff63e4d", "metadata": {}, "source": [ "For B-terms to work with the dependent variable, some special absorption rules are in place:" ] }, { "cell_type": "code", "execution_count": 15, "id": "b59d9006-48fd-4696-8e1a-46dec9e027a7", "metadata": {}, "outputs": [ { "data": { "text/html": [ "\$\\displaystyle B_{n \\ge 10}\\left(\\frac{3373}{1000} \\, {\\left| k^{2} \\right|} n^{-3}\\right)\$" ], "text/latex": [ "$\\displaystyle B_{n \\ge 10}\\left(\\frac{3373}{1000} \\, {\\left| k^{2} \\right|} n^{-3}\\right)$" ], "text/plain": [ "B(3373/1000*abs(k^2)*n^(-3), n >= 10)" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "AR.B(3 * k ^ 2 / n ^ 3, valid_from=10) + (1 - 2 * k + 3 * k ^ 2 - 4 * k ^ 3) / n ^ 5" ] }, { "cell_type": "markdown", "id": "473d61dc-c768-4201-8563-ef957a1cd555", "metadata": {}, "source": [ "We first estimate\n", "$$ \n", " \\bigg\\lvert\\frac{1 - 2k + 3k - 4k^3}{n^5}\\bigg\\rvert\n", " \\leq \\frac{(1 + 2 + 3 + 4) k ^3}{n^5} \n", " = \\frac{10 k^3}{n^5}.\n", "$$\n", "As the power of $k$ in this bound is larger than the maximal power of $k$ in the B-term,\n", "we may not yet proceed as above (otherwise we would increase the upper bound of the B-term,\n", "which we must avoid). Instead, we use first use $k\\leq n^{4/7}$, followed by\n", "$n\\geq 10$ to obtain\n", "$$ \n", " \\frac{10 k^3}{n^5} \\leq \\frac{10 k^2 n^{4/7}}{n^5} = \\frac{10 k^2}{n^{31/7}}\n", " \\leq \\frac{10 k^2}{10^{10/7} \\cdot n^3} = 10^{-3/7} \\cdot \\frac{k^2}{n^3},\n", "$$\n", "which the B-term now can absorb directly. Hence the value of the B-term constant\n", "is determined via $\\lceil (3 + 10^{-3/7}) \\cdot 10^3 \\rceil \\cdot 10^{-3} = \\frac{3373}{1000}$." ] }, { "cell_type": "markdown", "id": "9525c2b2-ce86-4bd8-8d47-ede884acd457", "metadata": {}, "source": [ "This module also ships with support for B-term bounded Taylor expansions. From a technical point of view, the error term (in its Lagrange form) is determined rigorously using interval arithmetic." ] }, { "cell_type": "code", "execution_count": 16, "id": "8aac3254-cf3d-42b7-adf1-36434d59e76f", "metadata": {}, "outputs": [ { "data": { "text/html": [ "\$\\displaystyle 1 + \\left({\\left(k + 1\\right)}^{2}\\right) n^{-2} + B_{n \\ge 10}\\left(\\left({\\left| \\frac{7351}{250} \\, k^{3} + 30 \\, k^{2} + 30 \\, k + 10 \\right|}\\right) n^{-3}\\right)\$" ], "text/latex": [ "$\\displaystyle 1 + \\left({\\left(k + 1\\right)}^{2}\\right) n^{-2} + B_{n \\ge 10}\\left(\\left({\\left| \\frac{7351}{250} \\, k^{3} + 30 \\, k^{2} + 30 \\, k + 10 \\right|}\\right) n^{-3}\\right)$" ], "text/plain": [ "1 + ((k + 1)^2)*n^(-2) + B((abs(7351/250*k^3 + 30*k^2 + 30*k + 10))*n^(-3), n >= 10)" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "arg = (1 + k) / n + AR.B(k ^ 3 / n ^ 3, valid_from=10)\n", "ex = dbt.taylor_with_explicit_error(\n", " lambda t: 1 / (1 - t ^ 2), arg, order=3, valid_from=10\n", ")\n", "ex" ] }, { "cell_type": "markdown", "id": "a39882bd-6241-4d73-aae3-05c538ab5ed9", "metadata": {}, "source": [ "Simplification (involving partial absorption of coefficients) may lead to situations where terms of technically weaker growth cannot (as a limitation of the current implementation) be absorbed by the constructed error term. In this case the returned expansion is still correct; just not as compact as it could be." ] }, { "cell_type": "code", "execution_count": 17, "id": "71d3d8a1-beef-4368-891e-e37ea15b70fc", "metadata": {}, "outputs": [ { "data": { "text/html": [ "\$\\displaystyle 1 + k^{2} n^{-2} + B_{n \\ge 10}\\left(\\left({\\left| \\frac{7351}{250} \\, k^{3} + 30 \\, k^{2} + 30 \\, k + 10 \\right|}\\right) n^{-3}\\right) + \\left(2 \\, k + 1\\right) n^{-2}\$" ], "text/latex": [ "$\\displaystyle 1 + k^{2} n^{-2} + B_{n \\ge 10}\\left(\\left({\\left| \\frac{7351}{250} \\, k^{3} + 30 \\, k^{2} + 30 \\, k + 10 \\right|}\\right) n^{-3}\\right) + \\left(2 \\, k + 1\\right) n^{-2}$" ], "text/plain": [ "1 + k^2*n^(-2) + B((abs(7351/250*k^3 + 30*k^2 + 30*k + 10))*n^(-3), n >= 10) + (2*k + 1)*n^(-2)" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dbt.simplify_expansion(ex)" ] }, { "cell_type": "markdown", "id": "296cbbc2-dabf-4858-a726-71fa64cbb1fb", "metadata": {}, "source": [ "By specifying `simplify_bterm_growth=True`, simplifying an expansion also collapses the dependent variables in all B-terms, resulting in an expansion involving a single, \"absolute\" B-term." ] }, { "cell_type": "code", "execution_count": 18, "id": "3bb1e9d1-9536-48f5-aec4-771467b9827a", "metadata": {}, "outputs": [ { "data": { "text/html": [ "\$\\displaystyle 1 + k^{2} n^{-2} + B_{n \\ge 10}\\left(\\frac{41441}{1000} n^{-\\frac{9}{7}}\\right)\$" ], "text/latex": [ "$\\displaystyle 1 + k^{2} n^{-2} + B_{n \\ge 10}\\left(\\frac{41441}{1000} n^{-\\frac{9}{7}}\\right)$" ], "text/plain": [ "1 + k^2*n^(-2) + B(41441/1000*n^(-9/7), n >= 10)" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dbt.simplify_expansion(ex, simplify_bterm_growth=True)" ] }, { "cell_type": "code", "execution_count": null, "id": "ee2a9d58-e8d9-4fe5-85fa-ef7567f02472", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "SageMath 10.3.beta8", "language": "sage", "name": "sagemath" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.12" } }, "nbformat": 4, "nbformat_minor": 5 }