{ "cells": [ { "cell_type": "code", "execution_count": 1, "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": 2, "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": 3, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "(n^2 + 2*n + 1)*Sn - 16*n^2 - 16*n - 4" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Recurrence relation for binomial(2n,n)^2\n", "# Sn denotes the shift operator n -> n+1\n", "rec = (n+1)^2*Sn - 4*(2*n+1)^2\n", "show(rec)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1, 4, 36, 400, 4900, 63504, 853776, 11778624, 165636900, 2363904400]\n" ] }, { "data": { "text/plain": [ "0: A002894: a(n) = binomial(2n, n)^2." ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# We verify the first terms match what we expect and check with OEIS\n", "LST = rec.to_list([1,4],10)\n", "print LST\n", "oeis(LST)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "[16^n*n^(-1)*(1 - 1/4*n^(-1) + 1/32*n^(-2) + 1/128*n^(-3) - 5/2048*n^(-4) + O(n^(-5)))]" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "show(rec.generalized_series_solutions())" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "(-16*z^3 + z^2)*Dz^3 + (-80*z^2 + 3*z)*Dz^2 + (-68*z + 1)*Dz - 4" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# The ore_algebra package can convert the recurrence into a differential equation\n", "annDiff = rec.to_D(Diff)\n", "show(annDiff)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "(-16*z^2 + z)*Dz^2 + (-32*z + 1)*Dz - 4" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "" ], "text/plain": [ "(z^2 - 1/16*z)*Dz^2 + (2*z - 1/16)*Dz + 1/4" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# It is not the minimal order differential equation\n", "# It has the following annihilating differential equation as a right factor\n", "annDiff2 = z*(1-16*z)*Dz^2+(1-32*z)*Dz-4\n", "show(annDiff2)\n", "show(annDiff.gcrd(annDiff2))" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(-16) * (z - 1/16) * z" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Potential singularities at z=0 and z=1/16\n", "annDiff2.leading_coefficient().factor()" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "[log(z) + 4*z*log(z) + 8*z + 36*z^2*log(z) + 84*z^2 + 400*z^3*log(z) + 2960/3*z^3,\n", " 1 + 4*z + 36*z^2 + 400*z^3]" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Two solutions as the origin, one with logarithmic singularity\n", "show(annDiff2.local_basis_expansions(0, order = 4))" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "[log(z - 1/16) - 4*1/16*(16*z - 1)*log(z - 1/16) - 8*z - 1/16 + 36*1/256*(16*z - 1)^2*log(z - 1/16) + 84*1/256*(16*z - 1)^2,\n", " 1 - 4*z - 1/16 + 36*1/256*(16*z - 1)^2]" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Similar at z=1/16\n", "show(annDiff2.local_basis_expansions(1/16, order = 3))" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "# At the origin, the power series element is the generating function\n", "# So the generating function is expressed in the basis with coordinates\n", "ini = [0,1]" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "[ [+/- 5.98e-9] + [+/- 6.32e-10]*I [-0.318 +/- 3.10e-4] + [+/- 2.93e-10]*I]\n", "[ [-3.14 +/- 1.60e-3] + [+/- 1.89e-8]*I [+/- 8.78e-9] + [1.00 +/- 3e-8]*I]" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# This is the change of basis matrix (we take a low tolerance to fit on screen)\n", "show(annDiff2.numerical_transition_matrix([0,1/16],1e-2))" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "([-0.31830988618379067 +/- 2.19e-18] + [+/- 3.04e-38]*I, [+/- 2.05e-19] + [1.0000000000000000 +/- 7.09e-18]*I)" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "" ], "text/plain": [ "[log(z - 1/16) - 4*1/16*(16*z - 1)*log(z - 1/16) - 8*z - 1/16 + 36*1/256*(16*z - 1)^2*log(z - 1/16) + 84*1/256*(16*z - 1)^2,\n", " 1 - 4*z - 1/16 + 36*1/256*(16*z - 1)^2]" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "CBmatrix = annDiff2.numerical_transition_matrix([0,1/16],1e-1);\n", "bas = annDiff2.numerical_transition_matrix([0,1/16]) * vector(ini)\n", "loc = annDiff2.local_basis_expansions(1/16,3)\n", "show(bas)\n", "show(loc)" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[-0.31830988618379067 +/- 2.19e-18] + [+/- 3.04e-38]*I" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# The extraction [z^n]log(1-16z) = -16^n/n then implies\n", "bas[0]" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.318309886183791" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# The connection constant is 1/pi\n", "(1/pi).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 }