5

I have the following problem:

My c++ code can compute two functions

f1(i1,i2,i3,i4)

f2(j1,j2)

for every set of {i1,i2,i3,i4} I get some value of f1 and for every set of {j1,j2} I get some value of f2.

the sets {i1,i2,i3,i4} and {j1,j2} are given on a FIXED mesh with some constant discretisation step "h".

I need to calculate, in mathematical language, an Integral F3(x1,x3)=Integral[f1(x1,x2,x3,x4)*f2(x3,x4) dx3 dx4]

The the simple summation is not good enough, since f2 has many jumps.

Is there some c++ library which can do this sort of integration? Or some algorhithm which is easy to implement (I am not really good on c++)

many thanks

Johan Lundberg
  • 26,184
  • 12
  • 71
  • 97
Sankp
  • 51
  • 4

4 Answers4

3

If you only have the values at the mesh points and no further mathematical knowledge of the form of the curves there's nothing better you can do than the trivial summation.

There's no way else than to change the mesh or use completely other methods like http://en.wikipedia.org/wiki/Monte_Carlo_integration

Johan Lundberg
  • 26,184
  • 12
  • 71
  • 97
  • assumed I know the jumps on f2 and assumed that f1 is smooth enough. what can I do in that case? – Sankp Jan 25 '12 at 21:08
  • I meant, if for example your functions are analytical in regions, if they are step functions, splines, polynomial, etc. Then you can do exact integrals in those regions and sum the regions. Do you know the exact shape of f2 with higher accuracy than the binning? Then why can you not evaluate the functions at arbitrary locations for the integration? Are they extremely time consuming or do you get the functions from somewhere else in a tabulated form? Also have a look at the answer by @ElKamina. You would need a 2d version, like this http://math.fullerton.edu/mathews/n2003/simpsonsrule2dmod.html – Johan Lundberg Jan 25 '12 at 21:42
  • 1
    For Monte Carlo, look at the VEGAS and MISER good ol' fortran subroutines. You can perhaps find them at www.netlib.org, or as C versions in the Numerical Recipes' book. They may help cope with the fact that f2 is not regular. – Alexandre C. Jan 25 '12 at 22:04
1

You can use Simpson's rule ( http://en.wikipedia.org/wiki/Simpson%27s_rule ). But, as Johan mentioned, if f2 is steep and erratic decreasing step size h is the only solution. Another approach you might want to consider is variable h across the mesh. That is:

1. Start with a global common h
2. Divide the space into smaller subspaces
3. Calculate integral for each subspace
4. Recalculate integral for each subspace using step size h/2
5. For only subspaces where difference between integrals (h and h/2) is substantial repeat the above mentioned steps (From step 3)
ElKamina
  • 7,747
  • 28
  • 43
1

Integration is defined for functions of real arguments. Hence if you only know your functions on a fixed mesh, you need to supply an additional rule of how you define your function for arguments in between the mesh points. This really doesn't have much to do with programming, it's just math.

For example, if you know that your function is reasonably smooth, you san use linear interpolation. Of something more complicated, if you need to. But without some rule of this sort, the integration problem is simply not well defined.

Once you have such a rule -- which can only come from an underlying meaning of your functions --- you can start choosing an integration algorithm. For functions of four variables, I'd second the suggestion by Johan Lundberg to look into Monte Carlo integrators.

ev-br
  • 24,968
  • 9
  • 65
  • 78
0

You've mentioned that you know the jumps on f2, can't you separate f2 into f2 = f2a + f2b, where; f2a is a smooth function, on which, the conventional numerical integration methods will be sufficient, and f2b is a very simple function with jumps, which, you can calculate the area analytically since it's simple. Then you can just add the values since integration is a linear operation. It all depends on what you know about f2 really, I think.

enobayram
  • 4,650
  • 23
  • 36