sagemath / sage

Main repository of SageMath
https://www.sagemath.org
Other
1.48k stars 488 forks source link

Sage fails to parse libgiac output for an integral #38498

Open Newtech66 opened 3 months ago

Newtech66 commented 3 months ago

Steps To Reproduce

sage: a = var('a', domain='positive')
sage: x_0, x_1, x_2 = var('x_0 x_1 x_2')
sage: show(integrate(e^(((x_0 - x_1)^2 + (x_1 - x_2)^2)*a), x_1, -Infinity, +Infinity, hold=True)) # this displays the integral
sage: integrate(e^(((x_0 - x_1)^2 + (x_1 - x_2)^2)*a), x_1, -Infinity, +Infinity, algorithm='libgiac')

Expected Behavior

The integral in question is an approximation of the path integral of the free particle with 1 intermediate point. Giac successfully evaluates it as:

$$\frac{\sqrt{2} \sqrt{\pi} e^{\left(\frac{1}{2} a x{0}^{2} - a x{0} x{2} + \frac{1}{2} a x{2}^{2}\right)}}{2 \sqrt{-a}}$$

Actual Behavior

NotImplementedError: Unable to parse Giac output: -sqrt(pi)*erf((-2*Infinity*sqrt(-sageVARa)*sqrt(2)+sageVARx_0*sqrt(-sageVARa)*sqrt(2)+sageVARx_2*sqrt(-sageVARa)*sqrt(2))/2)*exp((sageVARa*sageVARx_0^2-2*sageVARa*sageVARx_0*sageVARx_2+sageVARa*sageVARx_2^2)/2)/(2*sqrt(-sageVARa)*sqrt(2))+sqrt(pi)*erf((2*Infinity*sqrt(-sageVARa)*sqrt(2)+sageVARx_0*sqrt(-sageVARa)*sqrt(2)+sageVARx_2*sqrt(-sageVARa)*sqrt(2))/2)*exp((sageVARa*sageVARx_0^2-2*sageVARa*sageVARx_0*sageVARx_2+sageVARa*sageVARx_2^2)/2)/(2*sqrt(-sageVARa)*sqrt(2))

Environment

- **OS**: WSL 2 running Ubuntu 22.04.4 LTS
- **Sage Version**: 10.5

Checklist

fchapoton commented 3 months ago

minimized problem, starting from the bug :

test = libgiac('erf((-Infinity*sqrt(-sageVARa)))')
test.sage()

This infinity looks strange..

fchapoton commented 3 months ago

and debugging leads to the following

RuntimeError: indeterminate expression: infinity * f(x) encountered.

fchapoton commented 3 months ago

maybe giac is using the positivity assumption on a and libgiac is not ?

Newtech66 commented 3 months ago
sage: a = var('a')
sage: integrate(e^(a*x^2), x, 0, infinity, algorithm='giac')
1/2*sqrt(pi)/sqrt(-a)
sage: integrate(e^(a*x^2), x, 0, infinity, algorithm='libgiac') # fails
(stacktrace)
NotImplementedError: Unable to parse Giac output: -sqrt(pi)*erf(-Infinity*sqrt(-sageVARa))/(2*sqrt(-sageVARa))
fchapoton commented 3 months ago

Thanks for this other simpler and similar bug. So this is not an issue about assumption.