10

Bug introduced in 12.3.1 or earlier and fixed in 13.1.0

This is a problem comes out in the discussion here and I think it's worth posting a new question for it. Consider the following sample:

With[{ψ = ψ[x, t]}, kge = D[ψ, {t, 2}] - D[ψ, {x, 2}] + ψ;
                    ic = {ψ == Sin[Pi x], D[ψ, t] == I Sin[π x]} /. t -> 0;
                    bc = ψ == 0 /. {{x -> 0}, {x -> 1}}];    
sol = NDSolveValue[{kge == 0, bc, ic}, ψ, {x, 0, 1}, {t, 0, 10}, 
   Method -> {MethodOfLines, SpatialDiscretization -> FiniteElement}];

dsol = Derivative[0, 1][sol];
J = sol[#1, #2] + dsol[#1, #2] &;


help[x_?NumericQ, t_] := Abs@J[x, t]
ListPlot[Quiet@Table[NIntegrate[help[x, t], {x, 0, 1}], {t, 0, 1, 1/20}]]

Mathematica graphics

ListPlot[Quiet@Table[NIntegrate[Abs@J[x, t], {x, 0, 1}], {t, 0, 1, 1/20}]]

Mathematica graphics

Integration based on trapezoid rule suggests the first result is correct:

ListPlot[
  Table[With[{dh = 1/20}, 
    dh (Total@# - 1/2 (First@# - Last@#)) &@Table[Abs@J[x, t], {x, 0, 1, dh}]], 
       {t, 0, 1, 1/20}]]

Why doesn't NIntegrate work properly in 2nd case? Is _?NumericQ the only possible fix here?

Tested in v12.3.1 and v13.0.0 (Wolfram Cloud).

user21
  • 39,710
  • 8
  • 110
  • 167
xzczd
  • 65,995
  • 9
  • 163
  • 468
  • 1
    This looks like a bug, try to swap t and x. t first, x second. Is that what you expect. Also, increase the AccuracyGoal -> 8, PrecisionGoal -> 8 or reduce the MaxStepSize. – user21 Feb 24 '22 at 09:00
  • 2
    I have just merged a fix for this. The fix will be available in version 13.1. Thanks for letting me know. – user21 Feb 28 '22 at 16:14

1 Answers1

5

This is a bad bug of NIntegrate. Both intergrations use the same data as you can see using EvaluationMonitor. Not to produce too many data I only integrate a small piece:

NIntegrate[help[x, 1], {x, 0, .001}, EvaluationMonitor :> Print[{x, help[x, 1]}]]
NIntegrate[Abs@J[x, 1], {x, 0, .001}, EvaluationMonitor :> Print[{x, Abs@J[x, 1]}]]

The output from EvaluationMonitor is the same in both cases:

enter image description here

However, the first result is 1.05095*10^-6 and the second: 2.66313*10^-6. It is easy to verify that the first is correct. Therefore, something must go wrong in adding up the pieces.

Please report this to support@wolfram.com

AsukaMinato
  • 9,758
  • 1
  • 14
  • 40
Daniel Huber
  • 51,463
  • 1
  • 23
  • 57
  • I doubt this is a bug in NIntegrate. – user21 Feb 24 '22 at 10:03
  • Please explain. – Daniel Huber Feb 24 '22 at 10:05
  • The interpolating function already returns wrong results. – user21 Feb 24 '22 at 10:09
  • In this case the EvaluationMonitor would not display the same values. – Daniel Huber Feb 24 '22 at 10:34
  • "Both intergrations use the same data…" This statement isn't accurate, I'm afraid. The expression in EvaluationMonitor is actually arbitrary. Try e.g. NIntegrate[help[x, 1], {x, 0, .001}, EvaluationMonitor :> Print[{x, Abs@J[x, 1]}]] or even NIntegrate[help[x, 1], {x, 0, .001}, EvaluationMonitor :> Print[{x, Sin[x]}]]. Of course in both cases NIntegrate has sampled at same xs, but this doesn't mean the integrands are the same. – xzczd Feb 25 '22 at 13:16