Implementation of the recursive calculation of the coefficients of the series I, as explained in Remark 4.3¶

We start by defining the space of multivariate polynomials that we will work on (here in SageMath, this is faster than using symbolic variables):

In [1]:
P.<s,nuw,nub,t> = QQ[]

We then define the define the operator making up the right-hand side of equation (19):

In [2]:
def Lambda(J):
    return -2*t^2*nuw*(nub*nuw - 1)*diff(J, nub)- 2*t^2*(nub*nuw - 1)*diff(J, nuw) + 2*t^3*(nuw^2 + nub)*diff(J, t)
In [3]:
Qpap = (32*nuw^8 + 10*nub*nuw^6 + 20*nub^2*nuw^4 + 32*nub^3*nuw^2 + 118*nuw^5 + 32*nub^4 + 108*nub*nuw^3 + 74*nub^2*nuw + 64*nuw^2 + 22*nub)*s*t^3 + (nuw^5 + 2*nuw^2 + nub)*s + 2*(2*nuw^4 + nub*nuw^2 + 2*nub^2 + 3*nuw)^2*t^3 + 2*nuw*(2*nuw^4 + nub*nuw^2 + 2*nub^2 + 3*nuw);
In [4]:
def RHS(J):
    return s/12*Lambda(Lambda(Lambda(Lambda(J)))) +1/2*(Lambda(Lambda(J)))^2+ t*(nuw+2*t^3*(2*nuw^4+nub*nuw^2+2*nub^2+3*nuw))*Lambda(Lambda(J))+t^5*Qpap

We now write the recursion on the coefficients of I given by the non-homogeneous version of equation (21):

In [5]:
@cached_function
def RHSn(n):
    return RHS(IIser(n-3)).coefficient({t: n+2})

@cached_function
def coI(n,i,j):
    if n==0:
        return P.zero()
    elif n%3!=0:
        return 0
    elif ((i-j)%3)!=0:
        return P.zero()
    elif i<0 or j<0 or i+j>n:
        return P.zero()
    elif j==n:
        return 3/n/(3+n) * RHSn(n).coefficient({nuw: n+2})
    elif i==n:
        return coI(n,0,n)
    elif i<j:
        return coI(n,j,i)
    else: 
        return 3/2/(j-n)/(i-j-n)*(RHSn(n).coefficient({nub: i+1, nuw: j}))+ (i + 1 - n)*(i + 2*j - n)*coI(n,i + 1, j - 2)/(2*(j - n)*(i - j - n)) - (i + 2)*(2*i + 3 + j - 2*n)*coI(n,i + 2, j - 1)/((j - n)*(i - j - n)) + (j + 1)*coI(n,i + 1, j + 1)/(j - n) + 3*(i + 3)*(i + 2)*coI(n,i + 3, j)/(2*(j - n)*(i - j - n))

@cached_function
def coIn(n):
    if n==0:
        return P.zero()
    elif n%3!=0:
        return P.zero()
    else:
        r=P.zero()
        for j in range(n+1):
             for k in range(n+1):
                 if ((k-j)%3)==0 & k+j<=n:
                     r+=coI(n,k,j)*nuw^j*nub^k
        return r
        
@cached_function
def IIser(N):
    r=P.zero()
    for n in range(N+1):
        r+=coIn(n)*t^(n)
    return r

Checking the first few terms:

In [6]:
IIser(5) - (nub^3*(2/3 + s/6) + 1/3 + nuw^3*(2/3 + s/6) + nub*nuw + s/3)*t^3
Out[6]:
0

Clearing the cache to evaluate the performance of this implementation:

In [7]:
RHSn.clear_cache()
coI.clear_cache()
coIn.clear_cache()
IIser.clear_cache()
In [8]:
%time I=IIser(30)

CPU times: user 892 ms, sys: 15.7 ms, total: 908 ms
Wall time: 913 ms
In [9]:
RHSn.clear_cache()
coI.clear_cache()
coIn.clear_cache()
IIser.clear_cache()
In [10]:
%time I=IIser(45)

CPU times: user 22 s, sys: 71.1 ms, total: 22.1 s
Wall time: 22.1 s
In [11]:
RHSn.clear_cache()
coI.clear_cache()
coIn.clear_cache()
IIser.clear_cache()
In [12]:
%time I=IIser(54)

CPU times: user 2min 18s, sys: 474 ms, total: 2min 18s
Wall time: 2min 18s
In [13]:
RHSn.clear_cache()
coI.clear_cache()
coIn.clear_cache()
IIser.clear_cache()

This implementation is currently not as efficient as in Maple: it gets too big to terminate for n=57 and above on the author's laptop.

In [ ]:
%time I=IIser(57) #kernel stops before terminating
In [ ]: