Loading [MathJax]/jax/output/HTML-CSS/jax.js
Ask Your Question
2

How to calculate the double sine function via Sage to high precision?

asked 0 years ago

Archaon gravatar image

From a paper of Ruijsenaars, we can extract the following integral representation for the logarithm of the double sine function:

log(S2(w,z))=0(ey(1+wz)eyz(1ey)(1ewy)1y(2zw1w1))dyy

In Sage, it is straightforward to implement this function using the procedure "numerical_integral". It works just fine and produces results in perfect alignment with the known values of the double sine function. (E.g., we have S2(w,1)=w).

My problem is, that I don't know how to increase the number of digits I can get using "numerical_integral". I only get around 16 digits in the decimal expansion, but for my applications I need at least 32.

Question 1: How to increase the number of digits in Sage for the output of the procedure "numerical_integral"? In other words, how can I increase the precision of the procedure "numerical_integral"?

Here is my code:

numerical_integral(( (e^(-(1+w-z)*y)-e^(-z*y))/((1-e^(-y))*(1-e^(-w*y))) -1/y*(2*z/w-1/w-1) )/y,0,+Infinity,max_points=200)

Thank you very much in advance for your help!

PS. Other suggestions to obtain high precision results are welcome as well of course...

PS2. I would be fantastic to have a proper implementation of the double sine function in Sage!

Preview: (hide)

1 Answer

Sort by » oldest newest most voted
2

answered 0 years ago

achrzesz gravatar image

To use mpmath, you have to fix the parameters z, w :

from mpmath import *
mp.dps=100; mp.prett=True
z=1.0
w=2.0
f = lambda y:((e^(-(1+w-z)*y)-e^(-z*y))/((1-e^(-y))*(1-e^(-w*y))) -1/y*(2*z/w-1/w-1) )/y
ii=mp.quad(f,[0.0,mp.inf])
print(ii)
0.6868438210465139334339462432252338314613985881041879632836891912354724713670427615704250909648475892

To investigate more parameters you can use loops:

from mpmath import *
mp.dps=100; mp.prett=True
for w in [1,1.5,2]:
    for z in [2.,1.5,1.]:
        f = lambda y:((e^(-(1+w-z)*y)-e^(-z*y))/((1-e^(-y))*(1-e^(-w*y))) -1/y*(2*z/w-1/w-1) )/y
        ii=mp.quad(f,[0.0,mp.inf])
        print("%0.2f"%w,'\t',"%0.2f"%z,'\t',N(ii))  

1.00     2.00    238.130061186454
1.00     1.50    -0.403748148636294
1.00     1.00    0.000000000000000
1.50     2.00    -1.32360608384651
1.50     1.50    -5.27628283903327e88
1.50     1.00    1.31907070975832e88
2.00     2.00    -0.686843821046514
2.00     1.50    0.000000000000000
2.00     1.00    0.686843821046514
Preview: (hide)
link

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower

Stats

Asked: 0 years ago

Seen: 101 times

Last updated: Dec 21 '24