{ "cells": [ { "cell_type": "code", "execution_count": 1, "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": 2, "metadata": {}, "outputs": [], "source": [ "# Write a function to comptute the number of walks resticted to the quadrant and ending anywhere\n", "@CachedFunction \n", "def WalksIJ(i,j,n):\n", " if i<0 or j<0:\n", " return 0\n", " elif n==0 and i==0 and j==0:\n", " return 1\n", " elif n==0:\n", " return 0\n", " else:\n", " return WalksIJ(i-1,j,n-1) + WalksIJ(i+1,j,n-1) + WalksIJ(i,j-1,n-1) + WalksIJ(i,j+1,n-1)\n", " \n", "def Walks(n): \n", " return sum([WalksIJ(i,j,n) for i in range(n+1) for j in range(n+1)])" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1, 2, 6, 18, 60, 200, 700, 2450, 8820, 31752]\n" ] } ], "source": [ "LST = [Walks(k) for k in range(20)]\n", "print LST[0:10]" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "# Look the sequence up in OEIS\n", "OE = oeis(LST)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0: A005566: Number of walks of length n on square lattice, starting at origin, staying in first quadrant.\n", "0: a(n) is the number of involutions of length 2n which are invariant under the reverse-complement map and have no decreasing subsequences of length 5. - _Eric S. Egge_, Oct 21 2008\n" ] } ], "source": [ "print OE\n", "print OE[0].comments()" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(16*t^4 - t^2)*Dt^3 + (128*t^3 + 8*t^2 - 6*t)*Dt^2 + (224*t^2 + 28*t - 6)*Dt + 64*t + 12\n" ] } ], "source": [ "# Guess an annihilating differential equation (need about 30 terms)\n", "LST = [Walks(k) for k in range(40)]\n", "diffWalk = guess(LST,Diff)\n", "print diffWalk" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(4*t - 1) * (4*t + 1) * t^2" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Potential singularities at +/- 1/4\n", "diffWalk.leading_coefficient().factor()" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "[t^(-2) - 4*log(t)/t - 8*log(t) - 16*t*log(t) - 8*t, 1/t, 1 + 2*t]" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Get basis with expansions at 0, and get coefficients of GF in the basis\n", "show(diffWalk.local_basis_expansions(0, order = 4))\n", "ini = [0,0,1]" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "[log(t - 1/4) - 6*1/4*(4*t - 1)*log(t - 1/4) + 31*1/16*(4*t - 1)^2*log(t - 1/4),\n", " 1 - 14*1/16*(4*t - 1)^2,\n", " t - 1/4 - 15/2*1/16*(4*t - 1)^2]" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Expansions of basis at 1/4\n", "show(diffWalk.local_basis_expansions(1/4, order = 3))" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "([-1.2732395447351627 +/- 2.78e-17] + [+/- 2.45e-18]*I, [-2.3906971441245563 +/- 2.29e-17] + [4.0000000000000000 +/- 3.09e-17]*I, [12.890661954217663 +/- 2.66e-16] + [-24.000000000000000 +/- 1.23e-16]*I)" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Find connection coefficients\n", "bas = diffWalk.numerical_transition_matrix([0,1/4]) * vector(ini)\n", "show(bas)" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[1.2732395447351627 +/- 2.78e-17] + [+/- 2.45e-18]*I" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# The first basis element has leading singular term log(t-1/4) and the second is analytic\n", "# Thus the extraction [z^n]log(1-4z) = -4^n/n implies our sequence grows as 16^n/n * C where C is\n", "-bas[0]" ] }, { "cell_type": "code", "execution_count": 38, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1.27323954473516" ] }, "execution_count": 38, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Leading constant is 4/pi\n", "(4/pi).n()" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "[1 - 4*1/16*(4*t + 1)^2*log(t + 1/4),\n", " t + 1/4 + 1/16*(4*t + 1)^2*log(t + 1/4),\n", " 1/16*(4*t + 1)^2]" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Dominant singular term for expansions at t=-1/4 is (1+4t)log(1+4t)/4\n", "# Thus, asymptotically negligible\n", "show(diffWalk.local_basis_expansions(-1/4, order = 3))" ] }, { "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": 2 }