2
$\begingroup$

How do I evaluate the integral c?

u = -x^2 + 10 x
phiD6[x_] = WaveletPhi[DaubechiesWavelet[6], x]

c = Integrate[u phiD6[x], {x, 0, 10}]
$\endgroup$
  • $\begingroup$ Thank you Prof. José Antonio Díaz Navas. This is the perfect solution. $\endgroup$ – DHIRAJ UPADHYAY Dept of Mechan Feb 20 '18 at 13:12
6
$\begingroup$

Here's how to get Integrate[] to work.

As you may have already noticed, the result of WaveletPhi[DaubechiesWavelet[6], x] is a Piecewise[] function containing an InterpolatingFunction[]. Since the integration interval is well within the domain of validity of the InterpolatingFunction[], we can just extract that part and work with it:

dw6 = Head[Simplify[WaveletPhi[DaubechiesWavelet[6], x], 0 <= x <= 11]];

Let's check something out:

dw6["InterpolationOrder"]
   {1}

which tells us that the InterpolatingFunction[] is a piecewise linear interpolant, which makes things very easy. We just need to extract the points first:

pts = Transpose[{Flatten[dw6["Grid"]], dw6["ValuesOnGrid"]}];

From there, the strategy is to convert this into an explicit Piecewise[] object, which Integrate[] is equipped to handle:

pw[u_] = Piecewise[{InterpolatingPolynomial[#, u], #[[1, 1]] <= u < #[[2, 1]]} &
                    /@ Partition[pts, 2, 1]];

Now, we can compute the integral:

Integrate[(10 x - x^2) pw[x], {x, 0, 10}]
   11.911233499644986

This takes a while, tho. What we can do instead is to directly integrate over the pieces, and then add everything up at the end:

Total[Integrate[(10 x - x^2) InterpolatingPolynomial[#, x], {x, #[[1, 1]], #[[2, 1]]}] & /@
      Partition[Select[pts, 0 <= #[[1]] <= 10 &], 2, 1], 
      Method -> "CompensatedSummation"]
   11.911233499644995

Or, we can directly construct the closed-form expression for the integral of a single piece:

quad[{{x0_, y0_}, {x1_, y1_}}] = 
Simplify[Integrate[(10 x - x^2) InterpolatingPolynomial[{{x0, y0}, {x1, y1}}, x],
                   {x, x0, x1}]]

Total[quad /@ Partition[Select[pts, 0 <= #[[1]] <= 10 &], 2, 1],
      Method -> "CompensatedSummation"]
   11.911233499644984

which turns out to be the fastest method of all.


If, after all that, you still want to use NIntegrate[], the right Method setting to use is "InterpolationPointsSubdivision", which will split the InterpolatingFunction[] at its grid points, and integrate over all the pieces. Let's check something out first:

Length[pts] - 1
   2816

The reason I checked the number of points is that we need to set an option that will ensure the integration is necessarily done piece-by-piece:

Options[NIntegrate`InterpolationPointsSubdivision]
   {"MaxSubregions" -> 1000, "Method" -> "GlobalAdaptive", "SymbolicProcessing" -> 5}

So, we need to set "MaxSubregions" to be as high as the number of pieces, or greater. In addition, we can turn off "SymbolicProcessing", since it's a polynomial at each piece, and we can use a quadrature method that is exact for polynomials of degree $2+1=3$. With all these considerations:

NIntegrate[(10 x - x^2) dw6[x], {x, 0, 10}, 
           Method -> {"InterpolationPointsSubdivision",
                      "MaxSubregions" -> 3000,
                      Method -> {"GaussKronrodRule", "Points" -> 2}, 
                      "SymbolicProcessing" -> 0}]
   11.911233499644988

which is more accurate than the result from brute-forcing it with "GlobalAdaptive".

$\endgroup$
5
$\begingroup$

Try it with NIntegrate, as you want to obtain a numeric result. Further, phiD6[x]is a Piecewise function:

u = -x^2 + 10 x;
phiD6[x_] := WaveletPhi[DaubechiesWavelet[6], x];

NIntegrate[u phiD6[x], {x, 0, 10}, Method -> "GaussKronrodRule", AccuracyGoal -> 6, 
           MaxRecursion -> 20, WorkingPrecision -> 10]

(* 11.91115887 *)

The GaussKronrodRule, AccuracyGoal, MaxRecursion, and WorkingPrecision have been set to obtain a reliable and convergent numeric result. Moreover, an Automatic method could be used as well. However, you should maintain the other options, particularly, MaxRecursion (even with a lower value such as 10).

$\endgroup$
5
$\begingroup$

When integrating an interpolating function, I think it's simplest to just use NDSolveValue:

NDSolveValue[{g'[x] == (-x^2 + 10 x) phiD6[x], g[0] == 0}, g[10], {x, 0, 10}]

11.9112

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.