{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## Example 7.4 (Two Critical Points with Positive Coordinates)\n", "We run our algorithm and compute the Kronecker representation \"by hand\" using Sage's capability to compute Gröbner bases. \n", "*Requirements: None*" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# We will need our procedure to get the Hessian matrix appearing in asymptotics\n", "# This code copied from Chapter 5\n", "def getHes(H,vars):\n", " dd = len(vars)\n", " V = zero_vector(SR,dd)\n", " U = matrix(SR,dd)\n", " M = matrix(SR,dd-1)\n", "\n", " for j in range(dd):\n", " V[j] = 1\n", " for i in range(dd):\n", " U[i,j] = vars[i]*vars[j]*diff(H,vars[i],vars[j])/vars[-1]/diff(H,vars[-1])\n", " for i in range(dd-1):\n", " for j in range(dd-1):\n", " M[i,j] = V[i]*V[j] + U[i,j] - V[j]*U[i,-1] - V[i]*U[j,-1] + V[i]*V[j]*U[-1,-1]\n", " if i == j: M[i,j] = M[i,j] + V[i]\n", " return M" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[(x + 40*y - 20)*(x + y - 1) - 1,\n", " (2*x + 41*y - 21)*x - λ,\n", " (41*x + 80*y - 60)*y - λ,\n", " (t*x + 40*t*y - 20)*(t*x + t*y - 1) - 1]" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Set basic info for this example\n", "var('x,y,λ,t,u')\n", "vars = [x,y,λ,t]\n", "H = (1-x-y)*(20-x-40*y)-1\n", "\n", "# Form the extended critical point system\n", "C = [H] + [v * diff(H,v) - λ for v in vars[0:-2]]\n", "S = C + [H.subs([v==t*v for v in vars])]\n", "S" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "u^8 - 198726679/6563778*u^7 + 173956613387/511974684*u^6 - 73067829295337/39934025352*u^5 + 8519230060835315/1557426988728*u^4 - 2508167123647177/259571164788*u^3 + 655159497358069/64892791197*u^2 - 158493265458143/27323280504*u + 170155274912099/119802076056" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Define the multivariate polynomial ring and compute a Gröbner basis\n", "R = PolynomialRing(QQbar,5,'x,y,λ,t,u',order='lex')\n", "S = S + [u - (x+t)]\n", "GB = R.ideal([R(k) for k in S]).groebner_basis()\n", "\n", "# The Kronecker is representation parametrized by the roots of the following polynomial\n", "P = QQ[u](GB[-1])\n", "P" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# The polynomial P is square-free: its gcd with its derivative is zero\n", "# (Note: we need to cast between polynomial and symbolic rings to compute desired quantities)\n", "QQ[u](P).gcd(QQ[u](SR(P).diff(u)))" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\verb|In|\\phantom{\\verb!x!}\\verb|the|\\phantom{\\verb!x!}\\verb|Kronecker|\\phantom{\\verb!x!}\\verb|rep.|\\phantom{\\verb!x!}\\verb|the|\\phantom{\\verb!x!}\\verb|variable|\\phantom{\\verb!x!}\\verb|\t| x \\phantom{\\verb!x!}\\verb|is|\\phantom{\\verb!x!}\\verb|represented|\\phantom{\\verb!x!}\\verb|by\t| \\frac{899}{39} u^{7} - \\frac{123578457923}{255987342} u^{6} + \\frac{1313142566253}{369759494} u^{5} - \\frac{10209620955321937}{778713494364} u^{4} + \\frac{590457682083167}{21630930399} u^{3} - \\frac{335415361956118}{10246230189} u^{2} + \\frac{2753820801687041}{129785582394} u - \\frac{1502824746811243}{259571164788} \\phantom{\\verb!x!}\\verb|divided|\\phantom{\\verb!x!}\\verb|by\t| 8 u^{7} - \\frac{1391086753}{6563778} u^{6} + \\frac{173956613387}{85329114} u^{5} - \\frac{365339146476685}{39934025352} u^{4} + \\frac{8519230060835315}{389356747182} u^{3} - \\frac{2508167123647177}{86523721596} u^{2} + \\frac{1310318994716138}{64892791197} u - \\frac{158493265458143}{27323280504}$ " ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\verb|In|\\phantom{\\verb!x!}\\verb|the|\\phantom{\\verb!x!}\\verb|Kronecker|\\phantom{\\verb!x!}\\verb|rep.|\\phantom{\\verb!x!}\\verb|the|\\phantom{\\verb!x!}\\verb|variable|\\phantom{\\verb!x!}\\verb|\t| x \\phantom{\\verb!x!}\\verb|is|\\phantom{\\verb!x!}\\verb|represented|\\phantom{\\verb!x!}\\verb|by\t| \\frac{899}{39} u^{7} - \\frac{123578457923}{255987342} u^{6} + \\frac{1313142566253}{369759494} u^{5} - \\frac{10209620955321937}{778713494364} u^{4} + \\frac{590457682083167}{21630930399} u^{3} - \\frac{335415361956118}{10246230189} u^{2} + \\frac{2753820801687041}{129785582394} u - \\frac{1502824746811243}{259571164788} \\phantom{\\verb!x!}\\verb|divided|\\phantom{\\verb!x!}\\verb|by\t| 8 u^{7} - \\frac{1391086753}{6563778} u^{6} + \\frac{173956613387}{85329114} u^{5} - \\frac{365339146476685}{39934025352} u^{4} + \\frac{8519230060835315}{389356747182} u^{3} - \\frac{2508167123647177}{86523721596} u^{2} + \\frac{1310318994716138}{64892791197} u - \\frac{158493265458143}{27323280504}\n", "\\end{math}" ], "text/plain": [ "'In the Kronecker rep. the variable \\t' x ' is represented by\\t' 899/39*u^7 - 123578457923/255987342*u^6 + 1313142566253/369759494*u^5 - 10209620955321937/778713494364*u^4 + 590457682083167/21630930399*u^3 - 335415361956118/10246230189*u^2 + 2753820801687041/129785582394*u - 1502824746811243/259571164788 ' divided by\\t' 8*u^7 - 1391086753/6563778*u^6 + 173956613387/85329114*u^5 - 365339146476685/39934025352*u^4 + 8519230060835315/389356747182*u^3 - 2508167123647177/86523721596*u^2 + 1310318994716138/64892791197*u - 158493265458143/27323280504" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Compute the Q polynomials for the Kronecker representation from the Gröbner basis\n", "Pd = diff(P,u)\n", "Q = {}\n", "for v in vars:\n", " [Qv] = list(filter(lambda p: SR(p).has(v), GB))\n", " _,Q[v] = QQ[u](Pd * v.subs(solve(SR(Qv),v))).quo_rem(P)\n", "sys = [SR(Q[v]/Pd) for v in vars]\n", "\n", "# For instance, here is the first variable in the Kronecker representation\n", "show(\"In the Kronecker rep. the variable \\t\",vars[0], \" is represented by\\t\",Q[vars[0]],\" divided by\\t\",Pd)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[7, 7, 7, 7]" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# All polynomials in the Kronecker representation have degree one less than the degree of P\n", "[Q[v].degree() for v in vars]" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[5, 4, 5, 4]" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# They have coefficients that are at most 5 digits long\n", "# As expected, this is an *order of magnitude* smaller than the constants\n", "# appearing in the Gröbner basis, which have constants with ~60 digits\n", "[Q[v].norm(infinity).log(10).round() for v in vars]" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(u^4 - 1211/78*u^3 + 172525/3042*u^2 - 242111/3042*u + 9161/234) * (u^4 - 48409909/3281889*u^3 + 27671511239/511974684*u^2 - 19048879957/255987342*u + 18573875659/511974684)" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Note that P factors\n", "P.factor()" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "u^4 - 1211/78*u^3 + 172525/3042*u^2 - 242111/3042*u + 9161/234" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# The critical points, where t=1, are given by the values of u which are roots of \n", "Pt = gcd(P,Pd-Q[t])\n", "Pt" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The root u = 1.548232473567013? gives the critical point\n", "[x,y] = [0.548232473567013?, 0.3099773361396642?] \n", "\n", "The root u = 10.997110519821022? gives the critical point\n", "[x,y] = [9.997110519821022?, 0.2527749731647175?] \n", "\n", "The root u = 1.490149016126495? - 0.2807916065035742?*I gives the critical point\n", "[x,y] = [0.4901490161264949? - 0.2807916065035742?*I, 0.5808033325272964? + 0.1623547816338436?*I] \n", "\n", "The root u = 1.490149016126495? + 0.2807916065035742?*I gives the critical point\n", "[x,y] = [0.4901490161264949? + 0.2807916065035742?*I, 0.5808033325272964? - 0.1623547816338436?*I] \n", "\n" ] } ], "source": [ "# Now there are two critical point with positive coordinates\n", "for uval in Pt.roots(QQbar, multiplicities=False):\n", " print(\"The root u = {} gives the critical point\".format(uval))\n", " print(\"[x,y] =\",[p.subs(u == uval) for p in sys[:-2]],\"\\n\")" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "At the solutions of the system, t takes the values [1.000000000000000?, 1.709936749104775?, 0.09218565523266332?, 1.000000000000000?, 0.7114300338237230? - 0.10463158979488784?*I, 0.7114300338237230? + 0.10463158979488784?*I, 1.000000000000000? + 0.?e-27*I, 1.000000000000000? + 0.?e-27*I]\n", "Here are the values lying in (0,1): [0.09218565523266332?]\n" ] } ], "source": [ "# We can compute the algebraic values of t at solutions of the system\n", "tvals = [sys[-1].subs(u == uval) for uval in P.roots(QQbar, multiplicities=False)]\n", "print(\"At the solutions of the system, t takes the values {}\".format(tvals))\n", "\n", "# Because one of the positive points is not minimal, there *is* a solution with t in (0,1)\n", "num_in_zo = list(filter(lambda c: c.is_real() and c>0 and c<1, tvals))\n", "print(\"Here are the values lying in (0,1): {}\".format(num_in_zo))" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The ray from the origin to the critical point cp1 = [0.548232473567013?, 0.3099773361396642?] contains 1.000000000000000?*cp1 and 1.709936749104775?*cp1 \n", "\n", "The ray from the origin to the critical point cp2 = [9.997110519821022?, 0.2527749731647175?] contains 0.09218565523266332?*cp2 and 1.000000000000000?*cp2\n" ] } ], "source": [ "# Define the critical points with positive real coordinates \n", "[u1,u2] = filter(lambda r:SR(r).is_real(), Pt.roots(QQbar, multiplicities=False))\n", "cp1 = [r.subs(u==u1) for r in sys[:2]]\n", "cp2 = [r.subs(u==u2) for r in sys[:2]]\n", "\n", "# Find the solutions of the extended system which represent the same critical points\n", "alg_sys = [[p.subs(u=uval) for p in sys] for uval in P.roots(QQbar, multiplicities=False)]\n", "cp1t = list(filter(lambda pt: cp1 == pt[:2], alg_sys))\n", "cp2t = list(filter(lambda pt: cp2 == pt[:2], alg_sys))\n", "\n", "# Print results\n", "print(\"The ray from the origin to the critical point cp1 = {} contains {}*cp1 and {}*cp1 \\n\".format(cp1,cp1t[0][-1],cp1t[1][-1]))\n", "print(\"The ray from the origin to the critical point cp2 = {} contains {}*cp2 and {}*cp2\".format(cp2,cp2t[0][-1],cp2t[1][-1]))" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[0.548232473567013?,\n", " 9.997110519821022?,\n", " 0.5648805044366720?,\n", " 0.5648805044366720?]" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# No other critical point has the same coordinate-wise modulus as the minimal critical point\n", "[abs(sys[0].subs(u==uval)) for uval in Pt.roots(QQbar, multiplicities=False)]" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The Kronecker representation for the critical points (where t=1 in the extended system) is given by\n" ] }, { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left[x = \\frac{35061 \\, u^{3} - 203363 \\, u^{2} + 381283 \\, u - 234261}{12168 \\, u^{3} - 141687 \\, u^{2} + 345050 \\, u - 242111}, y = \\frac{20982 \\, u^{3} - 263923 \\, u^{2} + 644528 \\, u - 446047}{4 \\, {\\left(12168 \\, u^{3} - 141687 \\, u^{2} + 345050 \\, u - 242111\\right)}}\\right]$ " ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left[x = \\frac{35061 \\, u^{3} - 203363 \\, u^{2} + 381283 \\, u - 234261}{12168 \\, u^{3} - 141687 \\, u^{2} + 345050 \\, u - 242111}, y = \\frac{20982 \\, u^{3} - 263923 \\, u^{2} + 644528 \\, u - 446047}{4 \\, {\\left(12168 \\, u^{3} - 141687 \\, u^{2} + 345050 \\, u - 242111\\right)}}\\right]\n", "\\end{math}" ], "text/plain": [ "[x == (35061*u^3 - 203363*u^2 + 381283*u - 234261)/(12168*u^3 - 141687*u^2 + 345050*u - 242111),\n", " y == 1/4*(20982*u^3 - 263923*u^2 + 644528*u - 446047)/(12168*u^3 - 141687*u^2 + 345050*u - 242111)]" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "where u satisfies u^4 - 1211/78*u^3 + 172525/3042*u^2 - 242111/3042*u + 9161/234=0.\n" ] } ], "source": [ "# We can simplify the system to the points where t = 1\n", "# First, compute reciprocal of the new derivative modulo the new P\n", "newP = Pt\n", "newPd = diff(newP,u)\n", "_,invPd,_ = Pd.xgcd(newP)\n", "\n", "# Then update the Q polynomials using polynomial division\n", "newQ = {}\n", "for v in vars[0:-1]:\n", " _,newQ[v] = QQ[u](Q[v] * newPd * invPd).quo_rem(newP)\n", "\n", "# Print the new representation\n", "print(\"The Kronecker representation for the critical points (where t=1 in the extended system) is given by\")\n", "show([v == newQ[v]/newPd for v in vars[0:-2]])\n", "print(\"where u satisfies {}=0.\".format(newP))" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The Kronecker representation for the critical points (where t=1 in the extended system) is given by\n" ] }, { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left[T = \\frac{796653 \\, u^{3} - 5489818 \\, u^{2} + 11297717 \\, u - 7320092}{78 \\, {\\left(12168 \\, u^{3} - 141687 \\, u^{2} + 345050 \\, u - 242111\\right)}}, h = \\frac{21576789 \\, u^{3} - 92722412 \\, u^{2} + 118817977 \\, u - 40527214}{39 \\, {\\left(12168 \\, u^{3} - 141687 \\, u^{2} + 345050 \\, u - 242111\\right)}}, g = \\left(-1\\right)\\right]$ " ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left[T = \\frac{796653 \\, u^{3} - 5489818 \\, u^{2} + 11297717 \\, u - 7320092}{78 \\, {\\left(12168 \\, u^{3} - 141687 \\, u^{2} + 345050 \\, u - 242111\\right)}}, h = \\frac{21576789 \\, u^{3} - 92722412 \\, u^{2} + 118817977 \\, u - 40527214}{39 \\, {\\left(12168 \\, u^{3} - 141687 \\, u^{2} + 345050 \\, u - 242111\\right)}}, g = \\left(-1\\right)\\right]\n", "\\end{math}" ], "text/plain": [ "[T == 1/78*(796653*u^3 - 5489818*u^2 + 11297717*u - 7320092)/(12168*u^3 - 141687*u^2 + 345050*u - 242111),\n", " h == 1/39*(21576789*u^3 - 92722412*u^2 + 118817977*u - 40527214)/(12168*u^3 - 141687*u^2 + 345050*u - 242111),\n", " g == -1]" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "where u satisfies u^4 - 1211/78*u^3 + 172525/3042*u^2 - 242111/3042*u + 9161/234=0.\n" ] } ], "source": [ "# Compute the Hessian determinant Htilde for final asymptotic formula\n", "hes = vars[-3]*diff(H,vars[-3])*getHes(H,vars[0:-2])\n", "hesdet = hes.determinant().factor()\n", "\n", "# Add the Hessian, the product of the variables, and the negation of the numerator\n", "# to the Kronecker rep. We use a straightforward but inefficient approach by \n", "# computing another Gröbner basis!\n", "var('T,h,g')\n", "PolySys = [newP, T - prod(vars[:-2]), h - hesdet, g+1] + S\n", "R2 = PolynomialRing(QQ,8,'g,T,h,x,y,λ,t,u',order='lex')\n", "GB2 = R2.ideal([R2(str(k)) for k in PolySys]).groebner_basis()\n", "for v in [T,h,g]:\n", " [Qv] = list(filter(lambda p: SR(p).has(v), GB2))\n", " _,newQ[v] = QQ[u](newPd * v.subs(solve(SR(Qv),v))).quo_rem(newP)\n", " \n", "# Print the new elements of the representation\n", "print(\"The Kronecker representation for the critical points (where t=1 in the extended system) is given by\")\n", "show([v == newQ[v]/newPd for v in [T,h,g]])\n", "print(\"where u satisfies {}=0.\".format(newP))" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Asymptotics are given by evaluating\n" ] }, { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}-\\frac{\\sqrt{39} {\\left(12168 \\, u^{3} - 141687 \\, u^{2} + 345050 \\, u - 242111\\right)} \\left(\\frac{78 \\, {\\left(12168 \\, u^{3} - 141687 \\, u^{2} + 345050 \\, u - 242111\\right)}}{796653 \\, u^{3} - 5489818 \\, u^{2} + 11297717 \\, u - 7320092}\\right)^{n} \\sqrt{\\frac{588627 \\, u^{3} - 2845362 \\, u^{2} + 4562963 \\, u - 2409968}{21576789 \\, u^{3} - 92722412 \\, u^{2} + 118817977 \\, u - 40527214}}}{{\\left(588627 \\, u^{3} - 2845362 \\, u^{2} + 4562963 \\, u - 2409968\\right)} \\sqrt{\\pi n}}$ " ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}[1]{\\mathbf{#1}}-\\frac{\\sqrt{39} {\\left(12168 \\, u^{3} - 141687 \\, u^{2} + 345050 \\, u - 242111\\right)} \\left(\\frac{78 \\, {\\left(12168 \\, u^{3} - 141687 \\, u^{2} + 345050 \\, u - 242111\\right)}}{796653 \\, u^{3} - 5489818 \\, u^{2} + 11297717 \\, u - 7320092}\\right)^{n} \\sqrt{\\frac{588627 \\, u^{3} - 2845362 \\, u^{2} + 4562963 \\, u - 2409968}{21576789 \\, u^{3} - 92722412 \\, u^{2} + 118817977 \\, u - 40527214}}}{{\\left(588627 \\, u^{3} - 2845362 \\, u^{2} + 4562963 \\, u - 2409968\\right)} \\sqrt{\\pi n}}\n", "\\end{math}" ], "text/plain": [ "-sqrt(39)*(12168*u^3 - 141687*u^2 + 345050*u - 242111)*(78*(12168*u^3 - 141687*u^2 + 345050*u - 242111)/(796653*u^3 - 5489818*u^2 + 11297717*u - 7320092))^n*sqrt((588627*u^3 - 2845362*u^2 + 4562963*u - 2409968)/(21576789*u^3 - 92722412*u^2 + 118817977*u - 40527214))/((588627*u^3 - 2845362*u^2 + 4562963*u - 2409968)*sqrt(pi*n))" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "At the root 1.548232473567013? of the polynomial u^4 - 1211/78*u^3 + 172525/3042*u^2 - 242111/3042*u + 9161/234.\n", "This gives dominant asymptotic behaviour\n" ] }, { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\frac{0.0545997615256984 \\cdot 5.884442204019508?^{n}}{\\sqrt{n}}$ " ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\frac{0.0545997615256984 \\cdot 5.884442204019508?^{n}}{\\sqrt{n}}\n", "\\end{math}" ], "text/plain": [ "0.0545997615256984*5.884442204019508?^n/sqrt(n)" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Get final asymptotic constants\n", "d = 2\n", "A = newQ[g]/newQ[λ]\n", "B = newQ[λ]^(d-1)/newQ[h]*newPd^(2-d)\n", "C = newPd/newQ[T]\n", "\n", "# Find root of P that gives critical point with positive coordinates\n", "[posu] = filter(lambda uval: sys[0].subs(u == uval)==cp1[0], newP.roots(QQbar, multiplicities=False))\n", "\n", "# Get result of algorithm\n", "var('n')\n", "ASM = ((n*2*pi)^((1-d)/2) * A * sqrt(B) * SR(C)^n).simplify()\n", "\n", "# Parse Sage output to get numerical approximation of asymptotics\n", "def numeval(e):\n", " E = SR(e)\n", " if E.is_constant():\n", " return E.n()\n", " else:\n", " return E\n", "numASM = prod([numeval(k) for k in ASM.subs(u=posu).expand().operands()])\n", "\n", "# Print result\n", "print(\"Asymptotics are given by evaluating\")\n", "show(ASM)\n", "print(\"At the root {} of the polynomial {}.\".format(posu,newP))\n", "print(\"This gives dominant asymptotic behaviour\")\n", "show(numASM)" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "scrolled": true }, "outputs": [], "source": [ "# Because u is the root of a quadric polynomial, we could express asymptotics in radicals\n", "# Those brave enough to view the output can uncomment and run the next line of code\n", "# print(ASM.subs(u=posu.radical_expression()).simplify_full())" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "SageMath 9.2", "language": "sage", "name": "sagemath" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.5" } }, "nbformat": 4, "nbformat_minor": 2 }