wxMaxima中高阶k值下jk_prob1函数积分异常问题排查

问题概述 #

以下Maxima代码用于计算给定k对应的p值：jk_hermite生成Hermite多项式，jk_prob1通过两次积分计算目标概率值。根据参考文献，当k<27时计算结果正常，符合持续递减的预期；但k≥27时，结果出现振荡甚至负值，疑似int2积分环节出现问题。

核心代码 #

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
);

k=0到35的计算结果 #

(%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]]

运行环境 #

(%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")

问题根源分析 #

这个问题并非单纯机器算术误差，而是符号积分器对高次振荡函数的处理缺陷，叠加不必要的数值积分误差共同导致：

1. 归一化系数的数值积分引入误差：int1对应的是Hermite函数的模平方积分，有精确解析解sqrt(π)*2^k*k!，用数值积分计算反而会在k增大时累积误差，放大后续计算的不稳定性。
2. 高次Hermite多项式的振荡特性：k≥27时，Hermite多项式阶数极高，在区间[-sqrt(2k+1), sqrt(2k+1)]附近会出现剧烈振荡（龙格现象），Maxima默认的integrate函数对这类高振荡函数的数值积分处理能力不足，容易出现积分结果偏差，甚至让int2计算值大于1，导致p=1-int2出现负值。

修复建议 #

1. 替换归一化系数的数值积分为解析解：直接用Hermite多项式的正交性公式计算int1，避免数值误差。
2. 改用自适应数值积分函数：用Maxima专门处理振荡积分的quad_qags或quad_qag替代默认integrate，提升积分稳定性。

修改后的示例代码：

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

jk_prob1(k):=block( [],
    // 用解析解计算归一化积分，避免数值误差
    int1: sqrt(%pi)*2^k*factorial(k),
    b: sqrt(1/int1),
    // 用自适应数值积分计算int2，处理高振荡函数
    int2_result: quad_qags(abs(jk_hermite(k,x)*b*exp(-x^2/2))^2,x,-sqrt(2*k+1),sqrt(2*k+1)),
    // 取积分结果（quad_qags返回[积分值, 误差估计]）
    p: 1 - int2_result[1]
);

内容的提问来源于stack exchange，提问作者john