# Integrating a function containing a wavelet

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}]

• Thank you Prof. José Antonio Díaz Navas. This is the perfect solution. Feb 20 '18 at 13:12

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[NIntegrateInterpolationPointsSubdivision]
{"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".

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).

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