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:


where the right result should be as below:

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

Can someone please help me figure out what I am doing wrong here? Thanks a lot!

  • 1
    $\begingroup$ 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? $\endgroup$ – Akku14 Apr 19 '20 at 19:16
  • $\begingroup$ @Akku14 Yes you are right. Thank you! Here I just posted a trial, with the actual functions I use 1000 data points. $\endgroup$ – PhysicsLab1 Apr 19 '20 at 19:23

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

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

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

In[36]:= fgxint[2]

Out[36]= 2.07613

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

Out[37]= 2.07613
  • $\begingroup$ Thanks a lot! Now it works perfectly. $\endgroup$ – 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[1] == 0},
      z, x][y] /.
    Thread[{f', f} -> Rest@NestList[Derivative[-1], fdatg, 2]
     ] // Simplify;

(*  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:


enter image description here

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[1][f][x] - 2 y Derivative[1][f][x] + 
   x y Derivative[1][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]).


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