{ "cells": [ { "cell_type": "markdown", "id": "5a09a7f6", "metadata": {}, "source": [ "# Baker's algorithm" ] }, { "cell_type": "markdown", "id": "adef5768", "metadata": {}, "source": [ "Henri Padé, in his Ph.D. thesis [1], organized the Padé approximants in a $2\\times2$ infinite table similar to next one.\n", "\n", "$$\\begin{array}{c| c c c c c c c c c c c c }\n", "\t\t\t[p / q] & 0 & 1 & 2 & 3 &\\dots \\\\\n", "\t\t\t\\hline \\\\\n", "\t\t\t 0 & [0/0] & [0/1] & [0/2] & [0/3] &\\dots \\\\ \\\\\t\n", "\t\t\t 1 & [1/0] & [1/1] & [1/2] & [1/3] &\\dots \\\\ \t\\\\\n", "\t\t\t 2 & [2/0] & [2/1] & [2/2] & [2/3] &\\dots \\\\ \t\\\\\n", "\t\t\t 3 & [3/0] & [3/1] & [3/2] & [3/3] &\\dots \\\\ \t\t\n", "\t\t\t\\vdots& \\vdots &\\vdots &\\vdots&\\vdots & \\ddots \n", "\\end{array}$$ \n", "\n", "For different reasons we might be interested on particular blocks of the Padé table. On way of doing it is to use a recursive method. We now present Baker's recursive algorithm [2].\n", "\n", "Baker defined the sequences\n", "\\begin{equation*}\n", "\\begin{aligned}\n", "&\\frac{\\eta_{2j}(x)}{\\theta_{2j}(x)}=[p-j/j\\,]\\,, \\quad &\\frac{\\eta_{2j+1}(x)}{\\theta_{2j+1}(x)}=[p-j-1/j\\,]\\,\n", "\\end{aligned}\n", "\\end{equation*} \n", "with $j=0\\,,1\\,,\\dots\\,,p\\,,$ and derived the recursive relations: \n", "\\begin{equation*}\n", "\\begin{aligned} \n", "&\\frac{\\eta_{2j}(x)}{\\theta_{2j}(x)}=\\frac{[\\bar{\\eta}_{2j-1}\\,\\eta_{2j-2}(x)-x\\,\\bar{\\eta}_{2j-2}\\,\\eta_{2j-1}(x)]/\\bar{\\eta}_{2j-1}}{[\\bar{\\eta}_{2j-1}\\,\\theta_{2j-2}(x)-x\\,\\bar{\\eta}_{2j-2}\\,\\theta_{2j-1}(x)]/\\bar{\\eta}_{2j-1}}\\, \\\\ \\\\\n", "&\\frac{\\eta_{2j+1}(x)}{\\theta_{2j+1}(x)}=\\frac{[\\bar{\\eta}_{2j}\\,\\eta_{2j-1}(x)-\\bar{\\eta}_{2j-1}\\,\\eta_{2j}(x)]/(\\bar{\\eta}_{2j}-\\bar{\\eta}_{2j-1})}{[\\bar{\\eta}_{2j}\\,\\theta_{2j-1}(x)-\\bar{\\eta}_{2j-1}\\,\\theta_{2j}(x)]/(\\bar{\\eta}_{2j}-\\bar{\\eta}_{2j-1})}\\,\n", "\\end{aligned}\n", "\\end{equation*}\n", "\n", "Where $\\bar{\\eta}_{k}=$ is the coefficient of the higher power of $\\eta_{\\,k}$ polynomial; otherwise, $\\bar{\\eta}_{k}=0$. Both recursive relations are divided by a normalization constant $\\left(\\theta_{k}(0)=1\\,\\right)$. It requires the order of $p^2$ operations to derive an approximant. To start the recursive algorithm we input\n", "\\begin{equation*}\n", "\\begin{aligned}\n", "&\\eta_{0}(x)=\\sum_{n=0}^{p}a_n\\, x^n\\,, \\quad\n", "&\\theta_{0}(x)=1\\,,\\\\\n", "&\\eta_{1}(x)=\\sum_{n=0}^{p-1}a_n\\, x^n \\,, \\quad\n", "&\\theta_{1}(x)=1\\,,\n", "\\end{aligned}\n", "\\end{equation*} \n", "\n", "which follows the path in the Padé table illustrated below. The algorithm stops at $j=p\\,,$ with the calculation of $[0/p-1](x)\\,$ and $[0/p](x)\\,.$\n", "\n", "$$\\begin{array}{| l |l|l|l|l|l|}\n", "\t\\hline & & & & \\ \\ \\ \\ & \\ \\ \\\\\n", "\t\\hline & & & & &\\\\\n", "\t\\hline & & 5 \\rightarrow& 6 \\uparrow && \\\\\n", "\t\\hline & 3 \\rightarrow& 4 \\uparrow & & &\\\\\n", "\t\\hline 1 \\rightarrow& 2 \\uparrow & & & &\\\\\n", "\t\\hline 0 \\uparrow & & & & &\\\\\n", "\t\\hline\n", "\\end{array}$$\n", "\n", "For large values of $p$ the method is clearly more efficient than the presented direct algorithm. Although, the direct algorithm might be better choice when the Padé table is seriously non-normal [2].\n", "\n", "\n", "\n", "[1] Padé, H. (1892). Sur la représentation approchéee d'une fonction par des fractions rationnelles. Annales scientifiques\n", "de l' École normale supérieure, 9(9):3-93.\n", "\n", "[2] Baker, G. A. J. (1975). Essentials of Padé Approximants. Academic Press." ] }, { "cell_type": "markdown", "id": "e38352ab", "metadata": {}, "source": [ "## Examples" ] }, { "cell_type": "code", "execution_count": null, "id": "28869892", "metadata": {}, "outputs": [], "source": [ "import sympy as sp\n", "from padepy import baker_algorithm as ba" ] }, { "cell_type": "code", "execution_count": 8, "id": "3c3818e7", "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$\\displaystyle e^{x}$" ], "text/plain": [ "exp(x)" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "var = sp.Symbol(\"x\")\n", "func = sp.exp(var)\n", "p = 2\n", "q = 2\n", "func" ] }, { "cell_type": "code", "execution_count": 11, "id": "f2fafb2d", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Pade [2,2](x) stored at matrix index 4.\n" ] }, { "data": { "text/latex": [ "$\\displaystyle \\frac{\\frac{x^{2}}{12} + \\frac{x}{2} + 1}{\\frac{x^{2}}{12} - \\frac{x}{2} + 1}$" ], "text/plain": [ "(x**2/12 + x/2 + 1)/(x**2/12 - x/2 + 1)" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pade, path = ba.pade(p, q, var, func)\n", "pade" ] }, { "cell_type": "code", "execution_count": 12, "id": "7d43f619", "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$\\displaystyle \\left[\\begin{matrix}\\frac{x^{4}}{24} + \\frac{x^{3}}{6} + \\frac{x^{2}}{2} + x + 1\\\\\\frac{x^{3}}{6} + \\frac{x^{2}}{2} + x + 1\\\\\\frac{\\frac{x^{3}}{24} + \\frac{x^{2}}{4} + \\frac{3 x}{4} + 1}{1 - \\frac{x}{4}}\\\\\\frac{\\frac{x^{2}}{6} + \\frac{2 x}{3} + 1}{1 - \\frac{x}{3}}\\\\\\frac{\\frac{x^{2}}{12} + \\frac{x}{2} + 1}{\\frac{x^{2}}{12} - \\frac{x}{2} + 1}\\end{matrix}\\right]$" ], "text/plain": [ "Matrix([\n", "[ x**4/24 + x**3/6 + x**2/2 + x + 1],\n", "[ x**3/6 + x**2/2 + x + 1],\n", "[(x**3/24 + x**2/4 + 3*x/4 + 1)/(1 - x/4)],\n", "[ (x**2/6 + 2*x/3 + 1)/(1 - x/3)],\n", "[ (x**2/12 + x/2 + 1)/(x**2/12 - x/2 + 1)]])" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "path" ] }, { "cell_type": "code", "execution_count": 25, "id": "f8f70450", "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$\\displaystyle \\frac{z + 1}{\\sqrt{z^{2} + 1}}$" ], "text/plain": [ "(z + 1)/sqrt(z**2 + 1)" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "var = sp.Symbol(\"z\")\n", "func = (var + 1) / sp.sqrt(var**2 + 1)\n", "p = 3\n", "q = 5\n", "func" ] }, { "cell_type": "code", "execution_count": 28, "id": "fae86f08", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Pade [3,5](z) stored at matrix index 10.\n" ] }, { "data": { "text/latex": [ "$\\displaystyle \\frac{\\frac{147 z^{3}}{272} + \\frac{5 z^{2}}{8} + \\frac{147 z}{136} + 1}{- \\frac{z^{5}}{64} + \\frac{41 z^{4}}{272} + \\frac{5 z^{3}}{136} + \\frac{71 z^{2}}{68} + \\frac{11 z}{136} + 1}$" ], "text/plain": [ "(147*z**3/272 + 5*z**2/8 + 147*z/136 + 1)/(-z**5/64 + 41*z**4/272 + 5*z**3/136 + 71*z**2/68 + 11*z/136 + 1)" ] }, "execution_count": 28, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pade, full_path = ba.pade(p, q, var, func, 0, False)\n", "pade" ] }, { "cell_type": "code", "execution_count": 29, "id": "46db1dd4", "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$\\displaystyle \\left[\\begin{matrix}\\frac{35 z^{8}}{128} - \\frac{5 z^{7}}{16} - \\frac{5 z^{6}}{16} + \\frac{3 z^{5}}{8} + \\frac{3 z^{4}}{8} - \\frac{z^{3}}{2} - \\frac{z^{2}}{2} + z + 1\\\\- \\frac{5 z^{7}}{16} - \\frac{5 z^{6}}{16} + \\frac{3 z^{5}}{8} + \\frac{3 z^{4}}{8} - \\frac{z^{3}}{2} - \\frac{z^{2}}{2} + z + 1\\\\\\frac{- \\frac{75 z^{7}}{128} + \\frac{z^{6}}{64} + \\frac{45 z^{5}}{64} - \\frac{z^{4}}{16} - \\frac{15 z^{3}}{16} + \\frac{3 z^{2}}{8} + \\frac{15 z}{8} + 1}{\\frac{7 z}{8} + 1}\\\\\\frac{- \\frac{11 z^{6}}{16} + \\frac{7 z^{4}}{8} - \\frac{3 z^{2}}{2} + 1}{1 - z}\\\\\\frac{\\frac{z^{6}}{64} - \\frac{15 z^{5}}{352} - \\frac{z^{4}}{16} + \\frac{15 z^{3}}{44} + \\frac{3 z^{2}}{8} + \\frac{45 z}{44} + 1}{\\frac{75 z^{2}}{88} + \\frac{z}{44} + 1}\\\\\\frac{- \\frac{z^{5}}{24} - \\frac{z^{4}}{24} + \\frac{z^{3}}{3} + \\frac{z^{2}}{3} + z + 1}{\\frac{5 z^{2}}{6} + 1}\\\\\\frac{- \\frac{41 z^{5}}{704} + \\frac{z^{4}}{16} + \\frac{41 z^{3}}{88} + \\frac{3 z^{2}}{4} + \\frac{123 z}{88} + 1}{\\frac{5 z^{3}}{16} + \\frac{75 z^{2}}{88} + \\frac{35 z}{88} + 1}\\\\\\frac{- \\frac{17 z^{4}}{56} - \\frac{5 z^{2}}{7} + 1}{- \\frac{11 z^{3}}{14} + \\frac{11 z^{2}}{14} - z + 1}\\\\\\frac{\\frac{z^{4}}{16} + \\frac{41 z^{3}}{68} + \\frac{3 z^{2}}{4} + \\frac{41 z}{34} + 1}{\\frac{41 z^{4}}{272} + \\frac{11 z^{3}}{68} + \\frac{71 z^{2}}{68} + \\frac{7 z}{34} + 1}\\\\\\frac{\\frac{z^{3}}{2} + \\frac{z^{2}}{2} + z + 1}{\\frac{z^{4}}{8} + z^{2} + 1}\\\\\\frac{\\frac{147 z^{3}}{272} + \\frac{5 z^{2}}{8} + \\frac{147 z}{136} + 1}{- \\frac{z^{5}}{64} + \\frac{41 z^{4}}{272} + \\frac{5 z^{3}}{136} + \\frac{71 z^{2}}{68} + \\frac{11 z}{136} + 1}\\\\\\frac{1 - \\frac{23 z^{2}}{22}}{\\frac{17 z^{5}}{88} - \\frac{17 z^{4}}{88} - \\frac{5 z^{3}}{11} + \\frac{5 z^{2}}{11} - z + 1}\\\\\\frac{\\frac{5 z^{2}}{8} + \\frac{147 z}{92} + 1}{\\frac{147 z^{6}}{1472} - \\frac{85 z^{5}}{736} - \\frac{31 z^{4}}{368} + \\frac{25 z^{3}}{92} + \\frac{97 z^{2}}{184} + \\frac{55 z}{92} + 1}\\\\\\frac{z + 1}{\\frac{z^{6}}{16} - \\frac{z^{4}}{8} + \\frac{z^{2}}{2} + 1}\\\\\\frac{\\frac{179 z}{184} + 1}{- \\frac{5 z^{7}}{128} + \\frac{147 z^{6}}{1472} - \\frac{55 z^{5}}{1472} - \\frac{31 z^{4}}{368} - \\frac{15 z^{3}}{368} + \\frac{97 z^{2}}{184} - \\frac{5 z}{184} + 1}\\\\\\frac{1}{- \\frac{23 z^{7}}{16} + \\frac{23 z^{6}}{16} - \\frac{11 z^{5}}{8} + \\frac{11 z^{4}}{8} - \\frac{3 z^{3}}{2} + \\frac{3 z^{2}}{2} - z + 1}\\\\\\frac{1}{\\frac{179 z^{8}}{128} - \\frac{23 z^{7}}{16} + \\frac{23 z^{6}}{16} - \\frac{11 z^{5}}{8} + \\frac{11 z^{4}}{8} - \\frac{3 z^{3}}{2} + \\frac{3 z^{2}}{2} - z + 1}\\end{matrix}\\right]$" ], "text/plain": [ "Matrix([\n", "[ 35*z**8/128 - 5*z**7/16 - 5*z**6/16 + 3*z**5/8 + 3*z**4/8 - z**3/2 - z**2/2 + z + 1],\n", "[ -5*z**7/16 - 5*z**6/16 + 3*z**5/8 + 3*z**4/8 - z**3/2 - z**2/2 + z + 1],\n", "[ (-75*z**7/128 + z**6/64 + 45*z**5/64 - z**4/16 - 15*z**3/16 + 3*z**2/8 + 15*z/8 + 1)/(7*z/8 + 1)],\n", "[ (-11*z**6/16 + 7*z**4/8 - 3*z**2/2 + 1)/(1 - z)],\n", "[ (z**6/64 - 15*z**5/352 - z**4/16 + 15*z**3/44 + 3*z**2/8 + 45*z/44 + 1)/(75*z**2/88 + z/44 + 1)],\n", "[ (-z**5/24 - z**4/24 + z**3/3 + z**2/3 + z + 1)/(5*z**2/6 + 1)],\n", "[ (-41*z**5/704 + z**4/16 + 41*z**3/88 + 3*z**2/4 + 123*z/88 + 1)/(5*z**3/16 + 75*z**2/88 + 35*z/88 + 1)],\n", "[ (-17*z**4/56 - 5*z**2/7 + 1)/(-11*z**3/14 + 11*z**2/14 - z + 1)],\n", "[ (z**4/16 + 41*z**3/68 + 3*z**2/4 + 41*z/34 + 1)/(41*z**4/272 + 11*z**3/68 + 71*z**2/68 + 7*z/34 + 1)],\n", "[ (z**3/2 + z**2/2 + z + 1)/(z**4/8 + z**2 + 1)],\n", "[ (147*z**3/272 + 5*z**2/8 + 147*z/136 + 1)/(-z**5/64 + 41*z**4/272 + 5*z**3/136 + 71*z**2/68 + 11*z/136 + 1)],\n", "[ (1 - 23*z**2/22)/(17*z**5/88 - 17*z**4/88 - 5*z**3/11 + 5*z**2/11 - z + 1)],\n", "[ (5*z**2/8 + 147*z/92 + 1)/(147*z**6/1472 - 85*z**5/736 - 31*z**4/368 + 25*z**3/92 + 97*z**2/184 + 55*z/92 + 1)],\n", "[ (z + 1)/(z**6/16 - z**4/8 + z**2/2 + 1)],\n", "[(179*z/184 + 1)/(-5*z**7/128 + 147*z**6/1472 - 55*z**5/1472 - 31*z**4/368 - 15*z**3/368 + 97*z**2/184 - 5*z/184 + 1)],\n", "[ 1/(-23*z**7/16 + 23*z**6/16 - 11*z**5/8 + 11*z**4/8 - 3*z**3/2 + 3*z**2/2 - z + 1)],\n", "[ 1/(179*z**8/128 - 23*z**7/16 + 23*z**6/16 - 11*z**5/8 + 11*z**4/8 - 3*z**3/2 + 3*z**2/2 - z + 1)]])" ] }, "execution_count": 29, "metadata": {}, "output_type": "execute_result" } ], "source": [ "full_path" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3.9.5 64-bit", "language": "python", "name": "python3" }, "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.9.5" }, "vscode": { "interpreter": { "hash": "62ffe33b4b70e70e2ad7eabb350bd852d57ace81263fde7960940ad59587f0df" } } }, "nbformat": 4, "nbformat_minor": 5 }