# Numerical integration of an interpolating function with another analytic function [closed]

I have to numerically integrate over $$x$$ a complicated function that has this structure $$f(x,y)g(x)$$ with a parameter $$y$$. In my case the fastest way to deal with it is generate an interpolating function of $$g(x)$$ and then numerically integrate it with the rest. Though I tried to check this method on a much simpler function such as $$(x + y*x - 2*y) (x + Log[x] + Sqrt[x])$$ where in this case I call $$g(x)=(x + Log[x] + Sqrt[x]$$ and $$f(x,y)=(x + y*x - 2*y)$$ to see if it works but it actually gives wrong results if I try to cross check for a specific value of y. The integral has to be performed from $$1$$ to $$y$$ for any value of $$y$$,(I will need a plot wrt $$y$$ later)
Here it is what I am doing:

datg = Flatten[Table[{x, (x + Log[x] + Sqrt[x])}, {x, 1, 3, 0.0001}], 1];
fdatg = Interpolation@datg;
fg[x_, y_] := (x + y*x - 2*y) fdatg[x];
fgxint[y_?NumericQ] := Module[{x}, NIntegrate[fg[x, y], {x, 1, y}]];


For example for y=2 the code gives:

fgxint=1.14999


where the right result should be as below:

NIntegrate[(x + 2*x - 2*2) (x + Log[x] + Sqrt[x]), {x, 1, 2}]
out:2.07613


• Have you never tested your datg. Due to Flatten it mixes up x-values and y-values. Just leave the Flatten off. Next question; why do you need 20.001 data points? – Akku14 Apr 19 '20 at 19:16
• @Akku14 Yes you are right. Thank you! Here I just posted a trial, with the actual functions I use 1000 data points. – PhysicsLab1 Apr 19 '20 at 19:23

Your flattening of datg is wrong. This code produces agreeing results:

In:= datg = Table[{x, (x + Log[x] + Sqrt[x])}, {x, 1, 3, 0.0001}];
fdatg = Interpolation@datg;

In:= Clear[fg, fgxint];
fg[x_, y_] := (x + y*x - 2*y)*fdatg[x];
fgxint[y_?NumericQ] := NIntegrate[fg[x, y], {x, 1, y}];

In:= fgxint

Out= 2.07613

In:= NIntegrate[(x + 2*x - 2*2) (x + Log[x] + Sqrt[x]), {x, 1, 2}]

Out= 2.07613

• Thanks a lot! Now it works perfectly. – PhysicsLab1 Apr 19 '20 at 19:24

A semi-symbolic way that avoids calling NIntegrate every time fgxint[y] is executed:

datg = Table[{x, (x + Log[x] + Sqrt[x])}, {x, 1, 3, 0.0001}];
fdatg = Interpolation@datg;

ClearAll[fgxint, x, y, z, f];
fgxint[y_] = DSolveValue[
{z'[x] == (x + y*x - 2*y) f''[x], z == 0},
z, x][y] /.
Thread[{f', f} -> Rest@NestList[Derivative[-1], fdatg, 2]
] // Simplify;

fgxint
(*  2.07613  *)


Derivative[-1][fdatg] computes the antiderivative of the interpolating function fdatg and returns another interpolating function. This is done symbolically, since an InterpolatingFunction is a piecewise polynomial function and so is its antiderivative. The result is an expression in terms of interpolating function:

fgxint[y] This should be much more efficient than calling NIntegrate.

The other tricky bit is getting Integrate/DSolve to use integration by parts. This is why f'' is used to stand for fdatg. For example, using Integrate:

Integrate[(x + y*x - 2*y) f''[x], x]
(*
-f[x] - y f[x] + x Derivative[f][x] - 2 y Derivative[f][x] +
x y Derivative[f][x]
*)


If f'' === fdatg, then f' is the integral (Derivative[-1]) of fdatg and f is the integral of that (Derivative[-2][fdatg] == Derivative[-1]Derivative[-1][fdatg]).