Java - Simpson's method and error

2.8k views Asked by At

I am writing a Java program for Simpson's method. The basic program works as expected, although I cannot get the (absolute) error part to work.

I think I need to reference my while (absError < 0.000001) loop differently. What am I doing wrong?

First try

public static double function(double x, double s) {
    double sech = 1 / Math.cosh(x); // Hyperbolic cosecant
    double squared = Math.pow(sech, 2);
    return ((Math.pow(x, s)) * squared);
 }

 // Simpson's rule - Approximates the definite integral of f from a to b.
 public static double SimpsonsRule(double a, double b, double s, int n) {
    double dx, x, sum4x, sum2x;
    double absError = 1.0;
    double simpson = 0.0;
    double simpson2 = 0.0;

    dx = (b-a) / n;
    sum4x = 0.0;
    sum2x = 0.0;

    // 4/3 terms
    for (int i = 1; i < n; i += 2) {
        x = a + i * dx;
        sum4x += function(x,s);
    }

    // 2/3 terms
    for (int i = 2; i < n-1; i += 2) {
        x = a + i * dx;
        sum2x += function(x,s);
    }

    // Compute the integral approximation.
    simpson = function(a,s) + function(a,b);
    simpson = (dx / 3)*(simpson + 4 * sum4x + 2 * sum2x);
    while ( absError < 0.000001)
    {
        simpson2 = SimpsonsRule(a, b, s, n);
        absError = Math.abs(simpson2 - simpson) / 15;
        simpson = simpson2;
        n++;
    }
    System.out.println("Number of intervals is " + n + ".");
    return simpson2;
}

This doesn't work since I did not write

simpson2 = SimpsonsRule(a, b, s, n);

correctly.

I tried doing this a second way, but the solution ultimately also fails.

public static double function(double x, double s) {
    double sech = 1 / Math.cosh(x); // Hyperbolic cosecant
    double squared = Math.pow(sech, 2);
    return ((Math.pow(x, s)) * squared);
 }

 // Simpson's rule - Approximates the definite integral of f from a to b.
public static double SimpsonsRule(double a, double b, double s, int n) {
    double dx, x, sum4x, sum2x;
    double absError = 1.0;
    double simpson = 0.0;
    double simpson2 = 0.0;

    dx = (b-a) / n;
    sum4x = 0.0;
    sum2x = 0.0;

    // 4/3 terms
    for (int i = 1; i < n; i += 2) {
        x = a + i * dx;
        sum4x += function(x,s);
    }

    // 2/3 terms
    for (int i = 2; i < n-1; i += 2) {
        x = a + i * dx;
        sum2x += function(x,s);
    }

    // Compute the integral approximation.
    simpson = function(a,s) + function(a,b);
    simpson = (dx / 3)*(simpson + 4 * sum4x + 2 * sum2x);
    while ( absError < 0.000001)
    {
        n++;
        dx = (b-a) / n;
         // 4/3 terms
        for (int i = 1; i < n; i += 2) {
            x = a + i * dx;
            sum4x += function(x,s);
        }

        // 2/3 terms
        for (int i = 2; i < n-1; i += 2) {
            x = a + i * dx;
            sum2x += function(x,s);
        }
        simpson = function(a,s) + function(a,b);
        simpson2 = (dx / 3)*(simpson + 4 * sum4x + 2 * sum2x);
        absError = Math.abs(simpson2 - simpson) / 15;
        simpson = simpson2;
    }
    System.out.println("Number of intervals is " + n + ".");
    return simpson2;
}

I need to write the while loop differently. What is wrong with the way the error is referenced inside the while loop?

The java code up until

    while ( absError < 0.000001)
{
    simpson2 = SimpsonsRule(a, b, s, n);
    absError = Math.abs(simpson2 - simpson) / 15;
    simpson = simpson2;
    n++;
}
System.out.println("Number of intervals is " + n + ".");
return simpson2;

Works fine and correctly calculates Simpson's method.

2

There are 2 answers

0
pkozlov On BEST ANSWER

Looks like your Simpson's method implementation is not converging. The easiest thing you can do to avoid infinite cycle in while - you have to add another condition - maximum number of iterations.

Something like that:

int n = 0;
while (error < ACCURACY && n++ < MAX_ITERATIONS) {
    // while body
}

where ACCURACY is 0.000001 in your case (or 1e-6) and MAX_ITERATIONS is an integer constant, for example 100000 or 1e+6.

Why your algorithm is not converging - this is another question - look carefully on your formulas - use debugging tools. Good luck!

0
Axion004 On

I fixed it. Thanks for your help.

 // Simpson's rule - Approximates the definite integral of f from a to b.
 public static double SimpsonsRule(double a, double b, double s, int n) {
    double simpson, dx, x, sum4x, sum2x;

    dx = (b-a) / n;
    sum4x = 0.0;
    sum2x = 0.0;

    // 4/3 terms
    for (int i = 1; i < n; i += 2) {
        x = a + i * dx;
        sum4x += function(x,s);
    }

    // 2/3 terms
    for (int i = 2; i < n-1; i += 2) {
        x = a + i * dx;
        sum2x += function(x,s);
    }

    // Compute the integral approximation.
    simpson = function(a,s) + function(a,b);
    simpson = (dx / 3)*(simpson + 4 * sum4x + 2 * sum2x);
    return simpson;
}

// Handles the error for for f(x) = t^s * sech(t)^2. The integration is
// done from 0 to 100.
// Stop Simspson's Method when the relative error is less than 1 * 10^-6
public static double SimpsonError(double a, double b, double s, int n)
{
    double futureVal;
    double absError = 1.0;
    double finalValueOfN;
    double numberOfIterations = 0.0;
    double currentVal = SimpsonsRule(a,b,s,n);

    while (absError / currentVal > 0.000001) {
        n = 2*n;
        futureVal = SimpsonsRule(a,b,s,n);
        absError = Math.abs(futureVal - currentVal) / 15;
        currentVal = futureVal;
    }

    // Find the number of iterations. N starts at 8 and doubles every iteration.
    finalValueOfN = n / 8;
    while (finalValueOfN % 2 == 0) {
        finalValueOfN = finalValueOfN / 2;
        numberOfIterations++;
    }
    System.out.println("The number of iterations is " + numberOfIterations + ".");
    return currentVal;
}