Why is the code for the perturbation theory having inaccurate answers? Maple to Sage? (ANSWERED) [closed]
Hi there,
my professor uses Maple and he had posted a Maple file for use to use to compute some values, however he gave it to us as part of an online summer research project and I am off-campus so I don't have Maple, so I'm using Sage. The method used to get the values is perturbation theory used in the context of the confined quantum harmonic oscillator.
The problem is that the value for the function that my professor got for E(1, 1.0, 0.5) from the Maple file is: 9.902258637, whereas the one that is generated from this code gets close, but loses accuracy after the some decimals: 9.902277143
I have pasted my code below:
k = var('k')
l = var('l')
b = var('b')
alpha = var('alpha')
assume(k>0)
assume(l>0)
assume(k, 'integer')
assume(l, 'integer')
C(k, l, b) = -8*l*k*((b-1)*(-1)^(k+1)-b) / (pi^4*(k-l)^3*(k+l)^3)
expression = lambda l : (l^2-k^2)*C(k, l ,b)^2
S(k, b) = -pi^2*(sum([expression(l) for l in range(1,1001) if l != k]))
E(k, alpha, b) = k^2*pi^2 + alpha*((b-1/2)^2 + 1/12 - 1/(2*k^2*pi^2)) + alpha^2*S(k,b)
E(1,1,0.5).n(digits=10)
Then I noticed that the sum I used in S(k,b) actually outputs zero, so only the first order correction is used. And even when I attempt to make it not equal zero by removing the condition and splitting it into 2 sums, it actually becomes less accurate, E(k_value, 1, 0.5) = 9.901161373, for k_value = 1.
S(k, b) = -pi^2*(sum([expression(l) for l in range(1,k_value)]) + sum([expression(l) for l in range(k_value+1,1001)]))
Does anybody have any idea why this happens and why when the second order correction is 0 the answer is actually more accurate then when it isn't?
Note: I tried to upload an image of the Maple file but I cannot yet since I am new but if someone is interested I can attempt to give them a screenshot somehow or the Maple file.
EDIT: Here is the google drive link to both the maple file and the pdf version : Google Drive Link to Folder Equation 11 looks like the problem child. (Equations above 8 are just background calculations)
Got the link to the file on google drive at the bottom :) (and here)
In the definition of
C
the exponent of-1
should bek+l
instead ofk+1
. The numerical error probably comes from mixing symbolics (particularlypi
) and numerics, but I don't see exactly how (which an answer to the question should explain). A plain numerical approach works: