public interface ScalarFunction { float fun(float arg); } public class Interval { Interval (float fa, float fb) { a = fa; b = fb; } Interval devideLeft () { Interval i = new Interval(a,(a+b)*0.5); return i; } Interval devideRight () { Interval i = new Interval((a+b)*0.5,b ); return i; } float a; float b; } public class Knot implements Cloneable { Knot ( float ft, float fval) { t=ft; val=fval; } float t; float val; Object clone() { return new Knot(t,val); } } public class KnotCompare implements Comparator { boolean equals (Object a, Object b) { return (abs( (((Knot)a).val- ((Knot)b).val)/max(((Knot)a).val,((Knot)b).val)) < 1e-6); } int compare (Object a, Object b) { float f = ((Knot)a).val - ((Knot)b).val; if ( f <= 0) return -1; else return 1; } } public class GaussQuad { float adaptiveIntegrate (ScalarFunction f, float a, float b, float disp ) { intervals.clear (); results.clear(); intervals.add ( new Interval(a,b) ); results.add ( new Knot(a,disp) ); // subdevide at least 32 times /*int iters = 5; for ( int i =0; i < iters; ++i) { ArrayList newIntervals = new ArrayList(); for ( int j=0; j < intervals.size(); ++ j) { Interval inter = (Interval) intervals.get(j); newIntervals.add (inter.devideLeft()); newIntervals.add (inter.devideRight()); } intervals = newIntervals; } */ while ( intervals.size() != 0) { Interval first = (Interval) intervals.remove(0); float left = integrateGauss (f, first.devideLeft() ); float right = integrateGauss (f, first.devideRight()); float both = integrateGauss (f,first); // prevent algorithm for endless loop, so check for too small intervals if ( (abs(left+right-both) < 1e-6 || (first.b-first.a) < 1e-5) ) { Knot lastKnot = (Knot) results.get ( results.size()-1); results.add ( new Knot(first.b, left+right+lastKnot.val ) ); } else { intervals.add(0, first.devideRight() ); intervals.add(0, first.devideLeft()); } } return ((Knot) results.get(results.size()-1)).val; } ArrayList getKnotVals() { ArrayList ret = new ArrayList (); for (int i=0, n = results.size(); i < n; ++i ) ret.add( ((Knot) results.get(i)).clone() ); return ret; } float subsParam (float t, Interval i ) { return ((i.b - i.a)*t + i.a + i.b)*0.5; } float integrateGauss (ScalarFunction f, Interval inter) { float sum =0; for (int i =0; i < 5; ++i) sum += f.fun( subsParam(knots[i], inter) ) * weights[i]; return (inter.b-inter.a)*0.5*sum; } float[] knots = { -0.906179845938664, -0.538469310105683, 0.0, 0.538469310105683, 0.906179845938664 }; float[] weights = { 0.236926885056189, 0.478628670499366, 0.568888888888889, 0.47862867049936, 0.236926885056189 }; ArrayList intervals = new ArrayList(); ArrayList results = new ArrayList(); }