# Accurate Integration of 2D interpolating function with weighting

I want to integrate a 2D interpolating function received from NDSolve[]. To illustrate my problem I will give a simple example:

lb = -3; rb = 3;
testFun=FunctionInterpolation[Sin[x^2] + Cos[1/y], {x, lb, rb}, {y, 1, 6}];
normalInt = Integrate[testFun[x, y] , {x, lb, rb}]
weightningInt = Integrate[testFun[x, y] x ^2, {x, lb, rb}]


I created some 2D testfunction. Integrate can deal with it an gives solution for an interpolating function. However with a weighning factor x^2 the integration is no longer carried out. This even happens if I only multiply a numerical factor 2 to testfun.

This problem can be solved by applying function interpolation again,

solutionInt = Integrate[
FunctionInterpolation[testFun[x, y] x^2, {x, lb, rb}, {y, 1, 6}][x, y], {x, lb, rb}]


but the accuracy in my case is not sufficient. Is there any way to increase the accuracy of the integration or the second interpolation? Or can I make mathematica integrate the expression with a polynomial weightning factor times an Interpolatingfunction (should be no problem with integration by parts)?

If FunctionInterpolation is only applied once (without weightning x^2), the accuracy is much higher as you can see in the plots.

(*Exact solutions*)
weightningErg = Integrate[x^2 (Sin[x^2] + Cos[1/y]), {x, lb, rb}];
normalErg = Integrate[Sin[x^2] + Cos[1/y], {x, lb, rb}];
(*Plots*)
Plot[{weightningErg/solutionInt}, {y, 1, 6}, PlotRange -> Full]
Plot[{normalErg/normalInt}, {y, 1, 6}, PlotRange -> Full]


Although the documentation is not clear on this point, FunctionInterpolation has several options,

Options[FunctionInterpolation]
(* {InterpolationOrder -> 3, InterpolationPrecision -> Automatic,
AccuracyGoal -> Automatic, PrecisionGoal -> Automatic,
InterpolationPoints -> 11, MaxRecursion -> 6} *)


The default number of InterpolationPoints is quite small, often leading to inaccurate results. For instance,

Plot[{weightningErg/solutionInt}, {y, 1, 6}, PlotRange -> Full]


calculated by the code in the question yields, which, as noted in the question, is not very good. However, if

testFun = FunctionInterpolation[Sin[x^2] + Cos[1/y], {x, lb, rb}, {y, 1, 6},
InterpolationPoints -> 110]
solutionInt = Integrate[FunctionInterpolation[testFun[x, y] x^2, {x, lb, rb}, {y, 1, 6},
InterpolationPoints -> 110][x, y], {x, lb, rb}]


are used instead, the overall accuracy is greatly improved. 