{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# The latest version of ore_algebra is needed\n", "# git clone http://marc.mezzarobba.net/code/ore_algebra-analytic.git/" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "from ore_algebra import *\n", "Pols. = PolynomialRing(QQ)\n", "Diff. = OreAlgebra(Pols)\n", "\n", "Ind. = PolynomialRing(QQ)\n", "Shift. = OreAlgebra(Ind)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "annDiff = w^2*(256*w-81)*(16777216*w^2-5971968*w+216513)*(8388608*w^2+12773376*w-2381643)*Dw^4+3*w*(126100789566373888*w^5+145702882866364416*w^4-144989442627600384*w^3+37699919662546944*w^2-3508882176945024*w+69613650565965)*Dw^3+(12*(83879543059775488*w^5+127612480615612416*w^4-91463559904493568*w^3+18580583017193472*w^2-1333811087342316*w+13922730113193))*Dw^2+(72*(9147936743096320*w^4+18313981068312576*w^3-9089964406996992*w^2+1351389647038464*w-58491439678683))*Dw+50665495807918080*w^3+141873283866820608*w^2-42412899154526208*w+4304856990449664" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[1,\n", " 136/27,\n", " 12296/243,\n", " 123829568/177147,\n", " 18061677320/1594323,\n", " 25936424350144/129140163,\n", " 118279908284579776/31381059609,\n", " 20782338563398585856/282429536481,\n", " 3759288229575285225368/2541865828329,\n", " 1520786375427958515905082560/50031545098999707]" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "rec = annDiff.to_S(Shift)\n", "\n", "# Verify initial terms\n", "rec.to_list([1, 136/27, 12296/243, 123829568/177147, 18061677320/1594323, 25936424350144/129140163],10)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[(24.40775175129386?)^n*n^(-3/2)*(1 + O(n^(-1))),\n", " (5.954754895325534?)^n*n^(-5)*(1 + O(n^(-1))),\n", " (3.174739831197724?)^n*n^(-3/2)*(1 + O(n^(-1))),\n", " (256/81)^n*n^(-3/2)*(1 + O(n^(-1))),\n", " (-0.5914926431743928?)^n*n^(-5)*(1 + O(n^(-1)))]" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Asymptotics are C-linear combination of\n", "rec.generalized_series_solutions(1)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "(36028797018963968) * (w - 81/256) * w^2 * (w^2 - 729/2048*w + 216513/16777216) * (w^2 + 6237/4096*w - 2381643/8388608)" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "" ], "text/plain": [ "[-1.690638102670645?,\n", " 0,\n", " 0.04097059041691498?,\n", " 0.1679330245456446?,\n", " 0.3149864408330850?,\n", " 0.31640625000000000?]" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Get singularities of solutions to annDiff\n", "lc = annDiff.leading_coefficient()\n", "show(lc.factor())\n", "rts = annDiff.leading_coefficient().roots(QQbar, multiplicities=False); \n", "show(rts)" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "[1/2*log(w)^2, log(w), 1, w]" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "" ], "text/plain": [ "[1/2*log(w)^2 + 68/27*w*log(w)^2 + 248/27*w*log(w),\n", " log(w) + 136/27*w*log(w),\n", " 1,\n", " w]" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Generating function starts 1 + 136/27*w + O(w^2)\n", "# Set initial conditions for our solution in expansion at origin\n", "show(annDiff.local_basis_monomials(0))\n", "show(annDiff.local_basis_expansions(0, order = 2))\n", "ini = [0,0,1,136/27]" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "# Analytically continue to first singularity, and find new basis\n", "bas = annDiff.numerical_transition_matrix([0,rts[2]]) * vector(ini)\n", "loc = annDiff.local_basis_expansions(rts[2],1)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "[1, sqrt(w - 0.04097059041691498?), 0, 0]" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "" ], "text/plain": [ "[+/- 1.71e-26] + [4.864657531666138 +/- 1.13e-16]*I" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Only one singular basis element\n", "show(loc)\n", "show(bas[1])" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "4.86465753166614*I*sqrt(w - 0.04097059041691498?)" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Asymptotic contribution of this singularity is coefficient of z^n in\n", "C1 = 4.864657531666138*I;\n", "show(C1 * sqrt(w-rts[2]))" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "-1/2/sqrt(pi)*(1/zeta)^n*n^(-3/2) + O((1/zeta)^n*n^(-5/2))" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Asymptotics of [w^n] sqrt(1-w/zeta)\n", "var('zeta,n')\n", "asm = asymptotic_expansions.SingularityAnalysis('n', alpha=-1/2, zeta=zeta, precision=1)\n", "show(asm)\n", "asm = -1/2/sqrt(pi)*(1/zeta)^n*n^(-3/2)" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "0.492332365958782*24.40775175129386?^n/(sqrt(pi)*n^(3/2))" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# This matches our computationally obtained asymptotics \n", "ASM1 = C1*sqrt(-rts[2])*asm.subs({zeta:rts[2]})\n", "show(ASM1)" ] }, { "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 }