Simple definite integral calculator

Question

I made this code to compute definite integrals in Processing. It works by getting the rectangle between the maximum value between the previous function value and the current function value, and then I calculate the area by multiplying the maximum value described before and the step. And finally, the integral is just summing all the areas.

I will do a version based on trapezoids, but that could get really slow because I would need to use the square root function witch is computationally intensive to compute.

Any suggestions are welcome.

float step = 0.0001;
int iterations = 0;

// y = x
float f1(float x) {
  return 2 * x;
}

float integral(int a, int b) {
  float sum = 0.0;
  float prevI = a;

  for (float i = a; i < b; i += step) {
    float h = max(f1(prevI), f1(i));    
    sum += step * h;

    prevI = i;
    iterations++;
  }
  return sum;
}

void setup() {
  println(integral(0, 1));
  println(iterations);
}

Answer 1

• Using the max of two consecutive values is asymmetric.
  (I'd simply use the value sampled, avoids the following issue, too:)
• For most every "sample" (i), the code computes the value of f() twice (f(prevI)): keep prevF instead
• a double (64 bits in Processing, as in Java) may be a saner type designator than float (32 bits), especially for "the accumulator" sum. When the number of summands approaches the number representable by its mantissa, you need a more sophisticated strategy than adding in numbers one by one.
• the code multiplies each and every h by step:
  don't, instead, multiply sum by step after the loop
  (similarly, don't count trips/iterations, compute as (b-a)/step)
• consider making step a parameter
  (My pet idea would be setting a time limit and adding samples until time is up.
  Depending on sampling strategy, keeping track of weights (width of the trapezoids, step) is a challenge.)

Answer 2

This is a strange way to compute an integral. For each step, you're choosing the greater of the f(x) values given the current i or the previous i to calculate the area. I'd imagine it'd be better to choose a more conventional method such as Middle Riemann Sum, which is what I believe most modern calculators use. If you want to do a trapezoidal approximation, you can simply average a Left and Right Riemann sum (which is mathematically equivalent), but that would require you to cycle from a to b twice, which is slow, but maybe more accurate than a Middle Riemann Sum.

---

Your integral function ought to be a pure function: each input has the same output, it is unaffected by state, and it does not affect state.

Right now, the integral function is affected by the state step and f1 and also affects the state iterations.

Instead, your code would be much cleaner if you made step and function arguments to integral. As for iterations, that can be easily computed by (b-a)/step.