Message Boards Message Boards

0
|
438 Views
|
4 Replies
|
3 Total Likes
View groups...
Share
Share this post:

What does TrapezoidalRule at NIntegrate really do?

Hello everybody,

I have a question as regards the TrapezoidalRule at the implemented NIntegrate function of Mathematica.

If I calculate the numeric(!) integral of f(x)=x^2 from 0 to 3 with 5 trapezoids by hand, I get a result of 9.18. If I use the implemented method "TrapezoidalRule" with the option "Points" -> 5 it gives me exactly 9.0.

So I doubt, that I do not fully understand what Mathematica does when using this Rule. I have tried a lot of (changing the WorkingPrecision, setting the MaxRecursions to 0 etc...), but it always comes up with the "real" (analytical) result.

Can anyone please explain what's going on?

Thanks a lot,

JJJ

f[x_] := x^2;
a = 0; b = 3;
\[CapitalDelta]x[n_] := (b - a)/n;
left[n_] := \!\(
\*UnderoverscriptBox[\(\[Sum]\), \(i = 0\), \(n - 1\)]\(f[
     a + i*\[CapitalDelta]x[n]]*\[CapitalDelta]x[n]\)\);
right[n_] := \!\(
\*UnderoverscriptBox[\(\[Sum]\), \(i = 1\), \(n\)]\(f[
     a + i*\[CapitalDelta]x[n]]*\[CapitalDelta]x[n]\)\);
trap[n_] := (left[n] + right[n])/2;
trap2[n_] := \[CapitalDelta]x[n]/2 (f[a] + 2 \!\(
\*UnderoverscriptBox[\(\[Sum]\), \(i = 1\), \(n - 1\)]\(f[
        a + i*\[CapitalDelta]x[n]]\)\) + f[b]);
{left[5.], right[5.], trap[5.], trap2[5.]}
NIntegrate[f[x], {x, a, b}, 
 Method -> {"TrapezoidalRule", "Points" -> 3}, WorkingPrecision -> 10,
  MaxRecursion -> 0]
POSTED BY: Joachim Janezic
4 Replies

Wow, great. I will check that tomorrow.

Thank you very much! Your support is very much appreciated!

JJJ

POSTED BY: Joachim Janezic

Do you also get the warning message which says that the Integrate failed to converge to prescribed accuracy ...?

Yes. NIntegrate[] checks an error estimate. If the error estimate is larger than the error prescribed by PrecisionGoal and AccuracyGoal, then NIntegrate[] complains. This will happen often if you set MaxRecursion to zero. If you set PrecisionGoal -> -1, AccuracyGoal -> -1, then it doesn't complain about convergence (well, hardly ever).

What do you mean by....

Sorry, that was my mistake. I was in a hurry and didn't have time to check the documentation. Ignore "interior points." For some reason, "Points" often does not mean the actual number of points used in an integration rule. It's a parameter that affects the number of points. (For instance, the Romberg step needs $2n-1$ points, and the "Points" parameter is $n$.)

If you're trying to have simple demonstrations for students, one of the difficulties to face up to is that generally speaking, Mathematica is set up to give high-quality results, using sophisticated techniques under the hood. To get it to do the naive thing usually means setting a bunch of options that have a chance of confusing the student.

Of course, your demo code is simple and clear. I suppose you're just checking it against NIntegrate[]. But there you still run into a conflict between its sophisticated approach and the desired naive approach.

You probably don't need the following, but I have the code sitting around. It makes calling Nintegrate[] simple, but the user has to remember to define myTrapRule below. It's relatively easy to define any interpolatory integration rule in Mathematica. The internal data structure used to represent the rule has the following form:

NIntegrate`GeneralRule[{abscissas, integration_weights, error_weights}]

You can add a rule to NIntegrate's repertory by extending the definition of NIntegrate`InitializeIntegrationRule[] as follows (the unexplained arguments allow for greater sophistication, which is not needed here):

Options[myTrapRule] = {"Points" -> 2};
myTrapRule /: NIntegrate`InitializeIntegrationRule[
   myTrapRule, nfs_, ranges_, ruleOpts_, allOpts_] :=
  Module[{abscissas, weights, nPoints},
   nPoints = OptionValue[myTrapRule, ruleOpts, "Points"];
   abscissas = Subdivide[0, 1, nPoints - 1];
   weights = ConstantArray[1, nPoints]/(nPoints - 1);
   weights[[{1, -1}]] /= 2;
   NIntegrate`GeneralRule[{abscissas, weights, 0*weights}]
   ];

Then:

NIntegrate[x^2, {x, 0, 3}, Method -> {myTrapRule, "Points" -> 5}]
(*  9.28125  *)

NIntegrate[x^2, {x, 0, 3}, Method -> {myTrapRule, "Points" -> 6}]
(*  9.18  *)

The error weights are set to zero on purpose. There will be no recursion because the rule reports an error of zero on the first try. NIntegrate[] does not care if you lie to it. You can redefine the meaning of "Points" as you see fit. The abscissas are obtained by subdividing the unit interval [0, 1]. NIntegrate[] rescales an integral to be over [0, 1], so you do not need to worry about other intervals. (Infinite intervals usually get special handling that does not apply to integrals for which myTrapRule is an appropriate rule.)

POSTED BY: Michael Rogers

Thanks a lot.

Two more questions please:

  • Do you also get the warning message which says that the Integrate failed to converge to prescribed accuracy ...?
  • What do you mean by "interior points"? Let's say I have this interval [0,1] and I want to split it in 5 pieces (as it was in my example). Which number do I have to put at the Points-option? 4 (which seems to be obvious for "interior points" gives a bad result (9.125 instead of 9.18). In the documentation it says: "The option "Points"->k can be used how many coarse points are used. The number of points used by "TrapezoidalRule" is 2k-1" In my understanding, This means I cannot have 5 "pieces" (with 4 interior points), because there is no integer k which doubled - 1 = 4. Is this correct?

Thank you,

JJJ

POSTED BY: Joachim Janezic

It's using extrapolation I think. Set the suboption "RombergQuadrature" -> False, and you may as well set "SymbolicProcessing" -> None just in case. It's not needed for x^2, but I'm not sure when it might be necessary. Also, your n determines the number of intervals and has n+1 sample points. The "Points" option in the "TrapezoidalRule" specifies the number of interior points. [Correction: The "Points" -> n option in the "TrapezoidalRule" specifies the number of sample points to be 2n-1.]

{left[4], right[4], trap[4], trap2[4]} // N
NIntegrate[f[x], {x, a, b}, 
  Method -> {"TrapezoidalRule", "Points" -> 3,
   "RombergQuadrature" -> False,
   "SymbolicProcessing" -> None},
  WorkingPrecision -> 10,
   MaxRecursion -> 0]

Both give 9.28125.

Reference: https://reference.wolfram.com/language/tutorial/NIntegrateIntegrationRules.html#618158740

POSTED BY: Michael Rogers
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract