MidPointMatrixQuadrature.java
/*
* Copyright (C) 2023 Alberto Irurueta Carro (alberto@irurueta.com)
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package com.irurueta.numerical.integration;
import com.irurueta.algebra.Matrix;
import com.irurueta.algebra.WrongSizeException;
import com.irurueta.numerical.EvaluationException;
/**
* Implementation of matrix quadrature using mid-point algorithm.
* Mid-point algorithm is suitable for improper integrations, which consists of the following
* situations:
* - integrand goes to a finite limiting value at finite upper and lower limits, but cannot be
* evaluated right on one of those limits (e.g. sin(x)/x at x = 0)
* - the upper limit is positive infinity, or the lower limit is negative infinity.
* - integrand has an integrable singularity at either limit (e.g. x^-0.5 at x = 0)
* - integrand has an integrable singularity at a known place between its upper and lower limits.
* - integrand has an integrable singularity at an unknown place between its upper and lower limits.
* <p>
* Impossible integrals, such as those being infinite (e.g. integral from 1 to infinity of x^-1) or
* those that have no limiting sense (e.g. integral from -infinity to infinity of cos(x)) cannot
* be handled by this implementation.
*/
public class MidPointMatrixQuadrature extends MatrixQuadrature {
/**
* Lower limit of integration.
*/
private final double a;
/**
* Upper limit of integration.
*/
private final double b;
/**
* Number of rows of quadrature result.
*/
private final int rows;
/**
* Number of columns of quadrature result.
*/
private final int columns;
/**
* Current value of integral.
*/
private final Matrix s;
/**
* Listener to evaluate single dimension matrix functions at required points.
*/
protected final MatrixSingleDimensionFunctionEvaluatorListener listener;
/**
* Temporary value storing evaluation at mid-point.
*/
private final Matrix tmpMid;
/**
* Temporary value storing summation of evaluations.
*/
private final Matrix sum;
/**
* Temporary value storing evaluation at point x.
*/
private final Matrix tmpX;
/**
* Constructor.
*
* @param a Lower limit of integration.
* @param b Upper limit of integration.
* @param listener listener to evaluate a single dimension matrix function at required points.
* @throws WrongSizeException if size notified by provided listener is invalid.
*/
public MidPointMatrixQuadrature(
final double a, final double b, final MatrixSingleDimensionFunctionEvaluatorListener listener)
throws WrongSizeException {
this.n = 0;
this.a = a;
this.b = b;
rows = listener.getRows();
columns = listener.getColumns();
s = new Matrix(rows, columns);
tmpMid = new Matrix(rows, columns);
sum = new Matrix(rows, columns);
tmpX = new Matrix(rows, columns);
this.listener = listener;
}
/**
* Gets lower limit of integration.
*
* @return lower limit of integration.
*/
public double getA() {
return a;
}
/**
* Gets upper limit of integration.
*
* @return upper limit of integration.
*/
public double getB() {
return b;
}
/**
* Gets current value of integral.
*
* @return current value of integral.
*/
public Matrix getS() {
return s;
}
/**
* Returns the value of the integral at the nth stage of refinement.
*
* @param result instance where the value of the integral at the nth stage of refinement will
* be stored.
* @throws EvaluationException Raised if something failed during the evaluation.
*/
@SuppressWarnings("Duplicates")
@Override
public void next(Matrix result) throws EvaluationException {
try {
int it;
int j;
double x;
double tnm;
double del;
double ddel;
n++;
if (n == 1) {
// (s = (b - a) * func(0.5 * (a + b)))
func(0.5 * (a + b), tmpMid);
s.copyFrom(tmpMid);
s.multiplyByScalar(b - a);
result.copyFrom(s);
} else {
for (it = 1, j = 1; j < n - 1; j++) {
it *= 3;
}
tnm = it;
del = (b - a) / (3.0 * tnm);
// The added points alternate in spacing between del and ddel
ddel = del + del;
x = a + 0.5 * del;
sum.initialize(0.0);
for (j = 0; j < it; j++) {
func(x, tmpX);
sum.add(tmpX);
x += ddel;
func(x, tmpX);
sum.add(tmpX);
x += del;
}
// The new sum is combined with the old integral to give a refined integral
// s = (s + (b - a) * sum / tnm) / 3.0
sum.multiplyByScalar((b - a) / tnm);
s.add(sum);
s.multiplyByScalar(1.0 / 3.0);
result.copyFrom(s);
}
} catch (final WrongSizeException ex) {
throw new EvaluationException(ex);
}
}
/**
* Gets type of quadrature.
*
* @return type of quadrature.
*/
@Override
public QuadratureType getType() {
return QuadratureType.MID_POINT;
}
/**
* Gets number of rows of quadrature result.
*
* @return number of rows of quadrature result.
*/
@Override
protected int getRows() {
return rows;
}
/**
* Gets number of columns of quadrature result.
*
* @return number of columns of quadrature result.
*/
@Override
protected int getColumns() {
return columns;
}
/**
* Evaluates matrix function at x.
*
* @param x point where function is evaluated.
* @param result instance where result of evaluation is stored.
* @throws EvaluationException if evaluation fails.
*/
protected void func(final double x, final Matrix result) throws EvaluationException {
listener.evaluate(x, result);
}
}