Ask Your Question
1

numeric precision unexpectedly low

asked 2021-07-14 12:38:50 +0200

anonymous user

Anonymous

updated 2021-07-15 18:04:38 +0200

FrédéricC gravatar image

I try to integrate a probability density function over an fixed interval.

Can't give the like but it is a "relativistic Breit-Wigner Distribution" as in the English Wikipedia.

Gamma = 294000
m = 3686100
gamma = sqrt(m^2*(m^2+Gamma^2))
k = (2*sqrt(2)*m*Gamma*gamma)/(pi*sqrt(m^2+gamma))
BW = k/((x^2-m^2)^2+m^2*Gamma^2)
test = integrate(BW(x),x,m-20*Gamma,m+20*Gamma)

Till here I shouldn't have lost any precision. but when I convert it numerically with

n(test,2000)

I get:

1.003568780514681625312584816997550190541415598167450102773034501594600296432162707154049825083671452461581265874619355037287669117698123644188958803147725041330150928274806048238042399517154207573726037521427349383968841867335044816048064092169439887022730672979989338003256477687641391106254297483937686718675512005461848316884186358220468903611061147782405268579702217582451508116179810896698136112722502668649975423152214020083815884479176920571802765348290248949692543327981513253681778869531348343181499633971108206392658073359294594919457298263920263466439149236869382936665388139526661407807092

Since this is an integral over a probability density it must be less or equal to unity. To make things worse the test I used to catch this failed. When I print n(test,2000)-n(test,5000) I get:

1.548076368777034133563243683402305144019545879147318729129474813190957966982407215191933066512274489911906109674015771900520547143254593904332062734353178144359258512027924010153347526193868596270002650844292344493806124142870891743921764286481744814534440665865632030746674736110661609478466022255757297949839318763567151082721152805813677624313092475594710476219357188434129215184208643483256856632402761545583382935944646898012751255575740279260126952010697251536670906745564901191471261083975404910237568882848368855134234892724370503380053787304268124196831149676013601785994763640281219314046884e-517

So it seems as if 500 digits are correct. Is this a bug or am I not using this correctly? I'm using sage 9.0 as packaged in Ubuntu.

edit retag flag offensive close merge delete

1 Answer

Sort by » oldest newest most voted
3

answered 2021-07-14 21:37:46 +0200

tmonteil gravatar image

You should have a look at the expression test, wich is a very huge. Hence, when turned numerical, a lot of roundings accumulate. Instead, you could compute the numerical integral directly (the second part is the precision of the result):

sage: numerical_integral(BW,m-20*Gamma,m+20*Gamma)
(1.039353768756414, 7.667179586654595e-10)

As you can see, the integral is larger than 1. If you want more precision, you can use real ball field:

sage: CC = ComplexBallField(1000)
sage: f = fast_callable(BW,vars=[x],domain=CC)
sage: CC.integral(lambda x,_: f(x) ,CC(m-20*Gamma),CC(m+20*Gamma))
[1.0393537687564144005635175005493502562514240251865059997675473690905082074521147152980391279501557766285787071456276287443389069265579835090610752201107443966192623303858945758537114342432084054091316038569497030744623566434274060346642219779086790951559180203731682971294215371002868275171917298204 +/- 4.40e-299]
edit flag offensive delete link more

Comments

Thank You! If the integral is larger than one it cannot be a probability density function. In principle that is what the normalisation constant k should take care of. I need to do more checks on the formula. Thanks again!

just_tryin_to_use_this gravatar imagejust_tryin_to_use_this ( 2021-07-15 08:53:09 +0200 )edit

Great. This answer has just helped me a lot with a complex integral I could not solve before.

mwageringel gravatar imagemwageringel ( 2021-07-15 20:46:51 +0200 )edit

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: 2021-07-14 12:35:41 +0200

Seen: 322 times

Last updated: Jul 15 '21