Wednesday, August 12, 2015

Integral Approximation with Infinite Limits and Regular Partitions




Problem
  
There is a theorem in integral calculus that states that if a function f is integrable on [a, b], then its definite integral on [a, b] is equal to an infinite limit of the sum of f evaluated at regular points obtained with ever decreasing regular partitions. In other words:



Let us investigate how we can use this theorem to programmatically evaluate definite integrals.

Function Abstraction

We need to use functions as objects. Toward that result, we define a one-argument Function class that all other subclasses will override. Source code is here.

public class Function {
   
    public Function() {}
    public double v(double x) { return 0; }
    public double[] generateRangeInterval(double[] domain_values) {
        return null;
    }
   
}


The idea of regular partition can now be encapsulated as a subclass of the Function class (RegularPartition.java):

public class RegularPartitionSum extends Function {
   
    private double mA = 0;
    private double mB = 0;
    private Function mF = null;
   
    public RegularPartitionSum(Function f, double a, double b) {
        mA = a;
        mB = b;
        mF = f;
    }
   
    @Override
    public double v(double n) {
        final double step = (mB - mA)/n;
        final int in = (int)n;
        double sum = 0;
        for(int i = 0; i <= in; i++) {
            sum += mF.v(mA + i*step);
        }
        return step*sum;
    }
   
}


Given a function, we should be able to approximate its limit at infinity to any desirable number of steps, which we abstract into a FunctionLimit class (FunctionLimit.java).

public class FunctionLimit {

public static double limitAtInfinity(Function f, double start, double step, int num_steps) {
        double x = start;
        double y = 0;
        while ( num_steps >= 0 ) {
            y = f.v(x);
            x += step;
            num_steps--;
        }
        return y;

}

We proceed by implementing Riemann's sum (RiemannSum.java) with the method that uses FunctionLimit to estimate a definite integral with regular partitions and infinite limits. We compute the limit of the regular partition sums for a specific number of steps.

public class RiemannSum {
    
    public static double infiniteLimitOfRegularPartition(Function f, double a, double b, int num_steps) {
        Function rpF = new RegularPartitionSum(f, a, b);
        return FunctionLimit.limitAtInfinity(rpF, 1, 1, num_steps);
    }
}



Sample Functions

Let us define several functions and use the above implementation to evaluate their definite integrals. All of them come from James Stewart's "Calculus: Early Transendentals," 3rd edition. The name of the function class specifies the number of example and the page.

public class F_example_07_p343 extends Function {
   
    public F_example_07_p343() {}
   
    @Override
    public double v(double x) {
        return Math.sqrt(x);
    }
}


The function example_07_p343(), defined below, when placed in the main method, outputs "RS = 4.668915495400732," which is pretty close to the manual evaluation of the definite intergal.
 
public static void example_07_p343() {
        Function f = new F_example_07_p343();
        System.out.println("RS  = " + RiemannSum.infiniteLimitOfRegularPartition(f, 1, 4, 2000));

}

Here is another sample function:

public class F_ex09_p344 extends Function {
   
    public F_ex09_p344() {}
   
    @Override
    public double v(double x) {
        return x*x*x;
    }
}


We can evaluate this function in this method:

public static void ex_09_p344() {
        System.out.println("RS = " + RiemannSum.infiniteLimitOfRegularPartition(new F_ex09_p344(), 0, 5, 1000));

}

The output value is 156.56234375015595. More sample functions are implemented and tested here.
  

Midpoint Rule

The accuracy of our evaluations can be compared with the mid-point rule, which is a more standard method of evaluating definite integrals. Add these methods to RiemannSum and use them to evaluate the sample functions.

public static double midPointArea(Function f, double from, double upto, double step) {
        double area = 0;
        if ( from >= upto ) return 0;
        double currP = from;
        while ( currP <= upto - step ) {
            area += f.v(currP+step/2.0);
            currP += step;
        }
        return step*area;

}

public static double midPointRule(Function f, double a, double b, int n) {
        if ( b < a ) throw new IllegalArgumentException("midPointRule: b < a");
        if ( a == b ) return 0;
        final double step = (b - a)/n;
        return RiemannSum.midPointArea(f, a, b, step);

}