Ask Your Question
1

Why is the code for the perturbation theory having inaccurate answers? Maple to Sage? (ANSWERED) [closed]

asked 2024-01-30 07:23:48 +0200

SparC gravatar image

updated 2024-02-02 01:46:58 +0200

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)

edit retag flag offensive reopen merge delete

Closed for the following reason the question is answered, right answer was accepted by Max Alekseyev
close date 2024-02-02 02:12:02.858588

Comments

Got the link to the file on google drive at the bottom :) (and here)

SparC gravatar imageSparC ( 2024-02-01 03:31:07 +0200 )edit
1

In the definition of C the exponent of -1 should be k+l instead of k+1. The numerical error probably comes from mixing symbolics (particularly pi) and numerics, but I don't see exactly how (which an answer to the question should explain). A plain numerical approach works:

R = RealField(53)
PI = R(pi)
C = lambda k, l, b: -8*l*k*((b-1)*(-1)^(k+l)-b) / (PI^4*(k-l)^3*(k+l)^3)
S = lambda k, b: -PI^2*(sum((l^2-k^2)*C(k, l ,b)^2 for l in range(1,1001) if l != k))
E = lambda 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(R(1),R(1),R(1/2))
rburing gravatar imagerburing ( 2024-02-01 13:46:41 +0200 )edit

1 Answer

Sort by ยป oldest newest most voted
0

answered 2024-02-02 01:26:58 +0200

SparC gravatar image

updated 2024-02-02 01:46:26 +0200

Ok this question has been answered by rburing on his comment to my question, thank you very much :)

I replaced k+1 with k+l and not much happened, same with the PI, however once I switched from functions to lambdas everything magically fixed and went into the correct accuracy range that I was hoping for. I also played around with the PI and that simply seems to speed up the process and doesn't affect accuracy at all (odd).

Thank you very much for your answer, seems that the problem lay with me being unable to read the source material and use python code correctly lol

EDIT: It also seems to be more accurate than Maple in some cases which is a great result thanks

edit flag offensive delete link more

Question Tools

1 follower

Stats

Asked: 2024-01-30 07:23:48 +0200

Seen: 550 times

Last updated: Feb 02