Sunday, March 8, 2015

Function Approximations in Java

FunctionApproximations.java

A class containing methods to calculate the coefficients and values based on a set of initial values.

Supports:


Example: Approximating sine values based on 12 initial points.


int n = 11;

double[] stx = new double[n + 1];

double[] y = new double[n + 1];

double dif = (Math.PI * 2) / (n);

double start = -Math.PI;

for (int i = 0; i < n + 1; i++) {

    stx[i] = start;

    start += dif;

    y[i] = (double) Math.round(Math.sin(stx[i]) * 100000) / 100000;

}

double[] newtonPolCoef = FunctionApproximations.newtonPolCoef(stx, y);

double[] splineCoef = FunctionApproximations.splineCoef(stx, y, 3);

double[] lsCoef = FunctionApproximations.leastSqCoeff(stx, y, 5);

Random random=new Random();

double val,correct,newton,spline,leastsq;

for (int i = 0; i < 10; i++) {

    val=random.nextDouble();

    correct=Math.sin(val);

    newton=FunctionApproximations.newtonPolVal(newtonPolCoef, stx, getX(val));

    spline=FunctionApproximations.splineVal(splineCoef, stx, getX(val));

    leastsq=FunctionApproximations.leastSqVal(lsCoef, getX(val));

    System.out.println("sin("+val+")");

    System.out.println("Newton approximation :     "+newton+"    off by   "+(correct-newton));

    System.out.println("Cubic splines :            "+spline+"    off by   "+(correct-spline));

    System.out.println("Least squares 5th degree : " +leastsq+"    off by   "+(correct-leastsq));

    System.out.println("----");

}

static private double getX(double x) {

    if (Math.abs(x) <= Math.abs(Math.PI)) {

        return x;

    } else {

        if (x > 0) {

            x -= Math.round(x / (2 * Math.PI)) * (2 * Math.PI);


        } else {

            x += Math.round(Math.abs(x) / (2 * Math.PI)) * (2 * Math.PI);

        }

        return x;

    }

}

Sample output:

sin(0.8481863801681208)
Newton approximation :     0.7500825572784057    off by   -3.455657819895208E-7
Cubic splines :            0.7500853256850928    off by   -3.113972469104276E-6
Least squares 5th degree : 0.7446845168230829    off by   0.005397694889540783
----
sin(0.17370422573181987)
Newton approximation :     0.1728302183995918    off by   1.7900743105314643E-6
Cubic splines :            0.17285859735632642    off by   -2.6588882424072313E-5
Least squares 5th degree : 0.1704451410039368    off by   0.0023868674699655534
----
sin(0.9222664997998143)
Newton approximation :     0.7969736441928723    off by   -9.775787263022195E-7
Cubic splines :            0.7969127151051372    off by   5.995150900883761E-5
Least squares 5th degree : 0.7921460535443667    off by   0.004826613069779273
----
sin(0.42918641383160216)
Newton approximation :     0.4161283242433726    off by   2.818337875132304E-6
Cubic splines :            0.41601147561349894    off by   1.1966696774878827E-4
Least squares 5th degree : 0.4110357193164086    off by   0.005095423264839138
----
sin(0.5240316104273179)
Newton approximation :     0.5003722902694624    off by   2.5088399099315595E-6
Cubic splines :            0.5002095950629291    off by   1.6520404644326803E-4
Least squares 5th degree : 0.4946962081563806    off by   0.0056785909529917244
----
sin(0.9800507235185272)
Newton approximation :     0.8305269746540256    off by   -1.3510869794064462E-6
Cubic splines :            0.8303843741206196    off by   1.4124944642657233E-4
Least squares 5th degree : 0.8262637304226693    off by   0.00426189314437686
----
sin(0.9303877316203922)
Newton approximation :     0.8018527171875109    off by   -1.0374269545643244E-6
Cubic splines :            0.8017809837980641    off by   7.069596249231758E-5
Least squares 5th degree : 0.7970983401665966    off by   0.004753339593959738
----
sin(0.6236951478334793)
Newton approximation :     0.5840367621157957    off by   1.826036797991648E-6
Cubic splines :            0.5838957413617496    off by   1.428467908440867E-4
Least squares 5th degree : 0.5780582273554796    off by   0.005980360797114059
----
sin(0.5844538429625159)
Newton approximation :     0.5517418299721984    off by   2.1323006360596253E-6
Cubic splines :            0.5515834663608791    off by   1.6049591195543833E-4
Least squares 5th degree : 0.5458425471133572    off by   0.005901415159477286
----
sin(0.37727308863399756)
Newton approximation :     0.36838388119885634    off by   2.8261127495432525E-6
Cubic splines :            0.36831036229846276    off by   7.634501314313091E-5
Least squares 5th degree : 0.36372091601169265    off by   0.004665791299913236
----

Least squares method in Java

Wikipedia

Calculating the coefficients:

static public double[] leastSqCoeff(double[] x, double[] y, int degree) {

    if (degree < 0) {
        throw new IllegalArgumentException();
    }
    ++degree;
    int n = x.length;
    double[][] a = new double[n][degree];
    double[] b = new double[degree];
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < degree; j++) {
            a[i][j] = Math.pow(x[i], j);
        }
    }
    double[][] ata = new double[degree][degree];
    for (int i = 0; i < degree; i++) {
        for (int j = 0; j < degree; j++) {
            ata[i][j] = MatricesBasic.vectorM(MatricesBasic.getColumn(a, i), MatricesBasic.getColumn(a, j));
        }
        b[i] = MatricesBasic.vectorM(MatricesBasic.getColumn(a, i), y);
    }
    b = Matrices.gauss(ata, b);
    return b;

}

Approximating the value based on the coefficients:

static public double leastSqVal(double[] coef, double val) {

    double res = 0;
    for (int i = 0; i < coef.length; i++) {
        res += coef[i] * Math.pow(val, i);
    }
    return res;
}

This is part of the Function Approximations project.

Creating order 2 and 3 splines in Java

Wikipedia

Calculating the coefficients:

static public double[] splineCoef(double[] x, double[] y, int degree) {

    if (degree != 2 && degree != 3) {
        throw new IllegalArgumentException();
    }
    int n = (x.length - 1) * (degree + 1);
    double[][] a = new double[n][n];
    double[] b = new double[n];
    double[] coef;
//s
    int k = 0;
    for (int i = 0; i < (x.length - 1) * 2; i += 2, k++) {
        for (int j = 0; j < degree + 1; j++) {
            for (int l = 0; l < 2; l++) {
                a[i + l][j + ((i / 2) * (degree + 1))] = Math.pow(x[k + l], degree - j);
            }
        }
        for (int l = 0; l < 2; l++) {
            b[i + l] = y[k + l];
        }
    }
//s'
    int line = (x.length - 1) * 2;
    int column = 0;
    for (int i = 1; i < x.length - 1; i++) {
        for (int j = 0; j < degree; j++) {
            a[line + i - 1][column + j] = (degree - j) * Math.pow(x[i], degree - 1 - j);
            a[line + i - 1][(column + j) + degree + 1] = -a[line + i - 1][column + j];
        }
        column += degree + 1;
    }
//s"
    if (degree > 2) {
        line += x.length - 2;
        column = 0;
        for (int i = 1; i < x.length - 1; i++) {
            for (int j = 0; j < degree - 1; j++) {
                a[line + i - 1][column + j] = ((degree - j) * (degree - 1 - j)) * Math.pow(x[i], degree - 2 - j);
                a[line + i - 1][(column + j) + degree + 1] = -a[line + i - 1][column + j];
            }
            column += degree + 1;
        }
        if (degree == 3) {
            a[n - 2][0] = 6 * x[0];
            a[n - 2][1] = 2;
            a[n - 1][n - 4] = 6 * x[x.length - 1];
            a[n - 1][n - 3] = 2;
        }

    } else {
        a[n - 1][0] = 1;
    }
    Jama.Matrix A = new Jama.Matrix(a);
    Jama.Matrix c = new Jama.Matrix(b, n);
    Jama.Matrix sol = A.solve(c);
    coef = sol.getColumnPackedCopy();
    return coef;
}

Approximating the value based on the coefficients:

static public double splineVal(double[] coef, double[] x, double val) {

    if (val < x[0]) {
        return Double.NaN;
    }
    int degree = coef.length / x.length;
    for (int i = 1; i < x.length; i++) {
        if (val <= x[i]) {
            int indx = (i - 1) * (degree + 1);
            double s = 0;
            for (int j = 0; j < degree + 1; j++) {
                s += coef[indx + j] * Math.pow(val, degree - j);
            }
            return s;
        }
    }
    return Double.NaN;
}

This is part of the Function Approximations project.

Newton polynomial calculation in Java

Wikipedia

Calculating the coefficients:

static public double[] newtonPolCoef(double[] x, double[] y) {

    int n = x.length - 1;
    double dd[][] = new double[n][];
    double[] a = new double[n + 1];
    for (int i = 0; i < n; i++) {
        dd[i] = new double[n - i];
    }
    for (int j = 0; j < n; j++) {
        dd[0][j] = (y[j + 1] - y[j]) / (x[j + 1] - x[j]);
    }
    for (int i = 1; i < n; i++) {
        for (int j = 0; j < n - i; j++) {
            dd[i][j] = (dd[i - 1][j + 1] - dd[i - 1][j]) / (x[j + i + 1] - x[j]);
        }
    }
    a[0] = y[0];
    for (int i = 1; i < n + 1; i++) {
        a[i] = dd[i - 1][0];
    }
    return a;
}

Approximating the value based on the coefficients:

static public double newtonPolVal(double[] coef, double[] x, double val) {

    double r = coef[0];
    double mul = 1;
    for (int i = 1; i < x.length; i++) {
        mul *= (val - x[i - 1]);
        r += coef[i] * mul;
    }
    return r;
}

This is part of the Function Approximations project.