wxmaxima, why my function don't integrate well?

Question

jk_hermite(n,v):=expand((-1)^n * exp(v^2)*diff(exp(-v^2),v,n));

jk_prob1(k):=block( [],
 int1:integrate(abs(jk_hermite(k,x)*exp(-x^2/2))^2,x,minf,inf),
 b:sqrt(1/int1),
 int2:integrate(abs(jk_hermite(k,x)*b*exp(-x^2/2))^2,x,-sqrt(2*k+1),sqrt(2*k+1)),
 p:1-int2
);

This is a program (jk_prob1) that calculate a number (p) for given number (k). Working fine for k<27 according to bibliography, but for higher k there is some problem probably with int2 integration.

below we have a list [ [k,jk_prob(k)], ...] with the calculations for k in 0 to 35. And you can see that integration start oscillating after k=27.

Is that a problem with the arithmetic of my machine or a problem with integration() function?

(%i174) data1:makelist([k,jk_prob1(k)],k,0,35,1);
(%o174) [[0,0.1572992070502853],[1,0.1116102250947125],[2,0.09506943488572384],[3,0.08548286439326391],[4,0.07892635105633516],
[5,0.07403423706609669],[6,0.0701809281050576],[7,0.06703127879093984],[8,0.06438634028251189],[9,0.06211907732553879],[10,0.06014381448969242],
[11,0.05840025530396975],[12,0.05684448915503315],[13,0.05544362910545952],[14,0.0541724601132354],[15,0.05301126092266206],
[16,0.05194434169894979],[17,0.05095903331123497],[18,0.05004496080669607],[19,0.04919359561921899],[20,0.04839756630363612],
[21,0.04765147770002887],[22,0.04694842144159583],[23,0.04628695370710445],[24,0.04566365904566005],[25,0.04505081797459654],
[26,0.04457118540829141],[27,0.04373065855907632],[28,0.0441370340921643],[29,0.04112841762932007],[30,0.04805199290030093],
[31,0.0318955398934323],[32,0.0685550172107956],[33,-0.02410444428168423],[34,0.1631563092724104],[35,-0.2711053519616526]]

In graph you can see the results jk_prob1(k) according to k

We expect the results to decrease continuously but certainly not to oscillating.

(%i176) build_info();
(%o176) build_info(version=
"5.46.0",timestamp=
"2022-04-13 23:24:03",host=
"x86_64-w64-mingw32",
lisp_name="SBCL",
lisp_version="2.2.2",
,maxima_frontend=
"wxMaxima",
maxima_frontend_version=
"22.04.0_MSW")

Thank you in advance.

Interesting problem, I'll try to look at it later. Some things to look at. (1) Try plotting int1 and int2 separately. (2) Try plotting int1 and int2 for smaller increments of k (i.e. noninteger values). (3) Also plot jk_hermite. From the definition it seems possible it's an oscillating function, so it's at least plausible that its integrals are also oscillating. Something to look at, anyway. — Robert Dodier, Feb 27 '23 at 18:35
Hmm, I'm not seeing the oscillating behavior for k = 0, ..., 35, I see a smooth curve. Some of the terms are very large; maybe there is a problem with numerical evaluation of those. You could try `fpprec: nnn` where nnn = 100 or something, and then say `p: bfloat(1 - p)` at the end to get a bigfloat (variable precision float). — Robert Dodier, Feb 28 '23 at 02:28
Well, the bad behavior start at k=27, as you can see in the graph and of course negative results are unacceptable. For smaller numbers the program working well. Also k must be positive integer from the definition of jk_hermite() as the order of differential. I have try fpprec and bfloat but i have the same behavior. — john, Mar 01 '23 at 14:40
Right, after looking at it again, I see k must be an integer. I tried the function just as you showed it, and evaluating with k up to 35 doesn't show oscillations, just a monotone decrease. I wonder if there is something else in your Maxima session which is affecting the results. Did you, perhaps, enable the `keepfloat` flag or other global variables? Can you start a new session and copy the function definition as shown and then `makelist` and paste the results here? — Robert Dodier, Mar 01 '23 at 16:05
Thank you. I start a new session and the problem solved, as you suggested! After a few tries i realize that the flag ```numer:true``` have that bad effect to calculations. — john, Mar 03 '23 at 13:30
When I try `ev(jk_prob1(35), numer = true)`, I get some messages about converting rational numbers to floats. That in itself is okay (it is the expected function of `numer`), but I see some of the numbers are pretty large, e.g. 1e38 or so. If several large numbers are combining to make a result that's much smaller, there's a possibility of losing precision in the result. My first guess is that's what's happening here. (A secondary guess is that `numer` is causing the integration code to return incorrect results -- there are some known bugs having to do with float numbers in integration.) — Robert Dodier, Mar 03 '23 at 17:39

score 0 · Accepted Answer · answered Mar 01 '23 at 18:20

I can't reproduce the error you are seeing; I'm getting what appears to be a correct result. Maybe if you start over you will get the same result? If you start over and get a different result, maybe you can update your question to show the new session.

(%i2) display2d: false $

(%i3) jk_hermite(n,v):=expand((-1)^n * exp(v^2)*diff(exp(-v^2),v,n)) $

(%i4) jk_prob1(k):=block( [],
 int1:integrate(abs(jk_hermite(k,x)*exp(-x^2/2))^2,x,minf,inf),
 b:sqrt(1/int1),
 int2:integrate(abs(jk_hermite(k,x)*b*exp(-x^2/2))^2,x,-sqrt(2*k+1),sqrt(2*k+1)),
 p:1-int2
) $

(%i5) data1:makelist([k,jk_prob1(k)],k,0,35,1);

(%o5) [[0,1-erf(1)],
       [1,
        1-(%e^-3*(%e^3*sqrt(%pi)*erf(sqrt(3))-2*sqrt(3)))
          /sqrt(%pi)],
       [2,
        1-(%e^-5*(4*%e^5*sqrt(%pi)*erf(sqrt(5))-44*sqrt(5)))
          /(4*sqrt(%pi))],
       [3,
        1-(%e^-7*(24*%e^7*sqrt(%pi)*erf(sqrt(7))-1504*sqrt(7)))
          /(24*sqrt(%pi))],
       [4,
        1-(%e^-9*(192*%e^9*sqrt(%pi)*erf(3)-217584))
          /(192*sqrt(%pi))],
       [5,
        1-(%e^-11*(1920*%e^11*sqrt(%pi)*erf(sqrt(11))
                  -4548160*sqrt(11)))
          /(1920*sqrt(%pi))],
       [6,
        1-(%e^-13*(23040*%e^13*sqrt(%pi)*erf(sqrt(13))
                  -351666496*sqrt(13)))
          /(23040*sqrt(%pi))],
       [7,
        1-(%e^-15*(322560*%e^15*sqrt(%pi)*erf(sqrt(15))
                  -2156467968*15^(3/2)))
          /(322560*sqrt(%pi))],
       [8,
        1-(%e^-17*(5160960*%e^17*sqrt(%pi)*erf(sqrt(17))
                  -3450490725632*sqrt(17)))
          /(5160960*sqrt(%pi))],
       [9,
        1-(%e^-19*(92897280*%e^19*sqrt(%pi)*erf(sqrt(19))
                  -418814087689216*sqrt(19)))
          /(92897280*sqrt(%pi))],
       [10,
        1-(%e^-21*(1857945600*%e^21*sqrt(%pi)*erf(sqrt(21))
                  -2714276435051520*21^(3/2)))
          /(1857945600*sqrt(%pi))],
       [11,
        1-(%e^-23*(40874803200*%e^23*sqrt(%pi)*erf(sqrt(23))
                  -8597150358710259712*sqrt(23)))
          /(40874803200*sqrt(%pi))],
       [12,
        1-(%e^-25*(980995276800*%e^25*sqrt(%pi)*erf(5)
                  -7116923021525797376000))
          /(980995276800*sqrt(%pi))],
       [13,
        1-(%e^-27*(25505877196800*%e^27*sqrt(%pi)*erf(3^(3/2))
                  -1056159977061224660992*3^(13/2)))
          /(25505877196800*sqrt(%pi))],
       [14,
        1-(%e^-29*(714164561510400*%e^29*sqrt(%pi)*erf(sqrt(29))
                  -50060221449301592618909696*sqrt(29)))
          /(714164561510400*sqrt(%pi))],
       [15,
        1-(%e^-31*(21424936845312000*%e^31*sqrt(%pi)
                                    *erf(sqrt(31))
                  -10502935872905789447643136000*sqrt(31)))
          /(21424936845312000*sqrt(%pi))],
       [16,
        1-(%e^-33*(685597979049984000*%e^33*sqrt(%pi)
                                     *erf(sqrt(33))
                  -71470974634380182427966701568*33^(3/2)))
          /(685597979049984000*sqrt(%pi))],
       [17,
        1-(%e^-35*(23310331287699456000*%e^35*sqrt(%pi)
                                       *erf(sqrt(35))
                  -460766937322557657585603051520*35^(5/2)))
          /(23310331287699456000*sqrt(%pi))],
       [18,
        1-(%e^-37*(839171926357180416000*%e^37*sqrt(%pi)
                                        *erf(sqrt(37))
                  -143410601721679053521909949555539968
                   *sqrt(37)))
          /(839171926357180416000*sqrt(%pi))],
       [19,
        1-(%e^-39*(31888533201572855808000
                  *%e^39*sqrt(%pi)*erf(sqrt(39))
                  -988566037943992213347770624124125184
                   *39^(3/2)))
          /(31888533201572855808000*sqrt(%pi))],
       [20,
        1-(%e^-41*(1275541328062914232320000
                  *%e^41*sqrt(%pi)*erf(sqrt(41))
                  -10933914922854583659978082324446398382080
                   *sqrt(41)))
          /(1275541328062914232320000*sqrt(%pi))],
       [21,
        1-(%e^-43*(53572735778642397757440000
                  *%e^43*sqrt(%pi)*erf(sqrt(43))
                  -3262278905479206540619172723651100811460608
                   *sqrt(43)))
          /(53572735778642397757440000*sqrt(%pi))],
       [22,
        1-(%e^-45*(2357200374260265501327360000
                  *%e^45*sqrt(%pi)*erf(3*sqrt(5))
                  -4903257423120959495745281094360860593225728
                   *5^(9/2)))
          /(2357200374260265501327360000*sqrt(%pi))],
       [23,
        1-(%e^-47*(108431217215972213061058560000
                  *%e^47*sqrt(%pi)*erf(sqrt(47))
                  -334948078254017640663171030390911909706663460864
                   *sqrt(47)))
          /(108431217215972213061058560000*sqrt(%pi))],
       [24,
        1-(%e^-49*(5204698426366666226930810880000
                  *%e^49*sqrt(%pi)*erf(7)
                  -803416542526033681542640776803851745555237602066432))
          /(5204698426366666226930810880000*sqrt(%pi))],
       [25,
        1-(%e^-51*(260234921318333311346540544000000
                  *%e^51*sqrt(%pi)*erf(sqrt(51))
                  -804382628000720402412181989619600056578667593072640
                   *51^(3/2)))
          /(260234921318333311346540544000000*sqrt(%pi))],
       [26,
        1-(%e^-53*(13532215908553332190020108288000000
                  *%e^53*sqrt(%pi)*erf(sqrt(53))
                  -15268863879626398704434661565976027260923048925715234816
                   *sqrt(53)))
          /(13532215908553332190020108288000000*sqrt(%pi))],
       [27,
        1-(%e^-55*(730739659061879938261085847552000000
                  *%e^55*sqrt(%pi)*erf(sqrt(55))
                  -1953238901674385528202656418280853915490371302850560000
                   *55^(5/2)))
          /(730739659061879938261085847552000000*sqrt(%pi))],
       [28,
        1-(%e^-57*(40921420907465276542620807462912000000
                  *%e^57*sqrt(%pi)*erf(sqrt(57))
                  -41643536862895126218478556664899501086758124485946676084736
                   *57^(3/2)))
          /(40921420907465276542620807462912000000*sqrt(%pi))],
       [29,
        1-(%e^-59*(2373442412632986039472006832848896000000
                  *%e^59*sqrt(%pi)*erf(sqrt(59))
                  -988655575586727665988518643617635558359966422816622541847658496
                   *sqrt(59)))
          /(2373442412632986039472006832848896000000
           *sqrt(%pi))],
       [30,
        1-(%e^-61*(142406544757979162368320409970933760000000
                  *%e^61*sqrt(%pi)*erf(sqrt(61))
                  -426385466109686794974521830955494084352091033957137964359263191040
                   *sqrt(61)))
          /(142406544757979162368320409970933760000000
           *sqrt(%pi))],
       [31,
        1-(%e^-63*(8829205774994708066835865418197893120000000
                  *%e^63*sqrt(%pi)*erf(3*sqrt(7))
                  -237637173843068710884180931448387629810402892660900975740616441856
                   *7^(9/2)))
          /(8829205774994708066835865418197893120000000
           *sqrt(%pi))],
       [32,
        1-(%e^-65*(565069169599661316277495386764665159680000000
                  *%e^65*sqrt(%pi)*erf(sqrt(65))
                  -20743919730124955449185794118204935548012355489446073210608025600000
                   *65^(5/2)))
          /(565069169599661316277495386764665159680000000
           *sqrt(%pi))],
       [33,
        1-(%e^-67*(37294565193577646874314695526467900538880000000
                  *%e^67*sqrt(%pi)*erf(sqrt(67))
                  -41682427353927433909745163803213551387672794288391014517753878087599128576
                   *sqrt(67)))
          /(37294565193577646874314695526467900538880000000
           *sqrt(%pi))],
       [34,
        1-(%e^-69*(2536030433163279987453399295799817236643840000000
                  *%e^69*sqrt(%pi)*erf(sqrt(69))
                  -296226359950091579306352646020536576938264515830526239985823226493501702144
                   *69^(3/2)))
          /(2536030433163279987453399295799817236643840000000
           *sqrt(%pi))],
       [35,
        1-(%e^-71*(177522130321429599121737950705987206565068800000000
                  *%e^71*sqrt(%pi)*erf(sqrt(71))
                  -10324828777648860534610558615880230477741714899112291194458537009337241100615680
                   *sqrt(71)))
          /(177522130321429599121737950705987206565068800000000
           *sqrt(%pi))]]
(%i6) float (%);

(%o6) [[0.0,0.1572992070502852],[1.0,0.1116102250947126],
       [2.0,0.09506943488572639],[3.0,0.08548286439326436],
       [4.0,0.07892635105632473],[5.0,0.0740342370661331],
       [6.0,0.07018092810506227],[7.0,0.06703127879070192],
       [8.0,0.06438634028180634],[9.0,0.06211907732440791],
       [10.0,0.06014381449057227],[11.0,0.05840025529739457],
       [12.0,0.05684448917464546],[13.0,0.05544362909885237],
       [14.0,0.05417246007221554],[15.0,0.05301126098950437],
       [16.0,0.05194434166388195],[17.0,0.05095903209868335],
       [18.0,0.05004496694560001],[19.0,0.04919356800623786],
       [20.0,0.04839766284623237],[21.0,0.04765119897461811],
       [22.0,0.04694902640744292],[23.0,0.04628673000680339],
       [24.0,0.04566049861172194],[25.0,0.04506702174578892],
       [26.0,0.04450340725876845],[27.0,0.04396711504526329],
       [28.0,0.04345590324287807],[29.0,0.0429677842131182],
       [30.0,0.0425009882611056],[31.0,0.0420539335290957],
       [32.0,0.04162520085409083],[33.0,0.04121351264618334],
       [34.0,0.04081771504589515],[35.0,0.04043676277279806]]

I am working with a recent build on Linux, but I tried it with 5.46.0 on Windows and got the same result.

(%i9) build_info();
(%o9) 
Maxima version: "branch_5_46_base_739_ge5422f7f0"
Maxima build date: "2022-12-12 16:25:29"
Host type: "x86_64-pc-linux-gnu"
Lisp implementation type: "SBCL"
Lisp implementation version: "2.1.11.debian"
User dir: "/home/dodier/.maxima"
Temp dir: "/tmp"
Object dir: "/home/dodier/.maxima/binary/branch_5_46_base_739_ge5422f7f0/sbcl/2_1_11_debian"
Frontend: false

wxmaxima, why my function don't integrate well?

1 Answers1