{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## Example 7.5 (A Non-Combinatorial Series Expansion)\n", "We run our algorithm and compute the Kronecker representation \"by hand\" using Sage's capability to compute Gröbner bases. \n", "*Requirements: ore_algebra (for optional heuristic check of asymptotics at end)*" ] }, { "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": [], "source": [ "# Set basic info for this example\n", "var('x,y,z,λ,t,u')\n", "H = 1 - (x+y+z) + 81/8*x*y*z" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The ideal generated by [a1 - 1/2*lambdaR - 1/2, a2 - 1/2*lambdaR - 1/2, a3 - 1/2*lambdaR - 1/2, b1 - 1/2*lambdaI, b2 - 1/2*lambdaI, b3 - 1/2*lambdaI, lambdaR^3 + 3*lambdaR^2 + 49/27*lambdaR - 27/8*lambdaI^6 - 5*lambdaI^4 - 32/27*lambdaI^2 + 49/81, lambdaR*lambdaI - 9/8*lambdaI^5 - 5/3*lambdaI^3 + 49/81*lambdaI, lambdaI^7 + 16/9*lambdaI^5 + 64/81*lambdaI^3 - 3136/19683*lambdaI] has dimension 0.\n" ] } ], "source": [ "# Because we no longer have a combinatorial example, we must split the\n", "# variables into real/imaginary components and define the extended critical point system\n", "var('a1,a2,a3,b1,b2,b3,lambdaR,lambdaI')\n", "assume(a1,'real'); assume(a2,'real'); assume(a3,'real')\n", "assume(b1,'real'); assume(b2,'real'); assume(b3,'real')\n", "Hc = H.subs(x=a1+I*b1, y=a2+I*b2, z=a3+I*b3)\n", "Hr = Hc.real()\n", "Hi = Hc.imag()\n", "CPab1 = [diff(Hr,v)*v + diff(Hr,w)*w - lambdaR for (v,w) in zip([a1,a2,a3],[b1,b2,b3])]\n", "CPab2 = [diff(Hi,v)*v + diff(Hi,w)*w - lambdaI for (v,w) in zip([a1,a2,a3],[b1,b2,b3])]\n", "\n", "# Compute a Gröbner basis for the critical point equations in the a and b variables\n", "R = PolynomialRing(QQ,8,'a1,a2,a3,b1,b2,b3,lambdaR,lambdaI',order='lex')\n", "J = R.ideal([R(k) for k in [Hr,Hi] + CPab1 + CPab2])\n", "GB = R.ideal([R(k) for k in [Hr,Hi] + CPab1 + CPab2]).groebner_basis()\n", "\n", "# Note there are a finite number of (complex) solutions! \n", "print(\"The ideal generated by {} has dimension {}.\".format(GB,J.dimension()))" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "# To determine minimality, we add the extra equations using x and y variables\n", "# First, define xj + I*yj as a root of H\n", "# Then encode that the coordinate-wise moduli of xj + I*yj is a multiple of crit pt aj + I*bj\n", "# Finally, add the extra equations derived from critical point method\n", "var('x1,x2,x3,y1,y2,y3,ν')\n", "Hrxy = Hr.subs(a1==x1,a2==x2,a3==x3,b1==y1,b2==y2,b3==y3)\n", "Hixy = Hi.subs(a1==x1,a2==x2,a3==x3,b1==y1,b2==y2,b3==y3)\n", "mod = [X^2+Y^2 - t*(A^2+B^2) for (X,Y,A,B) in zip([x1,x2,x3],[y1,y2,y3],[a1,a2,a3],[b1,b2,b3])]\n", "nuxy = [(Y - ν*X)*diff(Hrxy,X) - (X + ν*Y)*diff(Hrxy,Y) for (X,Y) in zip([x1,x2,x3],[y1,y2,y3])]\n", "\n", "# Define the ideal with all these equations\n", "S = PolynomialRing(QQ,16,'ν,lambdaR,lambdaI,x1,x2,x3,y1,y2,y3,b1,b2,b3,a1,a2,a3,t',order='lex')\n", "Jxy = S.ideal([S(k) for k in list(GB) + [Hrxy,Hixy] + mod + nuxy])" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "scrolled": false }, "outputs": [], "source": [ "# Compute the corresponding Gröbner basis. This can take several minutes.\n", "# (Note: the same Gröbner computation in Maple is almost instantaneous)\n", "GBxy = Jxy.groebner_basis()" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(1/10460353203) * (t - 3) * (t - 1) * (3*t - 1) * (t^2 + 3) * (t^2 - 3*t + 9) * (t^2 - t + 1) * (t^2 + t + 1) * (t^2 + 3*t + 3) * (3*t^2 + 1) * (3*t^2 + 3*t + 1) * (9*t^2 - 3*t + 1) * (9*t^3 + 12*t^2 + 4*t - 27) * (81*t^3 + 36*t^2 + 4*t - 9) * (81*t^6 + 108*t^5 + 108*t^4 + 534*t^3 + 340*t^2 - 108*t + 729) * (729*t^6 - 972*t^5 + 540*t^4 - 144*t^3 - 308*t^2 + 108*t + 243)" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# As expected, there are a finite number of solutions\n", "# At these solutions, the t variable is a root of the following\n", "GBxy[-1].factor()" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[0.3333333333333334?, 0.3451548307798085?]" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# There are two values between 0 and 1 -- corresponding to non-minimal critical points\n", "tPol = QQbar[t](SR(GBxy[-1]))\n", "badt = list(filter(lambda r: SR(r).is_real() and r>0 and r<1, tPol.roots(multiplicities=False)))\n", "badt" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[-2/3, -2/3, -2/3]\n", "[-2/3, -2/3, -2/3]\n", "[-2/3, -2/3, -2/3]\n", "[-2/3, -2/3, -2/3]\n", "[-2/3, -2/3, -2/3]\n" ] } ], "source": [ "# Here we list all the non-minimal critical points\n", "# In fact, there is only one non-minimal critical point\n", "# This implies all other critical points are minimal\n", "S2 = PolynomialRing(QQbar,15,'ν,lambdaR,lambdaI,x1,x2,x3,y1,y2,y3,b1,b2,b3,a1,a2,a3',order='lex')\n", "for T in badt:\n", " J = S2.ideal([S2(SR(k.subs(t=T))) for k in GBxy])\n", " V = J.variety()\n", " for v in V:\n", " print([(v[j] + QQbar(I)*v[k]).radical_expression() for (j,k) in zip([a1,a2,a3],[b1,b2,b3])])" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "u^3 - 8/27*u + 8/81" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Now we find asymptotics in terms of the minimal critical points\n", "# Form the extended critical point system\n", "vars = [x,y,z,λ]\n", "C = [H] + [v * diff(H,v) - λ for v in vars[0:-1]]\n", "\n", "# Define the multivariate polynomial ring and compute a Gröbner basis\n", "R = PolynomialRing(QQbar,6,'x,y,z,λ,t,u',order='lex')\n", "C = C + [u - x]\n", "GB = R.ideal([R(k) for k in C]).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": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1" ] }, "execution_count": 10, "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": 11, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}{\\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{\\frac{16}{81} u - \\frac{8}{81}}{u^{2} - \\frac{8}{81}}\n", "\\end{math}" ], "text/plain": [ "'In the Kronecker rep. the variable \\t' x ' is represented by\\t' (16/81*u - 8/81)/(u^2 - 8/81)" ] }, "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, \" is represented by\\t\",Q[vars]/Pd)" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The root u = -0.6666666666666667? gives the critical point\n", "[x,y] = [-0.6666666666666667?, -0.6666666666666667?] \n", "\n", "The root u = 0.3333333333333334? - 0.1924500897298753?*I gives the critical point\n", "[x,y] = [0.3333333333333334? - 0.1924500897298753?*I, 0.3333333333333334? - 0.1924500897298753?*I] \n", "\n", "The root u = 0.3333333333333334? + 0.1924500897298753?*I gives the critical point\n", "[x,y] = [0.3333333333333334? + 0.1924500897298753?*I, 0.3333333333333334? + 0.1924500897298753?*I] \n", "\n" ] } ], "source": [ "# There are three minimal critical points: the non-minimal point (-2/3,-2/3,-2/3)\n", "# and two minimal critical points\n", "rts = P.roots(QQbar, multiplicities=False)\n", "for uval in rts:\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": 13, "metadata": {}, "outputs": [], "source": [ "# Select the minimal critical points \n", "[u1,u2] = filter(lambda U: not sys.subs(u==U)== -2/3, P.roots(QQbar, multiplicities=False))\n", "\n", "# Compute the Hessian determinant Htilde for final asymptotic formula\n", "hes = vars[-2]*diff(H,vars[-2])*getHes(H,vars[0:-1])\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 = [P, T - prod(vars[:-1]), h - hesdet, g+1] + C\n", "R2 = PolynomialRing(QQ,8,'g,T,h,x,y,z,λ,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", " _,Q[v] = QQ[u](Pd * v.subs(solve(SR(Qv),v))).quo_rem(P)" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Asymptotics are given by adding the evaluations of\n" ] }, { "data": { "text/html": [ "" ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}{\\mathbf{#1}}\\frac{\\sqrt{\\frac{1}{6}} {\\left(81 \\, u^{2} - 8\\right)} \\left(-\\frac{81 \\, {\\left(81 \\, u^{2} - 8\\right)}}{8 \\, {\\left(81 \\, u^{2} - 48 \\, u + 16\\right)}}\\right)^{n} \\sqrt{\\frac{6561 \\, u^{4} - 5184 \\, u^{3} + 2320 \\, u^{2} - 512 \\, u + 64}{162 \\, u^{4} - 81 \\, u^{3} - 16 \\, u^{2} + 8 \\, u}}}{4 \\, \\pi {\\left(81 \\, u^{2} - 32 \\, u + 8\\right)} n}\n", "\\end{math}" ], "text/plain": [ "1/4*sqrt(1/6)*(81*u^2 - 8)*(-81/8*(81*u^2 - 8)/(81*u^2 - 48*u + 16))^n*sqrt((6561*u^4 - 5184*u^3 + 2320*u^2 - 512*u + 64)/(162*u^4 - 81*u^3 - 16*u^2 + 8*u))/(pi*(81*u^2 - 32*u + 8)*n)" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "At u = -1/3*I*sqrt(1/3) + 1/3 and u = 1/3*I*sqrt(1/3) + 1/3\n" ] } ], "source": [ "# Get final asymptotic constants\n", "d = 3\n", "A = Q[g]/Q[λ]\n", "B = Q[λ]^(d-1)/Q[h]*Pd^(2-d)\n", "C = Pd/Q[T]\n", "\n", "# Get result of algorithm\n", "var('n')\n", "ASM = sqrt(n*2*pi)^((1-d)) * A * sqrt(B) * SR(C)^n\n", "print(\"Asymptotics are given by adding the evaluations of\")\n", "show(ASM)\n", "print(\"At u = {} and u = {}\".format(u1.radical_expression(),u2.radical_expression()))" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Explicitly, the dominant asymptotic behaviour is\n" ] }, { "data": { "text/html": [ "" ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}{\\mathbf{#1}}\\frac{3 \\, \\left(\\frac{81}{8} i \\, \\sqrt{3}\\right)^{n} \\sqrt{\\frac{1}{2} i \\, \\sqrt{3} + \\frac{1}{2}}}{4 \\, \\pi n} + \\frac{3 \\, \\left(-\\frac{81}{8} i \\, \\sqrt{3}\\right)^{n} \\sqrt{-\\frac{1}{2} i \\, \\sqrt{3} + \\frac{1}{2}}}{4 \\, \\pi n}\n", "\\end{math}" ], "text/plain": [ "3/4*(81/8*I*sqrt(3))^n*sqrt(1/2*I*sqrt(3) + 1/2)/(pi*n) + 3/4*(-81/8*I*sqrt(3))^n*sqrt(-1/2*I*sqrt(3) + 1/2)/(pi*n)" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# We break out the steps to help algebraic simplification\n", "maxima_calculus('algebraic: true;')\n", "A1 = QQbar(SR(A).subs(u==u1)).radical_expression()\n", "B1 = QQbar(SR(B).subs(u==u1)).radical_expression()\n", "C1 = QQbar(SR(C).subs(u==u1)).radical_expression()\n", "\n", "A2 = QQbar(SR(A).subs(u==u2)).radical_expression()\n", "B2 = QQbar(SR(B).subs(u==u2)).radical_expression()\n", "C2 = QQbar(SR(C).subs(u==u2)).radical_expression()\n", "\n", "ASMn = sqrt(n*2*pi)^((1-d))* QQbar(A1*sqrt(B1)).radical_expression() * C1^n\n", "ASMn = ASMn + sqrt(n*2*pi)^((1-d)) * QQbar(A2*sqrt(B2)).radical_expression() * C2^n\n", "\n", "# Print result\n", "print(\"Explicitly, the dominant asymptotic behaviour is\")\n", "show(ASMn)" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "# Using the ore_algebra package we can heuristically check our derived asymptotics\n", "# First, generate the first 20 terms of the diagonal (takes ~20 seconds)\n", "cfs = QQ[x,y,z]((1/H).taylor((x,0),(y,0),(z,0),58))\n", "LST = [cfs[k,k,k] for k in range(20)]" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "At n = 10, the ratio of asympt. to actual sequence is 0.961192243551694\n", "At n = 100, the ratio of asympt. to actual sequence is 0.996111746232913\n", "At n = 1000, the ratio of asympt. to actual sequence is 0.999611117301680\n" ] } ], "source": [ "# Next, use ore_algebra to guess a recurrence\n", "from ore_algebra import *\n", "Ind. = PolynomialRing(QQ); Shift. = OreAlgebra(Ind)\n", "rec = guess(LST,Shift)\n", "\n", "# Use this recurrence to quickly generate 1000 terms of diagonal\n", "LSTlong = rec.to_list(LST[:5],1001)\n", "\n", "# Print results\n", "print(\"At n = 10, the ratio of asympt. to actual sequence is\",(LSTlong/ASMn.subs(n=10)).n())\n", "print(\"At n = 100, the ratio of asympt. to actual sequence is\",(LSTlong/ASMn.subs(n=100)).n())\n", "print(\"At n = 1000, the ratio of asympt. to actual sequence is\",(LSTlong/ASMn.subs(n=1000)).n())" ] }, { "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 }