{ "cells": [ { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# The latest version of ore_algebra should be used, run the following command if you need it\n", "# sage -pip install git+https://github.com/mkauers/ore_algebra/" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# Create the differential algebra to encode linear differential equation\n", "# (and shift algebra which expresses linear recurrence for coefficients)\n", "from ore_algebra import *\n", "Pols. = PolynomialRing(QQ)\n", "Diff. = OreAlgebra(Pols)\n", "Ind. = PolynomialRing(QQ)\n", "Shift. = OreAlgebra(Ind)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "(16*z^4 - z^2)*Dz^3 + (128*z^3 + 8*z^2 - 6*z)*Dz^2 + (224*z^2 + 28*z - 6)*Dz + 64*z + 12" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Bostan and Kauers compute the following differential equation for walks on {N,S,E,W} in N^2\n", "annDiff = z^2*(4*z-1)*(4*z+1)*Dz^3+2*z*(4*z+1)*(16*z-3)*Dz^2+(2*(112*z^2+14*z-3))*Dz+64*z+12\n", "show(annDiff)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1, 2, 6, 18, 60, 200, 700, 2450, 8820, 31752]\n" ] }, { "data": { "text/plain": [ "0: A005566: Number of walks of length n on square lattice, starting at origin, staying in first quadrant." ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# The counting sequence satisfies a recurrence relation\n", "# We verify the first terms match what we expect and check with OEIS\n", "rec = annDiff.to_S(Shift)\n", "LST = rec.to_list([1,2],10)\n", "print LST\n", "oeis(LST)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "scrolled": false }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "[4^n*n^(-1)*(1 - 3/2*n^(-1) + O(n^(-2))),\n", " (-4)^n*n^(-3)*(1 - 9/2*n^(-1) + O(n^(-2)))]" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Asymptotically, the counting sequence is a C-linear combination of the following\n", "show(rec.generalized_series_solutions(n=2))" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[z^(-2) - 4*log(z)/z - 8*log(z) - 16*z*log(z) - 8*z, 1/z, 1 + 2*z]" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# There is a basis of solutions of the C-vectors space satisfying \n", "# annDiff with the following expansions at the origin, note the \n", "# final basis element is the generating function\n", "annDiff.local_basis_expansions(0, order = 4)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "([-1.2732395447 +/- 3.84e-11] + [+/- 5.69e-26]*I)*4^N/N" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# We will see how to use numerical analytic continuation to\n", "# numerically approximate the constants appearing in asymptotics of a_n\n", "var('N')\n", "(annDiff.numerical_transition_matrix([0,1/4],1e-10) * vector([0,0,1]))[0] * 4^N/N" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "#############################################\n", "# Now, we switch to walks on {S,NW,NE} in N^2\n", "#############################################\n", "# Bostan and Kauers guessed the following differential equation\n", "annDiff2 = 35389440*z^11-33030144*z^10-9289728*z^9-23841792*z^8-27162624*z^7-12344832*z^6-4327680*z^5-1108416*z^4-251136*z^3-48576*z^2-7632*z-288+z^3*(3*z-1)*(z+1)*(8*z^2-1)*(1536*z^7-1536*z^6-584*z^4-472*z^3-90*z^2-22*z-3)*(8*z^2+1)*Dz^5+(24*(7372800*z^13-5111808*z^12-3861504*z^11-3342528*z^10-5235328*z^9-2030496*z^8-353552*z^7-83716*z^6+12956*z^5+7144*z^4+2172*z^3+473*z^2+30*z-6))*Dz^2+(24*(7372800*z^12-5996544*z^11-2856960*z^10-4204416*z^9-5463040*z^8-2341824*z^7-629536*z^6-159376*z^5-19344*z^4-1714*z^3+810*z^2+291*z+36))*Dz+4*z*(14745600*z^13-8454144*z^12-9897984*z^11-4799616*z^10-9929216*z^9-3347904*z^8-151520*z^7-28552*z^6+86040*z^5+29732*z^4+3856*z^3+786*z^2-68*z-39)*Dz^3+2*z^2*(3686400*z^13-1671168*z^12-3059712*z^11-697344*z^10-2319872*z^9-654144*z^8+100864*z^7+19304*z^6+34644*z^5+9932*z^4-245*z^3-21*z^2-86*z-21)*Dz^4" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1, 1, 2, 3, 8, 15, 39, 77, 216, 459, 1265, 2739, 7842, 17641, 49854]\n" ] }, { "data": { "text/plain": [ "0: A151255: Number of walks within N^2 (the first quadrant of Z^2) starting at (0,0) and consisting of n steps taken from {(-1, -1), (-1, 1), (1, 0)}" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# The counting sequence satisfies a recurrence relation\n", "# We verify the first terms match what we expect and check with OEIS\n", "# Note the step set listed is rotated along y=x, but this is the same counting sequence\n", "rec2 = annDiff2.to_S(Shift)\n", "LST2 = rec2.to_list([1,1,2,3,8,15,39,77,216,459,1265,2739,7842],15)\n", "print LST2\n", "oeis(LST2)" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(-3*n^5 - 207*n^4 - 5703*n^3 - 78417*n^2 - 538110*n - 1474200)*Sn^13 + (-16*n^5 - 972*n^4 - 23416*n^3 - 279204*n^2 - 1644520*n - 3817632)*Sn^12 + (-37*n^5 - 1707*n^4 - 28237*n^3 - 181701*n^2 - 125590*n + 2053896)*Sn^11 + (-226*n^5 - 9530*n^4 - 144746*n^3 - 872614*n^2 - 852804*n + 6679728)*Sn^10 + (822*n^5 + 48634*n^4 + 1113518*n^3 + 12413318*n^2 + 67631692*n + 144363120)*Sn^9 + (3608*n^5 + 177528*n^4 + 3426488*n^3 + 32398152*n^2 + 149754224*n + 269735040)*Sn^8 + (2584*n^5 + 103208*n^4 + 1368248*n^3 + 5882104*n^2 - 10416432*n - 97310592)*Sn^7 + (19072*n^5 + 583168*n^4 + 5981184*n^3 + 17490944*n^2 - 64394368*n - 335913984)*Sn^6 + (-38784*n^5 - 1890048*n^4 - 35004288*n^3 - 311046144*n^2 - 1336300032*n - 2233147392)*Sn^5 + (-169984*n^5 - 6339584*n^4 - 92063744*n^3 - 653989888*n^2 - 2279353344*n - 3123947520)*Sn^4 + (-13824*n^5 - 1463808*n^4 - 27635712*n^3 - 210683904*n^2 - 721477632*n - 923074560)*Sn^3 + (-294912*n^5 - 6119424*n^4 - 50356224*n^3 - 205332480*n^2 - 414720000*n - 331776000)*Sn^2 + (-98304*n^5 - 2850816*n^4 - 27623424*n^3 - 119832576*n^2 - 238878720*n - 176947200)*Sn + 294912*n^5 + 4423680*n^4 + 25067520*n^3 + 66355200*n^2 + 80805888*n + 35389440\n" ] } ], "source": [ "# rec2 is a recurrence of order 13\n", "print rec2" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(119*n^8 + 5653*n^7 + 115861*n^6 + 1336183*n^5 + 9467140*n^4 + 42108436*n^3 + 114517520*n^2 + 173523456*n + 111670272)*Sn^6 + (-238*n^8 - 11306*n^7 - 230998*n^6 - 2647830*n^5 - 18589124*n^4 - 81634936*n^3 - 218192288*n^2 - 322785312*n - 200721024)*Sn^5 + (-357*n^8 - 18149*n^7 - 390129*n^6 - 4658007*n^5 - 33913714*n^4 - 154581604*n^3 - 431577656*n^2 - 675732384*n - 454791168)*Sn^4 + (3808*n^7 + 140912*n^6 + 2190688*n^5 + 18541136*n^4 + 92218096*n^3 + 269303408*n^2 + 426900960*n + 282623616)*Sn^3 + (-7616*n^8 - 279920*n^7 - 4369312*n^6 - 37687488*n^5 - 195457760*n^4 - 619368400*n^3 - 1156475776*n^2 - 1136454720*n - 426998016)*Sn^2 + (15232*n^8 + 571264*n^7 + 9101696*n^6 + 80378752*n^5 + 429980864*n^4 + 1425914944*n^3 + 2861835520*n^2 + 3178066176*n + 1495249920)*Sn + 22848*n^8 + 811200*n^7 + 12110976*n^6 + 98863488*n^5 + 480082752*n^4 + 1410943680*n^3 + 2430407424*n^2 + 2216884992*n + 804879360\n" ] } ], "source": [ "# In fact, b_n satisfies a lower order recurrence than the one defined by rec2\n", "rec3 = 22848*n^8+811200*n^7+12110976*n^6+98863488*n^5+480082752*n^4+1410943680*n^3+2430407424*n^2+2216884992*n+804879360+(64*(n+3))*(2+n)*(238*n^6+7736*n^5+102106*n^4+698972*n^3+2610955*n^2+5031314*n+3893880)*Sn-(16*(n+3))*(476*n^7+16067*n^6+224881*n^5+1680825*n^4+7173635*n^3+17189620*n^2+20710876*n+8895792)*Sn^2+(16*(238*n^7+8807*n^6+136918*n^5+1158821*n^4+5763631*n^3+16831463*n^2+26681310*n+17663976))*Sn^3-(n+6)*(357*n^7+16007*n^6+294087*n^5+2893485*n^4+16552804*n^3+55264780*n^2+99988976*n+75798528)*Sn^4-(2*(n+7))*(119*n^7+4820*n^6+81759*n^5+751602*n^4+4033348*n^3+12584032*n^2+21007920*n+14337216)*Sn^5+(n+4)*(n+7)*(119*n^4+2440*n^3+18089*n^2+56772*n+62316)*(n+8)^2*Sn^6\n", "print rec3" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1, 1, 2, 3, 8, 15, 39, 77, 216, 459, 1265, 2739, 7842, 17641, 49854]\n" ] }, { "data": { "text/plain": [ "(n^5 + 69*n^4 + 1901*n^3 + 26139*n^2 + 179370*n + 491400)*Sn^13 + (16/3*n^5 + 324*n^4 + 23416/3*n^3 + 93068*n^2 + 1644520/3*n + 1272544)*Sn^12 + (37/3*n^5 + 569*n^4 + 28237/3*n^3 + 60567*n^2 + 125590/3*n - 684632)*Sn^11 + (226/3*n^5 + 9530/3*n^4 + 144746/3*n^3 + 872614/3*n^2 + 284268*n - 2226576)*Sn^10 + (-274*n^5 - 48634/3*n^4 - 1113518/3*n^3 - 12413318/3*n^2 - 67631692/3*n - 48121040)*Sn^9 + (-3608/3*n^5 - 59176*n^4 - 3426488/3*n^3 - 10799384*n^2 - 149754224/3*n - 89911680)*Sn^8 + (-2584/3*n^5 - 103208/3*n^4 - 1368248/3*n^3 - 5882104/3*n^2 + 3472144*n + 32436864)*Sn^7 + (-19072/3*n^5 - 583168/3*n^4 - 1993728*n^3 - 17490944/3*n^2 + 64394368/3*n + 111971328)*Sn^6 + (12928*n^5 + 630016*n^4 + 11668096*n^3 + 103682048*n^2 + 445433344*n + 744382464)*Sn^5 + (169984/3*n^5 + 6339584/3*n^4 + 92063744/3*n^3 + 653989888/3*n^2 + 759784448*n + 1041315840)*Sn^4 + (4608*n^5 + 487936*n^4 + 9211904*n^3 + 70227968*n^2 + 240492544*n + 307691520)*Sn^3 + (98304*n^5 + 2039808*n^4 + 16785408*n^3 + 68444160*n^2 + 138240000*n + 110592000)*Sn^2 + (32768*n^5 + 950272*n^4 + 9207808*n^3 + 39944192*n^2 + 79626240*n + 58982400)*Sn - 98304*n^5 - 1474560*n^4 - 8355840*n^3 - 22118400*n^2 - 26935296*n - 11796480" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# To see that b_n satisfies rec3 note that there is a a solution B_n of rec3 that equals b_n\n", "# for the first 100 terms, and that b_n - B_n satisfies the following recurrence of order 13\n", "print rec3.to_list([1,1,2,3,8,15],15)\n", "rec3.lclm(rec2)" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "scrolled": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "3^n*n^(-1/2)*(1 + O(n^(-1)))\n", "(2.828427124746190?)^n*n^(-2)*(1 + O(n^(-1)))\n", "(0.?e-57 + 2.828427124746190?*I)^n*n^(-3)*(1 + O(n^(-1)))\n", "(0.?e-57 - 2.828427124746190?*I)^n*n^(-3)*(1 + O(n^(-1)))\n", "(-1)^n*n^(-3/2)*(1 + O(n^(-1)))\n", "(-2.828427124746190?)^n*n^(-2)*(1 + O(n^(-1)))\n" ] } ], "source": [ "# This, b_n has an asymptotic expansion as a C-linear combination of the following\n", "gs = rec3.generalized_series_solutions(n=1)\n", "for k in gs:\n", " print(k)" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2*sqrt(2)" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# The possible exponential growth factors are 3, -1, and 2\\sqrt{2} times a fourth root of unity\n", "gs[1].exponential_part().radical_expression()" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[log(z)/z^2 - 4*log(z)/z - 8*log(z) - 10*z*log(z) - 22*z^2*log(z) + 1/4*z^2, z^(-2), 1/z, 1 - z^2, z + 3*z^2]\n" ] } ], "source": [ "# There is a basis of solutions of the C-vectors space satisfying \n", "# annDiff with the following expansions at the origin, note the sum of\n", "# the last 2 basis elements give the generating function\n", "print annDiff2.local_basis_expansions(0, order = 5)\n", "ini = [0,0,0,1,1]" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[+/- 3.68e-106] + [+/- 1.43e-100]*I" ] }, "execution_count": 27, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# The constant lambda_1 appearing in front of the 3^n/sqrt(n) is numerically zero to many decimal places\n", "(annDiff2.numerical_transition_matrix([0,1/3],1e-100) * vector(ini))[0]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "SageMath 8.4", "language": "", "name": "sagemath" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.15" } }, "nbformat": 4, "nbformat_minor": 1 }