Ask Your Question

Revision history [back]

click to hide/show revision 1
initial version

answered 2023-01-01 20:51:50 +0100

achrzesz gravatar image

I'm rather sceptic about symbolic investigations of big, dense matrices, but it seems that for n=9 your determinant is computable

import time

P.<E>=QQ[]
energy=E

def expected_distance(n, l):
    if n <= -2:
        return 0
    if n == -1:
        return -2 * energy
    if n == 0:
        return 1
    return  (
        - 4 * (2 * (n + 1) - 1) * expected_distance(n - 1, l) \
        - n * ((n + 1) * (n - 1) - 4 * l * (l + 1)) * expected_distance(n - 2, l)
        ) / (8 * (n + 1) * energy)

l=2
rank=9
hankel = matrix.hankel([expected_distance(i + offset, l) for i in range(rank)],
        [expected_distance(i + offset + rank - 1, l) for i in range(1, rank)], P.fraction_field())

%time hd=hankel.det()

CPU times: user 3.39 s, sys: 3.91 ms, total: 3.39 s Wall time: 3.38 s

Finding positive values:

pos_sol=solve(SR(hd)>0,var('E'))

I'm rather sceptic sceptical about symbolic investigations of big, dense matrices, but it seems that for n=9 your determinant is computable

import time

P.<E>=QQ[]
energy=E

def expected_distance(n, l):
    if n <= -2:
        return 0
    if n == -1:
        return -2 * energy
    if n == 0:
        return 1
    return  (
        - 4 * (2 * (n + 1) - 1) * expected_distance(n - 1, l) \
        - n * ((n + 1) * (n - 1) - 4 * l * (l + 1)) * expected_distance(n - 2, l)
        ) / (8 * (n + 1) * energy)

l=2
rank=9
hankel = matrix.hankel([expected_distance(i + offset, l) for i in range(rank)],
        [expected_distance(i + offset + rank - 1, l) for i in range(1, rank)], P.fraction_field())

%time hd=hankel.det()

CPU times: user 3.39 s, sys: 3.91 ms, total: 3.39 s Wall time: 3.38 s

Finding positive values:

pos_sol=solve(SR(hd)>0,var('E'))

I'm rather sceptical about symbolic investigations of big, dense matrices, but it seems that for n=9 n=12 your determinant is computable

import time

P.<E>=QQ[]
energy=E

@cached_function
def expected_distance(n, l):
    if n <= -2:
        return 0
    if n == -1:
        return -2 * energy
    if n == 0:
        return 1
    return  (
        - 4 * (2 * (n + 1) - 1) * expected_distance(n - 1, l) \
        - n * ((n + 1) * (n - 1) - 4 * l * (l + 1)) * expected_distance(n - 2, l)
        ) / (8 * (n + 1) * energy)

offset=1
l=2
rank=9
hankel rank=12

hank = matrix.hankel([expected_distance(i + offset, l) for i in range(rank)],
        [expected_distance(i + offset + rank - 1, l) for i in range(1, rank)], rank)],
                     P.fraction_field())

%time hd=hankel.det()
hd=hank.det()

CPU times: user 3.39 s, 1min 48s, sys: 3.91 212 ms, total: 3.39 s 1min 48s Wall time: 3.38 s1min 47s

The value of hd is of the form (a*E^73+bE^72+...)/E^144, so numerical computations in usual precision are not adequate

Finding positive values:values works:

pos_sol=solve(SR(hd)>0,var('E'))

but its precision may be also to low. Using for example

sgn=[(r,sign(RR(hd(RealField(3000)(r))))) for r in srange(-0.5,0.5,0.001)]

one can check it numerically.

If you need faster method use pure numerical approach

import time
# I don't now if 3000 is a good choice,
# cf hd value below
R3000=RealField(3000)    
energy=R3000(-0.1)

@cached_function
def expected_distance(n, l):
    if n <= -2:
        return 0
    if n == -1:
        return -2 * energy
    if n == 0:
        return 1
    return  (
        - 4 * (2 * (n + 1) - 1) * expected_distance(n - 1, l) \
        - n * ((n + 1) * (n - 1) - 4 * l * (l + 1)) * expected_distance(n - 2, l)
        ) / (8 * (n + 1) * energy)

l=2
offset=1
rank=20
a1=[expected_distance(i + offset, l) for i in range(rank)];
a2=[expected_distance(i + offset + rank - 1, l) for i in range(1, rank)]
v1=vector(R3000,a1)
v2=vector(R3000,a2)
hank = matrix.hankel(v1,v2,R3000)
%time hd=hank.det()

CPU times: user 66.9 ms, sys: 0 ns, total: 66.9 ms Wall time: 66.4 ms

hd.n(20)

4.3908e394

I'm rather sceptical about symbolic investigations of big, dense matrices, but it seems that for n=12 your determinant is computable

import time

P.<E>=QQ[]
energy=E

@cached_function
def expected_distance(n, l):
    if n <= -2:
        return 0
    if n == -1:
        return -2 * energy
    if n == 0:
        return 1
    return  (
        - 4 * (2 * (n + 1) - 1) * expected_distance(n - 1, l) \
        - n * ((n + 1) * (n - 1) - 4 * l * (l + 1)) * expected_distance(n - 2, l)
        ) / (8 * (n + 1) * energy)

offset=1
l=2
rank=12
rank=20

hank = matrix.hankel([expected_distance(i + offset, l) for i in range(rank)],
        [expected_distance(i + offset + rank - 1, l) for i in range(1, rank)],
                     P.fraction_field())

# P and L have determinans 1
%time hd=hank.det()
P,L,U=hankel.LU(); d1=mul(U.diagonal())

CPU times: user 1min 48s, 1.19 s, sys: 212 ms, 0 ns, total: 1min 48s 1.19 s Wall time: 1min 47s1.19 s

The value of hd is of the form (a*E^73+bE^72+...)/E^144, (a*E^201+bE^200+...)/E^400, so numerical computations in usual precision are not adequate

Finding positive values works:

pos_sol=solve(SR(hd)>0,var('E'))

but its precision may be also to low. Using for example

sgn=[(r,sign(RR(hd(RealField(3000)(r))))) for r in srange(-0.5,0.5,0.001)]
nu=d1.numerator()
sol=nu.roots(ring=RealField(200))

one can check it numerically.with higher precision.

If you need faster method use pure numerical approach

import time
# I don't now if 3000 is a good choice,
# cf hd value below
R3000=RealField(3000)    
energy=R3000(-0.1)

@cached_function
def expected_distance(n, l):
    if n <= -2:
        return 0
    if n == -1:
        return -2 * energy
    if n == 0:
        return 1
    return  (
        - 4 * (2 * (n + 1) - 1) * expected_distance(n - 1, l) \
        - n * ((n + 1) * (n - 1) - 4 * l * (l + 1)) * expected_distance(n - 2, l)
        ) / (8 * (n + 1) * energy)

l=2
offset=1
rank=20
a1=[expected_distance(i + offset, l) for i in range(rank)];
a2=[expected_distance(i + offset + rank - 1, l) for i in range(1, rank)]
v1=vector(R3000,a1)
v2=vector(R3000,a2)
hank = matrix.hankel(v1,v2,R3000)
%time hd=hank.det()

CPU times: user 66.9 ms, sys: 0 ns, total: 66.9 ms Wall time: 66.4 ms

hd.n(20)

4.3908e394

I'm rather sceptical about symbolic investigations of big, dense matrices, but it seems that for n=12 your determinant is computable

import time

P.<E>=QQ[]
energy=E

@cached_function
def expected_distance(n, l):
    if n <= -2:
        return 0
    if n == -1:
        return -2 * energy
    if n == 0:
        return 1
    return  (
        - 4 * (2 * (n + 1) - 1) * expected_distance(n - 1, l) \
        - n * ((n + 1) * (n - 1) - 4 * l * (l + 1)) * expected_distance(n - 2, l)
        ) / (8 * (n + 1) * energy)

offset=1
l=2
rank=20

hank = matrix.hankel([expected_distance(i + offset, l) for i in range(rank)],
        [expected_distance(i + offset + rank - 1, l) for i in range(1, rank)],
                     P.fraction_field())

# P and L have determinans 1
%time P,L,U=hankel.LU(); d1=mul(U.diagonal())

CPU times: user 1.19 s, sys: 0 ns, total: 1.19 s Wall time: 1.19 s

The value of hd d1 is of the form (a*E^201+bE^200+...)/E^400, so numerical computations in usual precision are not adequate

Finding positive values works:

pos_sol=solve(SR(hd)>0,var('E'))
pos_sol=solve(SR(d1)>0,var('E'))

but its precision may be also to low. Using for example

nu=d1.numerator()
sol=nu.roots(ring=RealField(200))

one can check it with higher precision.

If you need faster method use pure numerical approach

import time
# I don't now if 3000 is a good choice,
# cf hd value below
R3000=RealField(3000)    
energy=R3000(-0.1)

@cached_function
def expected_distance(n, l):
    if n <= -2:
        return 0
    if n == -1:
        return -2 * energy
    if n == 0:
        return 1
    return  (
        - 4 * (2 * (n + 1) - 1) * expected_distance(n - 1, l) \
        - n * ((n + 1) * (n - 1) - 4 * l * (l + 1)) * expected_distance(n - 2, l)
        ) / (8 * (n + 1) * energy)

l=2
offset=1
rank=20
a1=[expected_distance(i + offset, l) for i in range(rank)];
a2=[expected_distance(i + offset + rank - 1, l) for i in range(1, rank)]
v1=vector(R3000,a1)
v2=vector(R3000,a2)
hank = matrix.hankel(v1,v2,R3000)
%time hd=hank.det()

CPU times: user 66.9 ms, sys: 0 ns, total: 66.9 ms Wall time: 66.4 ms

hd.n(20)

4.3908e394

I'm rather sceptical about symbolic investigations of big, dense matrices, but it seems that for n=12 n=20 your determinant is computable

import time

P.<E>=QQ[]
energy=E

@cached_function
def expected_distance(n, l):
    if n <= -2:
        return 0
    if n == -1:
        return -2 * energy
    if n == 0:
        return 1
    return  (
        - 4 * (2 * (n + 1) - 1) * expected_distance(n - 1, l) \
        - n * ((n + 1) * (n - 1) - 4 * l * (l + 1)) * expected_distance(n - 2, l)
        ) / (8 * (n + 1) * energy)

offset=1
l=2
rank=20

hank = matrix.hankel([expected_distance(i + offset, l) for i in range(rank)],
        [expected_distance(i + offset + rank - 1, l) for i in range(1, rank)],
                     P.fraction_field())

# P and L have determinans 1
%time P,L,U=hankel.LU(); d1=mul(U.diagonal())

CPU times: user 1.19 s, sys: 0 ns, total: 1.19 s Wall time: 1.19 s

The value of d1 is of the form (a*E^201+bE^200+...)/E^400, so numerical computations in usual precision are not adequate

Finding positive values works:

pos_sol=solve(SR(d1)>0,var('E'))

but its precision may be also to low. Using for example

nu=d1.numerator()
sol=nu.roots(ring=RealField(200))

one can check it with higher precision.

If you need faster method use pure numerical approach

import time
# I don't now if 3000 is a good choice,
# cf hd value below
R3000=RealField(3000)    
energy=R3000(-0.1)

@cached_function
def expected_distance(n, l):
    if n <= -2:
        return 0
    if n == -1:
        return -2 * energy
    if n == 0:
        return 1
    return  (
        - 4 * (2 * (n + 1) - 1) * expected_distance(n - 1, l) \
        - n * ((n + 1) * (n - 1) - 4 * l * (l + 1)) * expected_distance(n - 2, l)
        ) / (8 * (n + 1) * energy)

l=2
offset=1
rank=20
a1=[expected_distance(i + offset, l) for i in range(rank)];
a2=[expected_distance(i + offset + rank - 1, l) for i in range(1, rank)]
v1=vector(R3000,a1)
v2=vector(R3000,a2)
hank = matrix.hankel(v1,v2,R3000)
%time hd=hank.det()

CPU times: user 66.9 ms, sys: 0 ns, total: 66.9 ms Wall time: 66.4 ms

hd.n(20)

4.3908e394

I'm rather sceptical about symbolic investigations of big, dense matrices, but it seems that for n=20 your determinant is computable

import time

P.<E>=QQ[]
energy=E

@cached_function
def expected_distance(n, l):
    if n <= -2:
        return 0
    if n == -1:
        return -2 * energy
    if n == 0:
        return 1
    return  (
        - 4 * (2 * (n + 1) - 1) * expected_distance(n - 1, l) \
        - n * ((n + 1) * (n - 1) - 4 * l * (l + 1)) * expected_distance(n - 2, l)
        ) / (8 * (n + 1) * energy)

offset=1
l=2
rank=20

hank = matrix.hankel([expected_distance(i + offset, l) for i in range(rank)],
        [expected_distance(i + offset + rank - 1, l) for i in range(1, rank)],
                     P.fraction_field())

# P and L have determinans 1
 %time P,L,U=hankel.LU(); d1=mul(U.diagonal())
L, d = hank.indefinite_factorization('symmetric');d1=mul(d)

CPU times: user 1.19 s, 776 ms, sys: 0 ns, total: 1.19 s 776 ms Wall time: 1.19 s775 ms

The value of d1 is of the form (a*E^201+bE^200+...)/E^400, so numerical computations in usual precision are not adequate

Finding positive values works:

pos_sol=solve(SR(d1)>0,var('E'))

but its precision may be also to low. Using for example

nu=d1.numerator()
sol=nu.roots(ring=RealField(200))

one can check it with higher precision.

If you need faster method use pure numerical approach

import time
# I don't now if 3000 is a good choice,
# cf hd value below
R3000=RealField(3000)    
energy=R3000(-0.1)

@cached_function
def expected_distance(n, l):
    if n <= -2:
        return 0
    if n == -1:
        return -2 * energy
    if n == 0:
        return 1
    return  (
        - 4 * (2 * (n + 1) - 1) * expected_distance(n - 1, l) \
        - n * ((n + 1) * (n - 1) - 4 * l * (l + 1)) * expected_distance(n - 2, l)
        ) / (8 * (n + 1) * energy)

l=2
offset=1
rank=20
a1=[expected_distance(i + offset, l) for i in range(rank)];
a2=[expected_distance(i + offset + rank - 1, l) for i in range(1, rank)]
v1=vector(R3000,a1)
v2=vector(R3000,a2)
hank = matrix.hankel(v1,v2,R3000)
%time hd=hank.det()

CPU times: user 66.9 ms, sys: 0 ns, total: 66.9 ms Wall time: 66.4 ms

hd.n(20)

4.3908e394

Copyright Sage, 2010. Some rights reserved under creative commons license. Content on this site is licensed under a Creative Commons Attribution Share Alike 3.0 license.
about | faq | help | privacy policy | terms of service
Powered by Askbot version 0.7.59