Message Boards Message Boards

0
|
9486 Views
|
9 Replies
|
1 Total Likes
View groups...
Share
Share this post:

[?] Fastest way to compute integral in loops?

Posted 6 years ago

This got me concerned. I know that there is multiple way to avoid loops as they are so inefficient in Mathematica. But I could not find a way to avoid integration within a nested loops while saving the data to a file at the same time. An example of what I am trying to do is below:

For[i = 30, i <= 139, i = i + .1, For[j = 5, j <= 14, j = j + .1, PutAppend[NIntegrate[Func2[x, y], {x, 30, i}, {y, 5, j}], "C:\\FileToSaveTo.csv"]]]

This of course is very inefficient and takes forever. Even something simpler with one loop prove to be time consuming. For example:

For[i = 30, i <= 139, i = i + .001,  PutAppend[NIntegrate[Func2[x, y], {x, 30, i}, {y, 0, i}],  "C:\\FileToSaveTo.csv"]]

I would appreciate any alternative to what I wrote above and it is important to be able to save the output to a file.

Thanks,

POSTED BY: Abdullah Alhag
9 Replies

The table is two dimensional -- the rows vary in j and columns are changing i. You can flatten the array and get one column. I get the same answer using this technique as if I do the full integration. I suggest you start with a small example -- integrate from 30 to 31 by 1/10 and 5+1/10 to 6 by 1/10 and compare your answers and make sure everything is where you expect and modify the code to match.

If this approach is confusing, you can go back to your For[] loop and do the same idea of integrating to a (x,y) location saving the solution and then continuing the integration where you left off. I did it with some Mathematica shortcuts but you can do it either way and get a similar speedup. The loop could be integrate->return or Sow[] or Putappend solution->integrate further and add the previous solution->repeat.

Regards

POSTED BY: Neil Singer
Posted 6 years ago

So Neil, I have run the code on my function and the output came out to be not exactly what I expected. For once, I got multiple column withe the same values in each row. See attached for output. Here is the code I ran

ClearAll[integral3];
integralInY[nx_, ny_, delta_] := 
 integral3[nx, ny, delta] = 
  integral3[nx, ny - delta, delta] + 
   NIntegrate[Func2[x, y], {x, 30, nx}, {y, ny - delta, ny}]
integral3[x_, 5, 1/10] = 0;

And

AbsoluteTiming[
 Export["FileToSaveToFast.csv", 
  Table[integralInY[ii, jj, 1/10], {ii, 30, 139, 1/10}, {jj, 5 + 1/10,
     14, 1/10}]]]

My thought were that I would get the same output of running the following

For[i = 30, i <= 139, i = i + .1, 
 For[j = 5, j <= 14, j = j + .1, 
  PutAppend[NIntegrate[Func2[x, y], {x, 30, i}, {y, 5, j}], 
   "FileToSaveTo.csv"]]]

Any clue what went wrong.

Thanks, Abdullah

Attachments:
POSTED BY: Abdullah Alhag
Posted 6 years ago

Thanks a lot. I will go ahead and try playing around with it. Thanks again.

Regards, Abdullah

POSTED BY: Abdullah Alhag

Abdullah,

The extension to two dimensions is actually easier than the original, for a fixed x, you integrate in only y. Then you can change your x and integrate in y again:

ClearAll[integral3];
integralInY[nx_, ny_, delta_] := 
 integral3[nx, ny, delta] = 
  integral3[nx, ny - delta, delta] + 
   NIntegrate[func2[x, y], {x, 30, nx}, {y, ny - delta, ny}]
integral3[x_, 5, 1/10] = 0;

To call it and create an array of values:

Table[integralInY[ii, jj, 1/10], {ii, 30, 36, 1/10}, {jj, 5 + 1/10, 
  10, 1/10}]

In this example I go from x = 30 to 35 and y from 5+1/10 to 10 (you MUST start one delta increment larger than 5 based on the definition I chose for the integral. delta is 1/10 in this case. If you want more resolution you need to change the delta in the definition

POSTED BY: Neil Singer

Abdullah,

Your problem is the line you highlight.

 integral2[30, 30, 1] = 0;

should be

 integral2[30, 30, 1/10] = 0;

This line is the starting point. I created a recursive algorithm that starts at {30,30} and takes steps with the integral. By defining it with a delay of 1 instead of 1/10 you never trigger it. One way around having to remember to change this when you change your resolution (delta) is to remove delta from the function call and make delat global.

I don't think your change to make it work for i not equal j is correct. I'll look at that next...

Regards

POSTED BY: Neil Singer
Posted 6 years ago

Neil, thanks for being so responsive. First, I would like you to see the attached file for the output of running your second segment of code. Specifically the following

ClearAll[integral2];
integral2[nx_, ny_, delta_] := 
 integral2[nx, ny, delta] = 
  integral2[nx - delta, ny - delta, delta] + 
   NIntegrate[func2[x, y], {x, 30, nx - delta}, {y, ny - delta, ny}] +
    NIntegrate[func2[x, y], {x, nx - delta, nx}, {y, 0, ny}]
integral2[30, 30, 1] = 0;

And

AbsoluteTiming[
 Export["FileToSaveToFast.csv", 
  Table[integral2[i, i, 1/10], {i, 30, 139, 1/10}]]]

The output came to be the same for both my and your function. I am not sure how were you able to get a valid output, setting delta to 1/10 would always yield a similar output.

Second, regarding implementing your function for different x and y, below are my thoughts

ClearAll[integral2];
integral2[nx_, ny_, delta_] := 
 integral2[nx, ny, delta] = 
  integral2[nx - delta, ny - delta, delta] + 
   NIntegrate[func2[x, y], {x, 30, nx - delta}, {y, ny - delta, ny}] +
    NIntegrate[func2[x, y], {x, nx - delta, nx}, {y, 5, ny - delta}]
integral2[30, 30, 1] = 0;

And

AbsoluteTiming[
 Export["FileToSaveToFast.csv", 
  Table[integral2[i, j, 1/10], {i, 30, 139, 1/10},  {j, 5, 14, 1/10}]]]

Needles to say, this incorrect but I am not really sure for example what this line is for

 integral2[30, 30, 1] = 0;

Regards, Abdullah

Attachments:
POSTED BY: Abdullah Alhag

Abdullah,

I only did 10 entries -- from 30 to 40 by 1's. You can extend the algorithm to any size. The good news is that you can always do the integral in "chunks" because the program caches the intermediate values so you lose little by doing one range and then another.

If you want to do smaller intervals you will need to do exact values and not floats -- the machine precision will defeat the caching approach.

In my example I must use exact fractions (1/10 and "30", not "30.",etc.)

AbsoluteTiming[
 Export["FileToSaveToFast.csv", 
  Table[integral2[i, i, 1/10], {i, 30, 139, 1/10}]]]

This takes 97 seconds to do 1091 values of this integral from 30 to 139. You can divide it again (using 1/100) but at roughly 15-20 minutes you can go get lunch...

To do your nested integral you must modify my integral2 function to handle i not equal to j. The same concept applies. Start with a known integral which is delta away from your answer and integrate in both directions over the delta (1/10 or 1/100 in your case) and add the known answer to your new integrals.

Do you understand my integral2[] function? if so can you modify it to handle i not equal to j?

Regards

Neil

POSTED BY: Neil Singer
Posted 6 years ago

Thanks Neil. That is quite the improvement from what I had before. Would you mind however explaining your code a little bit especially how it could be adjusted it so the range of x is between 30 and 139 and the range of y is between 5 and 14 while delta remain same for both, say 0.001. When I run your code, I only get about 10 entries.

Here is what I naively tried doing but the output came out to be a mess

ClearAll[integral2];
integral2[nx_, ny_, delta_] := 
 integral2[nx, ny, delta] = 
  integral2[nx - delta, ny - delta, delta] + 
   NIntegrate[Func2[x, y], {x, 30, nx - delta}, {y, ny - delta, ny}] +
    NIntegrate[Func2[x, y], {x, nx - delta, nx}, {y, 5, ny}]
integral2[30, 30, .001] = 0;
AbsoluteTiming[
 Export["FileToSaveTo.csv", Table[integral2[i, 30, 139], {j, 5, 14}]]]

Thanks,

POSTED BY: Abdullah Alhag

Abdullah,

The more "Mathematica" way to do this is with a table and to do the write at the end rather than append to the file.

First, Lets define a function:

func2[x_Real, y_Real] := Sin[x]*Sin[y]

I use the _Real to make sure that Mathematica does not first do a symbolic integral.

your way:

AbsoluteTiming[
 For[i = 30, i <= 40, i = i + 1, 
  PutAppend[NIntegrate[func2[x, y], {x, 30, i}, {y, 0, i}], 
   "FileToSaveTo.csv"]]]

I get about 19 seconds

The table approach and write at the end:

AbsoluteTiming[
 Export["FileToSaveTo2.csv", 
  Table[NIntegrate[func2[x, y], {x, 30, i}, {y, 0, i}], {i, 30, 40}]]]

I get about 18 seconds for a modest savings.(although the savings will grow with problem size)

To really get some savings in time you should realize that you are computing the same areas over and over again (integrating from 0 to 30, then 0 to 31, etc.). If we save our most recent calculation and add on the small additional piece, we can get extreme computation savings. I make a function integral2[x,y,delta] which is my integral at {x,y} and I pass delta so the function knows my step size (in my example case, delta =1 )

ClearAll[integral2];
integral2[nx_, ny_, delta_] := 
 integral2[nx, ny, delta] = 
  integral2[nx - delta, ny - delta, delta] + 
   NIntegrate[func2[x, y], {x, 30, nx - delta}, {y, ny - delta, ny}] +
    NIntegrate[func2[x, y], {x, nx - delta, nx}, {y, 0, ny}]
integral2[30, 30, 1] = 0;

The integral2[] function saves the last answer by assigning it to integral2[] using this construct in the definition:

integral2[nx, ny, delta] = ...

Using this approach we solve the problem and write the answer to a file in 0.9 seconds!!

AbsoluteTiming[
 Export["FileToSaveToFast.csv", 
  Table[integral2[i, i, 1], {i, 31, 40}]]]

and because all the answers are now cached, running it a second time does the computation in 0.002 seconds! to see the cached answers type this:

?integral2

I hope this helps,

Regards,

Neil

POSTED BY: Neil Singer
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract