/*
* Open Source Physics software is free software as described near the bottom of this code file.
*
* For additional information and documentation on Open Source Physics please see:
* <http://www.opensourcephysics.org/>
*/
package org.opensourcephysics.numerics;
/**
* Class Integral defines various integration algorithms.
* This class cannot be subclassed or instantiated because all methods are static.
*
* @author Wolfgang Christian
*/
public final class Integral {
static final int MAX_ITERATIONS = 15;
private Integral() {} // prohibit instantiation because all methods are static
/**
* Integrates the function using the trapezoidal method.
*
* @param f the function
* @param start the first ordinate.
* @param stop the last ordinate.
* @param n the number of partitions
* @param tol relative tolerance
*
* @return the integral
*/
public static double trapezoidal(final Function f, double start, double stop, int n, final double tol) {
double step = stop-start;
int sign = (step<0) ? -1 : 1;
if(sign<1) {
step = -step;
double temp = start;
start = stop;
stop = temp;
}
int iterations = 0;
double sum = (f.evaluate(stop)+f.evaluate(start))*step*0.5;
double oldSum;
do {
oldSum = sum;
double x = start+0.5*step;
double newSum = 0;
while(x<stop) {
newSum += f.evaluate(x);
x += step;
}
sum = (step*newSum+sum)*0.5;
step *= 0.5;
iterations++;
n /= 2;
} while((n>0)||((iterations<MAX_ITERATIONS)&&(Util.relativePrecision(Math.abs(sum-oldSum), sum)>tol)));
return sign*sum;
}
/**
* Numerical integration using Simpson's rule.
*
* @param f a function.
* @param start the first ordinate.
* @param stop the last ordinate.
* @param n the number of partitions
* @return the integral
*/
public static double simpson(final Function f, final double start, final double stop, final int n) throws IllegalArgumentException {
if(n%2!=0) {
throw new IllegalArgumentException("Number of partitions must be even in Simpson's method."); //$NON-NLS-1$
}
double sumOdd = 0.0, sumEven = 0.0, x = start;
final double h = (stop-start)/(2*n);
for(int i = 0; i<n-1; i++) {
sumOdd += f.evaluate(x+h);
sumEven += f.evaluate(x+2*h);
x += 2.0*h;
}
sumOdd += f.evaluate(x+h);
return h/3.0*(f.evaluate(start)+4.0*sumOdd+2.0*sumEven+f.evaluate(stop));
}
/**
* Numerical integration using Simpson's rule.
*
* @param f the function
* @param start the first ordinate.
* @param stop the last ordinate.
* @param n minimum number of partitions
* @param tol relative tolerance
*
* @return the integral
*/
public static double simpson(final Function f, double start, double stop, int n, final double tol) {
double step = stop-start;
int sign = (step<0) ? -1 : 1;
if(sign<1) {
step = -step;
double temp = start;
start = stop;
stop = temp;
}
int iterations = 0;
double sum = (f.evaluate(stop)+f.evaluate(start))*step*0.5;
double result = sum;
double oldSum, oldResult = result;
do {
double x = start+0.5*step;
oldSum = sum;
double newSum = 0;
while(x<stop) {
newSum += f.evaluate(x);
x += step;
}
sum = (step*newSum+sum)*0.5;
step *= 0.5;
iterations++;
oldResult = result;
result = (4*sum-oldSum)/3.0;
n /= 2;
} while((n>0)||((iterations<MAX_ITERATIONS)&&(Util.relativePrecision(Math.abs(result-oldResult), result)>tol)));
return sign*result;
}
/**
* Integrates the function using Romberg's algorithm based on Richardson's deferred approach.
*
* @param f the function
* @param a
* @param b
* @param n
* @param tol tolerance
*
* @return the integral
*/
static public double romberg(Function f, double a, double b, int n, double tol) {
if(a==b) {
return(0);
}
if(tol<=0) {
return Double.NaN; // eps must be positive
}
double[] coef = new double[MAX_ITERATIONS];
double h = (b-a)/n; // starting value for step size
// first row
coef[0] = .5*(f.evaluate(a)+f.evaluate(b));
for(int k = 1; k<n; k++) {
coef[0] += f.evaluate(a+k*h);
}
coef[0] *= h;
for(int j = 1; j<MAX_ITERATIONS; j++) {
h /= 2;
double c0 = coef[0];
coef[0] = coef[j] = 0;
for(int k = 0; k<n; k++) { /* further quadrature */
coef[0] += f.evaluate(a+(2*k+1)*h);
}
coef[0] = .5*c0+h*coef[0];
int inc = 1;
for(int k = 1; k<=j; k++) {
inc *= 4;
double Lk = coef[k];
coef[k] = (inc*coef[k-1]-c0)/(inc-1);
c0 = Lk;
}
if(Util.relativePrecision(Math.abs(coef[j]-coef[j-1]), coef[j])<tol) {
// Math.abs(coef[j] - coef[j - 1]) is est error
return coef[j];
}
n *= 2;
}
return Double.NaN; /* accuracy not achieved */
}
/**
* Uses Simpson's rule to find the area of an array representing a
* function that's been evaluated at N intervals of size h, where N is
* an odd integer.
*
* Example of usage:
* int N=27;
* x=new double[N]; f=new double[N];
* double a=0, b=5, h=(b-a)/(N-1);
* for (int i=0; i< N;i++){
* x[i]=a+i*h;
* f[i]=x[i]*Math.exp(x[i]);
* }
* double sum=Simp.Simp(f,h);
*
* Results: sum=594.6615858178942
* @Author J E Hasbun 2007
*/
public static double simpson(double f[], double h) {
//Simpson's rule for numerical integration
//f is an odd array of evaluated functions in steps h
int ip = f.length; //must be an odd number
double sumOdd = 0.0, sumEven = 0.0;
for(int i = 1; i<ip-1; i += 2) {
sumOdd += f[i];
sumEven += f[i+1];
}
return(4.0*sumOdd+2.0*sumEven+f[0]-f[ip-1])*h/3.0;
}
/**
* Computes the integral of the function using an ODE solver.
*
* @param f the function
* @param start
* @param stop
* @param tol relative tolerance
*
* @return the integral
*/
static public double ode(final Function f, final double start, final double stop, final double tol) {
ODE ode = new FunctionRate(f, start);
RK45MultiStep ode_method = new RK45MultiStep(ode); // mimics a fixed size solver
// ODEMultistepSolver ode_method= new ODEMultistepSolver(ode); // mimics a fixed size solver
ode_method.setTolerance(tol);
ode_method.initialize(stop-start); // a fixed step size method
ode_method.step();
return ode.getState()[0];
}
/**
* Fills a data array with the integral of the given function.
*
* @param f Function to be integrated
* @param start double start of integral
* @param stop double end of integral
* @param tol double computation tolerance
* @param n int number of data points
* @return double[][]
*/
static public double[][] fillArray(final Function f, final double start, final double stop, final double tol, final int n) {
double[][] data = new double[2][n];
return fillArray(f, start, stop, tol, data);
}
/**
* Fills the given data array with the intgral of the given function.
* @param f Function to be integrated
* @param start double start of integral
* @param stop double end of integral
* @param tol double computation tolerance
* @param data double[][]
* @return double[][]
*/
static public double[][] fillArray(Function f, double start, double stop, double tol, double[][] data) {
ODE ode = new FunctionRate(f, start);
ODEAdaptiveSolver ode_method = new ODEMultistepSolver(ode); // must be a fixed step size algorithm
ode_method.setTolerance(tol);
double dx = 1;
int n = data[0].length;
if(n>1) {
dx = (stop-start)/(n-1);
}
ode_method.setStepSize(dx);
for(int i = 0; i<n; i++) {
data[0][i] = ode.getState()[1];
data[1][i] = ode.getState()[0];
ode_method.step();
}
return data;
}
static private final class FunctionRate implements ODE {
double state[];
Function f;
private FunctionRate(Function _f, double start) {
state = new double[2];
state[0] = 0; // integral
state[1] = start; // independent variable
f = _f;
}
public double[] getState() {
return state;
}
public void getRate(double[] state, double rate[]) {
rate[0] = f.evaluate(state[1]); // integral
rate[1] = 1; // indepenent variable
}
}
}
/*
* Open Source Physics software is free software; you can redistribute
* it and/or modify it under the terms of the GNU General Public License (GPL) as
* published by the Free Software Foundation; either version 2 of the License,
* or(at your option) any later version.
*
* Code that uses any portion of the code in the org.opensourcephysics package
* or any subpackage (subdirectory) of this package must must also be be released
* under the GNU GPL license.
*
* This software is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston MA 02111-1307 USA
* or view the license online at http://www.gnu.org/copyleft/gpl.html
*
* Copyright (c) 2007 The Open Source Physics project
* http://www.opensourcephysics.org
*/