{ "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": [ "(10*z^4 - 35*z^3 + 15*z^2)*Dz^4 + (118*z^3 - 278*z^2 + 72*z)*Dz^3 + (366*z^2 - 480*z + 36)*Dz^2 + (300*z - 132)*Dz + 36" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Start with differential equation\n", "annDiff = 5*z^2*(2*z-1)*(z-3)*Dz^4 + 2*z*(59*z^2-139*z+36)*Dz^3 + 6*(61*z^2-80*z+6)*Dz^2 + 12*(25*z-11)*Dz + 36\n", "show(annDiff)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "(15*n^4 + 102*n^3 + 237*n^2 + 222*n + 72)*Sn^2 + (-35*n^4 - 208*n^3 - 445*n^2 - 404*n - 132)*Sn + 10*n^4 + 58*n^3 + 122*n^2 + 110*n + 36" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "" ], "text/plain": [ "(15*n^2 + 57*n + 36)*Sn^2 + (-35*n^2 - 103*n - 66)*Sn + 10*n^2 + 28*n + 18" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Get recurrence on coefficients\n", "rec = annDiff.to_S(Shift)\n", "show(rec)\n", "rec = rec.primitive_part()\n", "show(rec)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "[2^n*n^(-1)*(1 - n^(-1) + n^(-2) + O(n^(-3))), (1/3)^n*(1 + O(n^(-3)))]" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# The following gives an asymptotic basis for solutions of the recurrence\n", "show(rec.generalized_series_solutions(n=3))" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "1/z" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "" ], "text/plain": [ "1 - 1/2*z^2 - 17/18*z^3 - 427/270*z^4" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "" ], "text/plain": [ "z + 11/6*z^2 + 53/18*z^3 + 1291/270*z^4" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "" ], "text/plain": [ "z^(6/5) + 11/6*z^(11/5) + 187/63*z^(16/5) + 1705/351*z^(21/5) + 20515/2511*z^(26/5) + 3421/243*z^(31/5)" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "" ], "text/plain": [ "[0, f0, f1, 0]" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Four solutions as the origin, two with singularitites\n", "# A power series solution starting f0 + f1*z + ... has the following coefficients\n", "for k in annDiff.local_basis_expansions(0, order = 6):\n", " show(k)\n", "var('f0 f1')\n", "ini = [0,f0,f1,0]\n", "show(ini)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(10) * (z - 3) * (z - 1/2) * z^2" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Potential singularities at z=0, z=3 and z=1/2\n", "annDiff.leading_coefficient().factor()" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "log(z - 1/2) - 2*1/2*(2*z - 1)*log(z - 1/2) + 4*1/4*(2*z - 1)^2*log(z - 1/2) - 8*1/8*(2*z - 1)^3*log(z - 1/2)" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "" ], "text/plain": [ "1 + 8/25*1/8*(2*z - 1)^3" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "" ], "text/plain": [ "z - 1/2 + 4/25*1/8*(2*z - 1)^3" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "" ], "text/plain": [ "1/4*(2*z - 1)^2 - 2*1/8*(2*z - 1)^3" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Expansions at z=1/2 - one has singularity\n", "for k in annDiff.local_basis_expansions(1/2, order = 4):\n", " show(k)" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "[ [+/- 2.43e-14] + [+/- 3.41e-15]*I [0.50000000000 +/- 2.8e-15] + [+/- 1.45e-16]*I [-1.5000000000 +/- 9e-15] + [+/- 4.25e-16]*I [-1.3789520923 +/- 3.90e-11] + [+/- 3.89e-16]*I]\n", "[ [2.0000000000 +/- 8.3e-14] + [+/- 8.88e-14]*I [2.1465735903 +/- 2.01e-11] + [-1.5707963268 +/- 5.12e-12]*I [-2.8397207708 +/- 4.00e-11] + [4.7123889804 +/- 1.54e-11]*I [-2.7321494258 +/- 1.19e-11] + [4.3321057627 +/- 4.09e-13]*I]\n", "[ [-4.0000000000 +/- 1.06e-13] + [+/- 1.55e-13]*I [0.026852819440 +/- 5.85e-14] + [3.1415926536 +/- 1.03e-11]*I [1.3594415417 +/- 2.02e-11] + [-9.4247779608 +/- 3.07e-11]*I [1.0516521563 +/- 4.10e-11] + [-8.6642115254 +/- 8.3e-13]*I]\n", "[ [8.0000000000 +/- 1.02e-13] + [+/- 2.94e-13]*I [1.6742943611 +/- 1.99e-11] + [-6.2831853072 +/- 2.05e-11]*I [-4.4468830834 +/- 4.04e-11] + [18.849555922 +/- 4.62e-10]*I [-4.0889953255 +/- 3.78e-11] + [17.328423051 +/- 2.02e-10]*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(annDiff.numerical_transition_matrix([0,1/2],1e-10))" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "scrolled": false }, "outputs": [], "source": [ "# The coefficient of the singular term is f0/2 -3f1/2\n", "bas = annDiff.numerical_transition_matrix([0,1/2],1e-10) * vector(ini)" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "([0.50000000000 +/- 2.8e-15] + [+/- 1.45e-16]*I)*f0 + ([-1.5000000000 +/- 9e-15] + [+/- 4.25e-16]*I)*f1" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "show(bas[0])" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "1/(z - 3)" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "" ], "text/plain": [ "1" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "" ], "text/plain": [ "z - 3" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "" ], "text/plain": [ "(z - 3)^2" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Expansions at z=3 - one has singularity\n", "for k in annDiff.local_basis_expansions(3, order = 4):\n", " show(k)" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "[ [+/- 7.47] + [+/- 7.47]*I [+/- 6.03] + [+/- 1.53]*I [+/- 7.21] + [+/- 2.71]*I [+/- 6.69] + [+/- 1.75]*I]\n", "[ [+/- 9.23] + [+/- 8.89]*I [+/- 1.96] + [+/- 2.09]*I [+/- 3.63] + [+/- 4.01]*I [+/- 2.33] + [+/- 2.81]*I]\n", "[ [+/- 2.25] + [+/- 2.14]*I [+/- 0.448] + [+/- 0.524]*I [+/- 0.807] + [+/- 1.04]*I [+/- 0.516] + [+/- 0.740]*I]\n", "[[+/- 0.767] + [+/- 0.730]*I [+/- 0.153] + [+/- 0.179]*I [+/- 0.274] + [+/- 0.352]*I [+/- 0.182] + [+/- 0.251]*I]" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# This is the change of basis matrix (we take a low tolerance to fit on screen)\n", "# Note -- the path for continuation must move off the real line to AVOID the singularity at z=1/2\n", "show(annDiff.numerical_transition_matrix([0,I,3],1e-1))" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "([-4.500000000 +/- 8.07e-10] + [+/- 8.07e-10]*I)*f0 + ([4.50000000 +/- 1.88e-9] + [+/- 1.88e-9]*I)*f1" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# The coefficient of the singular term is (9/2)(f1-f0)\n", "bas = annDiff.numerical_transition_matrix([0,I,3],1e-10) * vector(ini)\n", "show(bas[0])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Seems like f_n ~ (3f1 - f0)/2 * 2^n/n unless 3f1=f0 in which case f_n ~ (3/2)(f0-f1)*3^(-n)" ] } ], "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 }