Reputation: 33
I've been trying to implement trapezoid rule for double integral. I've tried many approaches, but I can't get it to work right.
static double f(double x) {
return Math.exp(- x * x / 2);
}
// trapezoid rule
static double trapezoid(double a, double b, int N) {
double h = (b - a) / N;
double sum = 0.5 * h * (f(a) + f(b));
for (int k = 1; k < N; k++)
sum = sum + h * f(a + h*k);
return sum;
}
I understand the method for a single variable integral, but I don't know how to do it for a 2D integral, say: x + (y*y). Could someone please explain it briefly?
Upvotes: 3
Views: 2068
Reputation: 11
Consider using the class jhplot.F2D
from the DataMelt Java program. You can integrate and visualize 2D functions doing something like:
f1=F2D("x*y",-1,1,-1,1) # define in a range
print f1.integral()
Upvotes: 0
Reputation: 1543
If you're intent on using the trapezoid rule then you would do it like so:
// The example function you provided.
public double f(double x, double y) {
return x + y * y;
}
/**
* Finds the volume under the surface described by the function f(x, y) for a <= x <= b, c <= y <= d.
* Using xSegs number of segments across the x axis and ySegs number of segments across the y axis.
* @param a The lower bound of x.
* @param b The upper bound of x.
* @param c The lower bound of y.
* @param d The upper bound of y.
* @param xSegs The number of segments in the x axis.
* @param ySegs The number of segments in the y axis.
* @return The volume under f(x, y).
*/
public double trapezoidRule(double a, double b, double c, double d, int xSegs, int ySegs) {
double xSegSize = (b - a) / xSegs; // length of an x segment.
double ySegSize = (d - c) / ySegs; // length of a y segment.
double volume = 0; // volume under the surface.
for (int i = 0; i < xSegs; i++) {
for (int j = 0; j < ySegs; j++) {
double height = f(a + (xSegSize * i), c + (ySegSize * j));
height += f(a + (xSegSize * (i + 1)), c + (ySegSize * j));
height += f(a + (xSegSize * (i + 1)), c + (ySegSize * (j + 1)));
height += f(a + (xSegSize * i), c + (ySegSize * (j + 1)));
height /= 4;
// height is the average value of the corners of the current segment.
// We can use the average value since a box of this height has the same volume as the original segment shape.
// Add the volume of the box to the volume.
volume += xSegSize * ySegSize * height;
}
}
return volume;
}
Hope this helps. Feel free to ask any questions you may have about my code (warning: The code is untested).
Upvotes: 3
Reputation: 4732
Many ways to do it.
If you already know it for 1d you could make it like this:
x
That way you can basically have as many dimensions as you like. Though it scales poorly. For larger problems it might be neccesary to use monte carlo integration.
Upvotes: 0