英文:
wxmaxima, why my function don't integrate well?
问题
This is a program (jk_prob1) that calculates a number (p) for a given number (k). It appears to be working fine for k < 27, according to the bibliography, but for higher values of k, there seems to be a problem with the int2 integration.
Below is a list of calculations for various values of k (from 0 to 35), along with the corresponding results:
[[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]]
You can see that the integration results start oscillating after k=27.
The issue might be related to the integration() function or potentially the arithmetic of your machine. It's challenging to determine the exact cause without further analysis. You may want to check the precision settings for the integration function and experiment with different options to see if that resolves the issue.
英文:
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.
答案1
得分: 0
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.
(%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(%
<details>
<summary>英文:</summary>
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)bexp(-x^2/2))^2,x,-sqrt(2k+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^3sqrt(%pi)erf(sqrt(3))-2sqrt(3)))
/sqrt(%pi)],
[2,
1-(%e^-5(4*%e^5sqrt(%pi)erf(sqrt(5))-44sqrt(5)))
/(4sqrt(%pi))],
[3,
1-(%e^-7*(24*%e^7sqrt(%pi)erf(sqrt(7))-1504sqrt(7)))
/(24sqrt(%pi))],
[4,
1-(%e^-9*(192*%e^9sqrt(%pi)erf(3)-217584))
/(192sqrt(%pi))],
[5,
1-(%e^-11(1920*%e^11sqrt(%pi)erf(sqrt(11))
-4548160sqrt(11)))
/(1920sqrt(%pi))],
[6,
1-(%e^-13*(23040*%e^13sqrt(%pi)erf(sqrt(13))
-351666496sqrt(13)))
/(23040sqrt(%pi))],
[7,
1-(%e^-15*(322560*%e^15sqrt(%pi)erf(sqrt(15))
-215646796815^(3/2)))
/(322560sqrt(%pi))],
[8,
1-(%e^-17*(5160960*%e^17sqrt(%pi)erf(sqrt(17))
-3450490725632sqrt(17)))
/(5160960sqrt(%pi))],
[9,
1-(%e^-19*(92897280*%e^19sqrt(%pi)erf(sqrt(19))
-418814087689216sqrt(19)))
/(92897280sqrt(%pi))],
[10,
1-(%e^-21*(1857945600*%e^21sqrt(%pi)erf(sqrt(21))
-271427643505152021^(3/2)))
/(1857945600sqrt(%pi))],
[11,
1-(%e^-23*(40874803200*%e^23sqrt(%pi)erf(sqrt(23))
-8597150358710259712sqrt(23)))
/(40874803200sqrt(%pi))],
[12,
1-(%e^-25*(980995276800*%e^25sqrt(%pi)erf(5)
-7116923021525797376000))
/(980995276800sqrt(%pi))],
[13,
1-(%e^-27(25505877196800*%e^27sqrt(%pi)erf(3^(3/2))
-10561599770612246609923^(13/2)))
/(25505877196800sqrt(%pi))],
[14,
1-(%e^-29*(714164561510400*%e^29sqrt(%pi)erf(sqrt(29))
-50060221449301592618909696sqrt(29)))
/(714164561510400sqrt(%pi))],
[15,
1-(%e^-31*(21424936845312000*%e^31sqrt(%pi)
erf(sqrt(31))
-10502935872905789447643136000sqrt(31)))
/(21424936845312000sqrt(%pi))],
[16,
1-(%e^-33*(685597979049984000*%e^33sqrt(%pi)
erf(sqrt(33))
-7147097463438018242796670156833^(3/2)))
/(685597979049984000sqrt(%pi))],
[17,
1-(%e^-35*(23310331287699456000*%e^35sqrt(%pi)
erf(sqrt(35))
-46076693732255765758560305152035^(5/2)))
/(23310331287699456000sqrt(%pi))],
[18,
1-(%e^-37*(839171926357180416000*%e^37sqrt(%pi)
erf(sqrt(37))
-143410601721679053521909949555539968
sqrt(37)))
/(839171926357180416000sqrt(%pi))],
[19,
1-(%e^-39(31888533201572855808000
%e^39sqrt(%pi)erf(sqrt(39))
-988566037943992213347770624124125184
39^(3/2)))
/(31888533201572855808000sqrt(%pi))],
[20,
1-(%e^-41(1275541328062914232320000
%e^41sqrt(%pi)erf(sqrt(41))
-10933914922854583659978082324446398382080
sqrt(41)))
/(1275541328062914232320000sqrt(%pi))],
[21,
1-(%e^-43(53572735778642397757440000
%e^43sqrt(%pi)erf(sqrt(43))
-3262278905479206540619172723651100811460608
sqrt(43)))
/(53572735778642397757440000sqrt(%pi))],
[22,
1-(%e^-45(2357200374260265501327360000
%e^45sqrt(%pi)erf(3sqrt(5))
-4903257423120959495745281094360860593225728
5^(9/2)))
/(2357200374260265501327360000sqrt(%pi))],
[23,
1-(%e^-47(108431217215972213061058560000
%e^47sqrt(%pi)erf(sqrt(47))
-334948078254017640663171030390911909706663460864
sqrt(47)))
/(108431217215972213061058560000sqrt(%pi))],
[24,
1-(%e^-49(5204698426366666226930810880000
%e^49sqrt(%pi)erf(7)
-803416542526033681542640776803851745555237602066432))
/(5204698426366666226930810880000sqrt(%pi))],
[25,
1-(%e^-51*(260234921318333311346540544000000
%e^51sqrt(%pi)erf(sqrt(51))
-804382628000720402412181989619600056578667593072640
51^(3/2)))
/(260234921318333311346540544000000sqrt(%pi))],
[26,
1-(%e^-53(13532215908553332190020108288000000
%e^53sqrt(%pi)erf(sqrt(53))
-15268863879626398704434661565976027260923048925715234816
sqrt(53)))
/(13532215908553332190020108288000000sqrt(%pi))],
[27,
1-(%e^-55(730739659061879938261085847552000000
%e^55sqrt(%pi)erf(sqrt(55))
-1953238901674385528202656418280853915490371302850560000
55^(5/2)))
/(730739659061879938261085847552000000sqrt(%pi))],
[28,
1-(%e^-57(40921420907465276542620807462912000000
%e^57sqrt(%pi)erf(sqrt(57))
-41643536862895126218478556664899501086758124485946676084736
57^(3/2)))
/(40921420907465276542620807462912000000sqrt(%pi))],
[29,
1-(%e^-59(2373442412632986039472006832848896000000
%e^59sqrt(%pi)*erf(sqrt(59))
-988655575586727665988518643617635558359966422816622541847658496
*sqrt(59)))
/(2373442412632986039472006832848896000000
sqrt(%pi))],
[30,
1-(%e^-61(142406544757979162368320409970933760000000
%e^61sqrt(%pi)*erf(sqrt(61))
-426385466109686794974521830955494084352091033957137964359263191040
*sqrt(61)))
/(142406544757979162368320409970933760000000
sqrt(%pi))],
[31,
1-(%e^-63(8829205774994708066835865418197893120000000
%e^63sqrt(%pi)erf(3sqrt(7))
-237637173843068710884180931448387629810402892660900975740616441856
*7^(9/2)))
/(8829205774994708066835865418197893120000000
sqrt(%pi))],
[32,
1-(%e^-65(565069169599661316277495386764665159680000000
%e^65sqrt(%pi)*erf(sqrt(65))
-20743919730124955449185794118204935548012355489446073210608025600000
*65^(5/2)))
/(565069169599661316277495386764665159680000000
sqrt(%pi))],
[33,
1-(%e^-67(37294565193577646874314695526467900538880000000
%e^67sqrt(%pi)*erf(sqrt(67))
-41682427353927433909745163803213551387672794288391014517753878087599128576
*sqrt(67)))
/(37294565193577646874314695526467900538880000000
sqrt(%pi))],
[34,
1-(%e^-69(2536030433163279987453399295799817236643840000000
%e^69sqrt(%pi)*erf(sqrt(69))
-296226359950091579306352646020536576938264515830526239985823226493501702144
*69^(3/2)))
/(2536030433163279987453399295799817236643840000000
sqrt(%pi))],
[35,
1-(%e^-71(177522130321429599121737950705987206565068800000000
%e^71sqrt(%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
</details>
通过集体智慧和协作来改善编程学习和解决问题的方式。致力于成为全球开发者共同参与的知识库，让每个人都能够通过互相帮助和分享经验来进步。
评论