Integral changes sign when changing log(%e^x^2/2) to x^2-log(2)

[rtoy]

Imported from SourceForge on 2024-07-04 01:50:46 Created by tomasriker on 2022-06-03 04:05:18 Original: https://sourceforge.net/p/maxima/bugs/3988

---

This integral has the wrong sign:

(%i1) float(integrate(%e^-x^2*log(%e^x^2/2), x, 0, inf));
(%o1) 0.1711722319875089

When we change log(%e^x^2/2) to the equivalent x^2-log(2), the integral is correct:

(%i2) float(integrate(%e^-x^2*(x^2-log(2)), x, 0, inf));
(%o2) -0.171172231987509

[rtoy]

Imported from SourceForge on 2024-07-04 01:50:47 Created by robert_dodier on 2022-06-03 13:33:44 Original: https://sourceforge.net/p/maxima/bugs/3988/#5516

---

trace(?rischint) shows that RISCHINT (Risch integrator) is called for %i1 but not %i2. Not sure what to make of it.

[rtoy]

Imported from SourceForge on 2024-07-04 01:50:51 Created by macrakis on 2022-06-03 14:52:28 Original: https://sourceforge.net/p/maxima/bugs/3988/#c99d

---

See #3989 for a possibly related bug in limit -- a simplified form of the antiderivative of the integrand.

[rtoy]

Imported from SourceForge on 2024-07-04 01:50:55 Created by robert_dodier on 2022-06-04 15:54:47 Original: https://sourceforge.net/p/maxima/bugs/3988/#7299

---

Here is some evidence, not conclusive, that the Risch code (via ANTIDERIV) is coming up with a valid antiderivative, so it's other stuff (either in defint or in limit) where the problem lies.

Integrand in question:

(%i2) ee: %e^-x^2*log(%e^x^2/2);
                                      2
                             2       x
                          - x      %e
(%o2)                   %e     log(----)
                                    2

When the lower limit of integration is greater than zero, defint seems to get a correct result (not shown except for lower limit = 1). First we see that we get the same result from ANTIDERIV whether lower limit = 0 or 1.

(%i3) trace (?antideriv) $
(%i4) display2d: false $
(%i5) integrate (ee, x, 0, inf);
1 Enter antideriv [(yx*%e^-yx)/sqrt(log(2*%e^yx))]
1 Exit  antideriv false
1 Enter antideriv [log(yx)/(yx^2*sqrt(log(2*yx)))]
1 Exit  antideriv 
(-(2*sqrt(%pi)*log(2*yx)-sqrt(%pi))*erf(sqrt(log(2*yx))))
 +2*sqrt(%pi)*log(yx)*erf(sqrt(log(2*yx)))-sqrt(log(2*yx))/yx
(%o5) -(sqrt(%pi)-2*sqrt(%pi)*log(2))/4
(%i6) integrate (ee, x, 1, inf);
1 Enter antideriv [log(yx)/(yx^2*sqrt(log(2*yx)))]
1 Exit  antideriv 
(-(2*sqrt(%pi)*log(2*yx)-sqrt(%pi))*erf(sqrt(log(2*yx))))
 +2*sqrt(%pi)*log(yx)*erf(sqrt(log(2*yx)))-sqrt(log(2*yx))/yx
(%o6) ((-2*sqrt(%pi)*erf(1)*log(%e/2))
 -2*sqrt(%pi)*log(2)+sqrt(%pi)*erf(1)+sqrt(%pi)+2*%e^-1)
 /4
(%i7) ee1a: (-(2*sqrt(%pi)*log(2*yx)-sqrt(%pi))*erf(sqrt(log(2*yx))))
 +2*sqrt(%pi)*log(yx)*erf(sqrt(log(2*yx)))-sqrt(log(2*yx))/yx $
(%i8) ee1b: (-(2*sqrt(%pi)*log(2*yx)-sqrt(%pi))*erf(sqrt(log(2*yx))))
 +2*sqrt(%pi)*log(yx)*erf(sqrt(log(2*yx)))-sqrt(log(2*yx))/yx $
(%i9) is (ee1a = ee1b);
(%o9) true

Now verify that simplifying the integrand yields the same numerical result:

(%i10) integrate (ev (ee, logexpand = super), x, 1, inf);
1 Enter antideriv [(x^2-log(2))*%e^-x^2]
1 Exit  antideriv 
(-(sqrt(%pi)*log(2)*erf(x))/2)+(sqrt(%pi)*erf(x))/4
                              -(x*%e^-x^2)/2
(%o10) (-(sqrt(%pi)*(2*log(2)-1))/4)
 +(sqrt(%pi)*erf(1)*log(2))/2-(sqrt(%pi)*erf(1))/4+%e^-1/2
(%i11) [%o6, %o10], numer;
(%o11) [0.1570144642250585,0.1570144642250585]

while the form of the expression is different (this is okay):

(%i12) is (%o6 = %o10);
(%o12) false

and the numerical result for lower limit = 1 agrees with quadrature applied to the simplified integrand:

(%i13) quad_qagi (ev (ee, logexpand = super), x, 1, inf);
(%o13) [0.1570144642250585,1.527495796238393e-11,165,0]

We have some evidence that ANTIDERIV came up with a valid result, but it's not as simple as just plugging it into the fundamental theorem of calculus; defint must be doing something at least a little more complicated:

(%i16) limit (ee1b, yx, inf) - ev (ee1b, yx = 1);
(%o16) (-(sqrt(%pi)-2*sqrt(%pi)*log(2))*erf(sqrt(log(2))))
 -2*sqrt(%pi)*log(2)+sqrt(log(2))+sqrt(%pi)
(%i17) %, numer;
(%o17) 0.6688921216552621

which doesn't match %o11.

[rtoy]

Imported from SourceForge on 2024-07-04 01:50:58 Created by robert_dodier on 2022-06-04 20:27:25 Original: https://sourceforge.net/p/maxima/bugs/3988/#4924

---

On looking at this some more, my working hypothesis at the moment is that the source of the problem is (maybe) actually farther upstream from ANTIDERIV and limit, namely in INTCV which makes a change of variable and then calls DEFINT again. It seems possible that the sign got switched by INTCV, and then RMCONST1 pulls out -1/4 instead of 1/4, so when INITIAL-ANALYSIS returns something, the final result has the sign switched.

Need to look more carefully at INTCV, DINTLOG (which calls INTCV), DEFINT (called by INTCV), and RMCONST1 (called by DEFINT).

[rtoy]

Imported from SourceForge on 2024-07-04 01:51:01 Created by robert_dodier on 2022-06-05 16:23:36 Original: https://sourceforge.net/p/maxima/bugs/3988/#bf58

---

I think I see it. The problem appears to be that in doing the change of variable exp(x^2)/2 --> yx, INTCV chooses the wrong branch when the lower limit is 0, and otherwise chooses the correct branch. Any change of variable which is not one to one is probably susceptible to an analogous bug.

Let ee be the integrand in question. INTCV calls TEST-INVERSE to test potential solutions found by SOLVE.

Lower limit = 1/10, INTCV chooses correct branch:

(%i37) integrate(ee, x, 1/100, inf);
                  2
                 x
               %e
1 Enter intcv [----, false]
                2
                           2
                          x
                        %e                               1
 1 Enter test\-inverse [----, x, - sqrt(log(2 yx)), yx, ---]
                         2                              100
 1 Exit  test\-inverse false
                           2
                          x
                        %e                             1
 1 Enter test\-inverse [----, x, sqrt(log(2 yx)), yx, ---]
                         2                            100
 1 Exit  test\-inverse true
                                             1/10000
                                   1       %e
1 Exit  intcv ((- 2 sqrt(%pi) erf(---) log(---------))
                                  100          2
                                            1
                        4999 sqrt(%pi) erf(---)
                                           100
 - 2 sqrt(%pi) log(2) - ----------------------- + sqrt(%pi)
                                 5000
     - 1/10000
   %e
 + -----------)/4
       50
                                      1/10000
                            1       %e
(%o37) ((- 2 sqrt(%pi) erf(---) log(---------))
                           100          2
                                            1
                        4999 sqrt(%pi) erf(---)
                                           100
 - 2 sqrt(%pi) log(2) - ----------------------- + sqrt(%pi)
                                 5000
     - 1/10000
   %e
 + -----------)/4
       50
(%i38) %, numer;
(%o38)                - 0.1642413245373725

Lower limit = 0, INTCV chooses wrong branch:

(%i39) integrate(ee, x, 0, inf);
                  2
                 x
               %e
1 Enter intcv [----, false]
                2
                           2
                          x
                        %e
 1 Enter test\-inverse [----, x, - sqrt(log(2 yx)), yx, 0]
                         2
 1 Exit  test\-inverse true
                sqrt(%pi) - 2 sqrt(%pi) log(2)
1 Exit  intcv - ------------------------------
                              4
                  sqrt(%pi) - 2 sqrt(%pi) log(2)
(%o39)          - ------------------------------
                                4
(%i40) %, numer;
(%o40)                  0.171172231987509

Not sure how to detect correct branch in this case or in general.

Incidentally changevar seems to choose the wrong branch whether the lower limit is 0 or nonzero. I'll open a bug report about that.