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.
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