{ "cells": [ { "cell_type": "code", "execution_count": null, "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": 5, "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": 9, "metadata": {}, "outputs": [], "source": [ "pp = z^2-1" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [], "source": [ "annDiff = Dz - 1" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [], "source": [ "annDiff.numerical_solution?" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[2.71828182845904523536028747135266249775724709369995957496696762772407663035354759457138217852516642742746639193200305992181741359662904357290033429526059563073813232862794349076323382988075319525101901157383418793070215408914993488416750924476146066808226480016847741185374234544243710753907774499206955170276183860626133138458300075204493382656029760673711320070932870912744374704723069697720931014169283681902551510865746377211125238978442505695369677078544996996794686445490598793163688923009879312773617821542499922957635148220826989519366803318252886939849646510582093923982948879332036250944311730123819706841614039701983767932068328237646480429531180232878250981945581530175671736133206981125099618188159304169035159888851934580727386673858942287922849989208680582574927961048419844436346324496848756023362482704197862320900216099023530436994184914631409343173814364054625315209618369088870701676839642437814059271456354906130310720851038375051011574770417189861068739696552126715468895703503540 +/- 3.40e-1002]" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "annDiff.numerical_solution(ini=[1], path=[0, 1], eps=1e-1000)" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2.7182818284590452353602874714" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Numeric approximation of exp(1)\n", "exp(1).n(100)" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [], "source": [ "annDiff = -2*(1-z)*Dz-1" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "(2*z - 2)*Dz - 1" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "show(annDiff)" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "[1 - 1/2*z - 1/8*z^2 - 1/16*z^3 - 5/128*z^4 + O(z^5)]" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "" ], "text/plain": [ "1 + (-1/2)*z + (-1/8)*z^2 + (-1/16)*z^3 + (-5/128)*z^4 + Order(z^5)" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "show(annDiff.generalized_series_solutions())\n", "show(sqrt(1-z).series(z,5))" ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[0.70710678119 +/- 4.25e-12]" ] }, "execution_count": 33, "metadata": {}, "output_type": "execute_result" } ], "source": [ "annDiff.numerical_solution(ini=[1], path=[0, 1/2], eps=1e-10)" ] }, { "cell_type": "code", "execution_count": 34, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0" ] }, "execution_count": 34, "metadata": {}, "output_type": "execute_result" } ], "source": [ "annDiff.numerical_solution(ini=[1], path=[0, 1], eps=1e-10)" ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [ { "ename": "ValueError", "evalue": "Step 0 --> 2 passes through or too close to singular point 1 (to compute the connection to a singular point, make it a vertex of the path)", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mValueError\u001b[0m Traceback (most recent call last)", "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mannDiff\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnumerical_solution\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mini\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mInteger\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mpath\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mInteger\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mInteger\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0meps\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mRealNumber\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'1e-10'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m", "\u001b[0;32m/Applications/SageMath-8.4.app/Contents/Resources/sage/local/lib/python2.7/site-packages/ore_algebra/ore_operator_1_1.pyc\u001b[0m in \u001b[0;36mnumerical_solution\u001b[0;34m(self, ini, path, eps, post_transform, **kwds)\u001b[0m\n\u001b[1;32m 2519\u001b[0m post_mat = matrix(1, dop.order(),\n\u001b[1;32m 2520\u001b[0m lambda i, j: ZZ(j).factorial()*post_transform[j])\n\u001b[0;32m-> 2521\u001b[0;31m \u001b[0mctx\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mancont\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mContext\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdop\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mpath\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0meps\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwds\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2522\u001b[0m \u001b[0mpairs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mancont\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0manalytic_continuation\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mctx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mini\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mini\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mpost\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mpost_mat\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2523\u001b[0m \u001b[0;32massert\u001b[0m \u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpairs\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;32m/Applications/SageMath-8.4.app/Contents/Resources/sage/local/lib/python2.7/site-packages/ore_algebra/analytic/analytic_continuation.pyc\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, dop, path, eps, keep, algorithm, assume_analytic)\u001b[0m\n\u001b[1;32m 36\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpath\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0minitial_path\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mPath\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpath\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdop\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 37\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0massume_analytic\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 38\u001b[0;31m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0minitial_path\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcheck_singularity\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 39\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mall\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mis_regular\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mx\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpath\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mvert\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 40\u001b[0m raise NotImplementedError(\"analytic continuation through irregular \"\n", "\u001b[0;32m/Applications/SageMath-8.4.app/Contents/Resources/sage/local/lib/python2.7/site-packages/ore_algebra/analytic/path.pyc\u001b[0m in \u001b[0;36mcheck_singularity\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m 733\u001b[0m \"\"\"\n\u001b[1;32m 734\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mstep\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 735\u001b[0;31m \u001b[0mstep\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcheck_singularity\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 736\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 737\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mcheck_convergence\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;32m/Applications/SageMath-8.4.app/Contents/Resources/sage/local/lib/python2.7/site-packages/ore_algebra/analytic/path.pyc\u001b[0m in \u001b[0;36mcheck_singularity\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m 568\u001b[0m \u001b[0;34m\"Step {} passes through or too close to singular point{} {} \"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 569\u001b[0m \u001b[0;34m\"(to compute the connection to a singular point, make it \"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 570\u001b[0;31m \"a vertex of the path)\".format(self, plural, sings))\n\u001b[0m\u001b[1;32m 571\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 572\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mcheck_convergence\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;31mValueError\u001b[0m: Step 0 --> 2 passes through or too close to singular point 1 (to compute the connection to a singular point, make it a vertex of the path)" ] } ], "source": [ "annDiff.numerical_solution(ini=[1], path=[0, 2], eps=1e-10)" ] }, { "cell_type": "code", "execution_count": 36, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[+/- 1.84e-14] + [-1.0000000000 +/- 2.46e-13]*I" ] }, "execution_count": 36, "metadata": {}, "output_type": "execute_result" } ], "source": [ "annDiff.numerical_solution(ini=[1], path=[0, I, 2], eps=1e-10)" ] }, { "cell_type": "code", "execution_count": 38, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[+/- 8.65e-12] + [-1.0000000000 +/- 8.68e-12]*I" ] }, "execution_count": 38, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Same as semi-circle with z = 1-e^(-i*theta) for theta from 0 to pi\n", "annDiff.numerical_solution(ini=[1], path=[0, 1+I, 2], eps=1e-10)" ] }, { "cell_type": "code", "execution_count": 39, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[+/- 1.84e-14] + [1.0000000000 +/- 2.46e-13]*I" ] }, "execution_count": 39, "metadata": {}, "output_type": "execute_result" } ], "source": [ "annDiff.numerical_solution(ini=[1], path=[0, -I, 2], eps=1e-10)" ] }, { "cell_type": "code", "execution_count": 42, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[-1.0000000000 +/- 4.37e-13] + [+/- 3.37e-13]*I" ] }, "execution_count": 42, "metadata": {}, "output_type": "execute_result" } ], "source": [ "annDiff.numerical_solution(ini=[1], path=[0, I, 2, -I, 0], eps=1e-10)" ] }, { "cell_type": "code", "execution_count": 43, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[1.0000000000 +/- 4.60e-12] + [+/- 4.57e-12]*I" ] }, "execution_count": 43, "metadata": {}, "output_type": "execute_result" } ], "source": [ "annDiff.numerical_solution(ini=[1], path=[0, I, 2, -I, 0, I, 2, -I, 0], eps=1e-10)" ] }, { "cell_type": "code", "execution_count": 44, "metadata": {}, "outputs": [], "source": [ "annDiff = (1-z)*Dz^2-Dz" ] }, { "cell_type": "code", "execution_count": 45, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[log(z - 1), 1]\n", "1*z + 1/2*z^2 + Order(z^3)\n" ] } ], "source": [ "print annDiff.local_basis_monomials(1)\n", "print log(1/(1-z)).series(z,3)" ] }, { "cell_type": "code", "execution_count": 50, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[+/- 4.01e-12] + [6.2831853072 +/- 2.51e-11]*I\n" ] } ], "source": [ "print annDiff.numerical_solution(ini=[0,1], path=[0, I, 2, -I, 0], eps=1e-10)" ] }, { "cell_type": "code", "execution_count": 51, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[+/- 2.97e-10] + [12.566370614 +/- 6.55e-10]*I\n", "[+/- 5.80e-9] + [18.84955592 +/- 7.34e-9]*I\n" ] } ], "source": [ "print annDiff.numerical_solution(ini=[0,1], path=[0, I, 2, -I, 0, I, 2, -I, 0], eps=1e-10)\n", "print annDiff.numerical_solution(ini=[0,1], path=[0, I, 2, -I, 0, I, 2, -I, 0, I, 2, -I, 0], eps=1e-10)" ] }, { "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 }