{ "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*(81*w^2+14*w+1)*Dw^3+3*w*(162*w^2+21*w+1)*Dw^2+(21*w+1)*(27*w+1)*Dw+81*w+3" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[1, -3, 9, -3, -279, 2997, -19431, 65853, 292329, -7202523]" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "rec = annDiff.to_S(Shift)\n", "\n", "# Verify initial terms\n", "rec.to_list([1,-3],10)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[(-7.000000000000000? + 5.656854249492380?*I)^n*n^(-3/2)*(1 + O(n^(-1))),\n", " (-7.000000000000000? - 5.656854249492380?*I)^n*n^(-3/2)*(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": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(81\\right) \\cdot w^{2} \\cdot (w^{2} + \\frac{14}{81} w + \\frac{1}{81})$ " ], "text/plain": [ "(81) * w^2 * (w^2 + 14/81*w + 1/81)" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left[0, -0.0864197530864198? - 0.06983770678385654? \\sqrt{-1}, -0.0864197530864198? + 0.06983770678385654? \\sqrt{-1}\\right]$ " ], "text/plain": [ "[0,\n", " -0.0864197530864198? - 0.06983770678385654?*I,\n", " -0.0864197530864198? + 0.06983770678385654?*I]" ] }, "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": 7, "metadata": {}, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left[\\frac{1}{2} \\, \\log\\left(w\\right)^{2}, \\log\\left(w\\right), 1\\right]$ " ], "text/plain": [ "[1/2*log(w)^2, log(w), 1]" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left[\\frac{1}{2}\\log\\left(w\\right)^{2} - \\frac{3}{2}w \\log\\left(w\\right)^{2} - 4w \\log\\left(w\\right), \\log\\left(w\\right) - 3w \\log\\left(w\\right) - 4w, 1 - 3w\\right]$ " ], "text/plain": [ "[1/2*log(w)^2 - 3/2*w*log(w)^2 - 4*w*log(w),\n", " log(w) - 3*w*log(w) - 4*w,\n", " 1 - 3*w]" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Generating function starts 1 - 3w + 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]" ] }, { "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[1]]) * vector(ini)\n", "loc = annDiff.local_basis_expansions(rts[1],1)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left[1, \\sqrt{w + 0.0864197530864198? + 0.06983770678385654? \\sqrt{-1}}, 0\\right]$ " ], "text/plain": [ "[1, sqrt(w + 0.0864197530864198? + 0.06983770678385654?*I), 0]" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\verb|[-3.5933098558743233|\\phantom{\\verb!x!}\\verb|+/-|\\phantom{\\verb!x!}\\verb|1.95e-17]|\\phantom{\\verb!x!}\\verb|+|\\phantom{\\verb!x!}\\verb|[-0.38132214909311386|\\phantom{\\verb!x!}\\verb|+/-|\\phantom{\\verb!x!}\\verb|2.47e-18]*I|$ " ], "text/plain": [ "[-3.5933098558743233 +/- 1.95e-17] + [-0.38132214909311386 +/- 2.47e-18]*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": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}-\\left(3.59330985587432 + 0.381322149093114i\\right) \\, \\sqrt{w + 0.0864197530864198? + 0.06983770678385654? \\sqrt{-1}}$ " ], "text/plain": [ "-(3.59330985587432 + 0.381322149093114*I)*sqrt(w + 0.0864197530864198? + 0.06983770678385654?*I)" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Asymptotic contribution of this singularity is coefficient of z^n in\n", "C1 = (-3.5933098558743233 - 0.3813221490931139*I)\n", "show(C1 * sqrt(w-rts[1]))" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}-\\frac{1}{2 \\, \\sqrt{\\pi}} \\left(\\frac{1}{\\zeta}\\right)^{n} n^{-\\frac{3}{2}} + O\\!\\left(\\left(\\frac{1}{\\zeta}\\right)^{n} n^{-\\frac{5}{2}}\\right)$ " ], "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": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\frac{\\left(0.543449606382202 + 0.259547320313100i\\right) \\, \\left(-7.000000000000000? + 5.656854249492380? \\sqrt{-1}\\right)^{n}}{\\sqrt{\\pi} n^{\\frac{3}{2}}}$ " ], "text/plain": [ "(0.543449606382202 + 0.259547320313100*I)*(-7.000000000000000? + 5.656854249492380?*I)^n/(sqrt(pi)*n^(3/2))" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "ASM1 = C1*sqrt(-rts[1])*asm.subs({zeta:rts[1]})\n", "show(ASM1)" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(0.543449606382202 + 0.259547320313100*I)*(-7.000000000000000? + 5.656854249492380?*I)^n/(sqrt(pi)*n^(3/2))" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ASM1" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left[1, \\sqrt{w + 0.0864197530864198? - 0.06983770678385654? \\sqrt{-1}}, 0\\right]$ " ], "text/plain": [ "[1, sqrt(w + 0.0864197530864198? - 0.06983770678385654?*I), 0]" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\verb|[-3.5933098558743233|\\phantom{\\verb!x!}\\verb|+/-|\\phantom{\\verb!x!}\\verb|1.95e-17]|\\phantom{\\verb!x!}\\verb|+|\\phantom{\\verb!x!}\\verb|[0.38132214909311386|\\phantom{\\verb!x!}\\verb|+/-|\\phantom{\\verb!x!}\\verb|2.47e-18]*I|$ " ], "text/plain": [ "[-3.5933098558743233 +/- 1.95e-17] + [0.38132214909311386 +/- 2.47e-18]*I" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Repeat for second singular point\n", "bas = annDiff.numerical_transition_matrix([0,rts[2]]) * vector(ini)\n", "loc = annDiff.local_basis_expansions(rts[2],1)\n", "show(loc); show(bas[1])" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\frac{\\left(0.543449606382202 - 0.259547320313100i\\right) \\, \\left(-7.000000000000000? - 5.656854249492380? \\sqrt{-1}\\right)^{n}}{\\sqrt{\\pi} n^{\\frac{3}{2}}}$ " ], "text/plain": [ "(0.543449606382202 - 0.259547320313100*I)*(-7.000000000000000? - 5.656854249492380?*I)^n/(sqrt(pi)*n^(3/2))" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "C2 = (-3.5933098558743233 + 0.3813221490931139*I)\n", "ASM2 = C2*sqrt(-rts[2])*asm.subs({zeta:rts[2]})\n", "show(ASM2)" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\frac{\\left(0.543449606382202 + 0.259547320313100i\\right) \\, \\left(-7.000000000000000? + 5.656854249492380? \\sqrt{-1}\\right)^{n}}{\\sqrt{\\pi} n^{\\frac{3}{2}}} + \\frac{\\left(0.543449606382202 - 0.259547320313100i\\right) \\, \\left(-7.000000000000000? - 5.656854249492380? \\sqrt{-1}\\right)^{n}}{\\sqrt{\\pi} n^{\\frac{3}{2}}}$ " ], "text/plain": [ "(0.543449606382202 + 0.259547320313100*I)*(-7.000000000000000? + 5.656854249492380?*I)^n/(sqrt(pi)*n^(3/2)) + (0.543449606382202 - 0.259547320313100*I)*(-7.000000000000000? - 5.656854249492380?*I)^n/(sqrt(pi)*n^(3/2))" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "show(ASM1 + ASM2)" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "-0.351180716688941" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# This matches our computationally obtained asymptotics \n", "ASM = ASM1 + ASM2\n", "(ASM/9^n*n^(3/2)).subs({n:10000}).n()" ] }, { "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 }