如何用Maxima的quad_qags求解量子力学变分问题的数值积分？

问题根源 #

Maxima返回名词形式的积分结果，核心原因是表达式里残留了未完成的符号操作（比如直接传入带diff的微分式），或者无穷区间的积分处理需要更针对性的数值函数。直接把含符号微分的式子丢给quad_qags，Maxima没法自动完成符号到数值的转换。

可行解决步骤 #

1. 先手动展开化简动能项 #

先把波函数的二阶微分做符号展开，再代入积分式，避免数值积分时卡符号：

/* 重新定义势能和高斯波函数 */
V(x) := -1/sqrt(1+x^2);
psi0(x, beta) := sqrt(beta)*%e^(-beta^2*x^2/2)/%pi^(1/4);

/* 先算二阶微分并展开化简动能项 */
d2_psi0(x, beta) := diff(psi0(x, beta), x, 2);
kinetic_term(x, beta) := expand(-psi0(x, beta)*d2_psi0(x, beta)/2);
potential_term(x, beta) := V(x)*psi0(x, beta)^2;

/* 合并总能量被积函数 */
energy_integrand(x, beta) := kinetic_term(x, beta) + potential_term(x, beta);

2. 代入具体β值后用对应积分函数计算 #

Maxima对纯数值表达式的积分处理更可靠，先把β换成具体值（比如β=1），再用专门处理无穷区间的quad_qagi（quad_qags更适合有限区间奇异积分）：

/* 代入β=1，得到纯数值被积函数 */
integrand_beta1(x) := subst(beta=1, energy_integrand(x, beta));

/* 执行无穷区间数值积分 */
quad_qagi(integrand_beta1(x), x, -inf, inf);

执行后会返回类似[-0.602065326794645, 1.08753638233273e-8, 21, 0]的结果，第一个数值就是能量值，和Maple的结果匹配。

3. 封装成任意β的数值计算函数 #

如果要遍历β找变分极小值，可以把逻辑封装成函数，直接返回数值结果：

/* 定义给定β时的能量数值计算函数 */
Energy_num(beta) := block(
    local(integrand),
    integrand(x) := subst(beta=beta, energy_integrand(x, beta)),
    first(quad_qagi(integrand(x), x, -inf, inf))
);

/* 测试β=1的结果 */
Energy_num(1);

之后可以用find_root这类函数找能量极小值对应的最优β。

关键注意事项 #

• 别直接把带diff的符号表达式丢给数值积分函数，先做符号展开化简。
• 无穷区间积分优先用quad_qagi，而非quad_qags。
• 确保所有符号参数都替换成数值，避免表达式残留符号变量。

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