org.fxyz3d.geometry.GaussianQuadrature Maven / Gradle / Ivy
The newest version!
/**
* GaussianQuadrature.java
*
* Copyright (c) 2013-2016, F(X)yz
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
* * Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.
* * Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
* * Neither the name of F(X)yz, any associated website, nor the
* names of its contributors may be used to endorse or promote products
* derived from this software without specific prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
* WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
* DISCLAIMED. IN NO EVENT SHALL F(X)yz BE LIABLE FOR ANY
* DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
* (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
* LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
* ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
* SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
package org.fxyz3d.geometry;
import java.util.ArrayList;
import java.util.List;
import java.util.function.Function;
/**
*
* @author jpereda
*/
public class GaussianQuadrature {
private final double a;
private final double b;
private final List gauss=new ArrayList<>();
public GaussianQuadrature(){
this(5,-1d,1d);
}
public GaussianQuadrature(int points){
this(points,-1d,1d);
}
public GaussianQuadrature(int points, double a, double b){
this.a=a;
this.b=b;
switch(points){
case 2: gauss.add(new Pair(-0.5773502691896258, 1d));
gauss.add(new Pair(0.5773502691896258, 1d));
break;
case 3: gauss.add(new Pair(-0.7745966692414834, 0.5555555555555554));
gauss.add(new Pair(0d, 0.8888888888888888));
gauss.add(new Pair(0.7745966692414834, 0.5555555555555554));
break;
case 4: gauss.add(new Pair(-0.8611363115940526, 0.34785484513745396));
gauss.add(new Pair(-0.33998104358485626, 0.6521451548625462));
gauss.add(new Pair(0.33998104358485626, 0.6521451548625462));
gauss.add(new Pair(0.8611363115940526, 0.34785484513745396));
break;
case 5: gauss.add(new Pair(-0.906179845938664, 0.23692688505618922));
gauss.add(new Pair(-0.538469310105683, 0.4786286704993666));
gauss.add(new Pair(0, 0.5688888888888889));
gauss.add(new Pair(0.538469310105683, 0.4786286704993666));
gauss.add(new Pair(0.906179845938664, 0.23692688505618922));
break;
case 6: gauss.add(new Pair(-0.9324695142031519, 0.17132449237917097));
gauss.add(new Pair(-0.6612093864662645, 0.3607615730481388));
gauss.add(new Pair(-0.23861918608319727, 0.4679139345726911));
gauss.add(new Pair(0.23861918608319727, 0.4679139345726911));
gauss.add(new Pair(0.6612093864662645, 0.3607615730481388));
gauss.add(new Pair(0.9324695142031519, 0.17132449237917097));
break;
default: gauss.add(new Pair(-0.9491079123427586, 0.12948496616886915));
gauss.add(new Pair(-0.7415311855993945, 0.27970539148927703));
gauss.add(new Pair(-0.40584515137739713, 0.38183005050511903));
gauss.add(new Pair(0d, 0.4179591836734694));
gauss.add(new Pair(0.40584515137739713, 0.38183005050511903));
gauss.add(new Pair(0.7415311855993945, 0.27970539148927703));
gauss.add(new Pair(0.9491079123427586, 0.12948496616886915));
break;
}
}
public double NIntegrate(Function f){ // f[t]
return gauss.parallelStream()
.mapToDouble(p->p.getWeight()*((double)f.apply((b-a)/2d*p.getAbscissa()+(b+a)/2d)))
.reduce(Double::sum).getAsDouble()*(b-a)/2d;
}
// // Test
// public static void main(String[] args) {
// GaussianQuadrature gauss=new GaussianQuadrature(7,0,Math.PI);
// System.out.println(""+gauss.NIntegrate(x->Math.sin(x)));
// }
private class Pair {
private final Double abscissa;
private final Double weight;
public Pair(double abscissa, double weight){
this.abscissa=abscissa;
this.weight=weight;
}
public double getAbscissa() { return abscissa; }
public double getWeight() { return weight; }
}
}