Here's a program I made for estimating a double integral over a general region. It's pretty accurate, but it's VERY, VERY slow. (Getting an accurate enough result takes about 30 seconds)

Here's the code; the lines are indented just for readability.

Input "Surface: ",Str0
Input "y-Bound 1: ",Str0
Input "y-Bound 2: ",Str0
Input "X1: ",A
Input "X2: ",B
Input "Num of intv: ",N
DelVar DDelVar XDelVar YDelVar Str0
Disp "
Disp "Volume in R: ",V

Variables explained:

Str0: Arbitrary (input)
Y₈: Equation of surface f(x,y)
Y₉: Equation of one bound of region R
Y₀: Equation of other bound of region R
A: x-value of left bound of region R (input)
B: x-value of right bound of region R (input)
N: Number of Δx intervals (input)
D: Size of Δx interval
X: x-value of midpoint of Δx interval
V: Volume under the surface f(x,y) bound by the region R (output)

Basically, how I'm tackling it is I'm taking a slice of the region R on the x-y plane (Δx is the width) and finding its area (fnInt(abs(Y₀-Y₉),X,A+KD,A+(K+1)D)). Then I find the x-value that's the middle of that slice, and I find the average z-value along that x (fnInt(Y₈,Y,Y₀(X),Y₉(X))/(Y₀(X)-Y₉(X))), which I assume to be a reasonable estimate for the height of that entire slice of volume. Then I run a For( loop a bunch of times and add up all of the slices together to find the total volume.

But the problem is that I need lots of intervals (about 50-100) to get an accurate enough answer, which takes a ridiculous amount of time.

I already tried to optimize it a little by getting around all the expr('s I used to have, but it's still very very slow.

Any suggestions on optimizing?


Here's the important part:


I suggest you replace the fourth line with:


Here's the reason: You should have thin enough slices that the volume in a single slice of width D can be approximated by DfnInt(Y₈,Y,Y₀,Y₉); that's more important than getting the area exactly right. This change will increase your overall speed by about 3x.

Once you've done that, you can drop the For( loop entirely, and use sum(seq( instead. This only helps the performance a little; your bottleneck is still the fnInt(. However it's better style.

Another trick: since your precision is limited more by the Midpoint rule calculation along the x axis than the accurate fnInt( along the Y axis, you may want to increase the tolerance of fnInt( to, say, .005. Then your entire loop becomes:


(This uses the optional 'step' argument to seq(.

Another possibility is to use the Midpoint rule on both axes instead of fnInt(. I'm not sure how much this would help speed or decrease accuracy, but I suspect the accuracy wouldn't be too much worse than using the Midpoint rule on one axis.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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