In this example, the order of the Polyharmonic Spline, k =2. The function is sampled at 200 random points and the error is computed at 10000 random points.
using System; using Imsl.Math; public class RadialBasisEx2 { // The function to approximate static double fcn(double[] x) { return Math.Exp((x[1]) / 2.0) * Math.Sin(x[0]) - Math.Cos((x[1]) / 2.0); } public class PolyHarmonicSpline : RadialBasis.IFunction { virtual public int Order { get { return order; } } virtual public bool EvenOrder { get { return isEven; } } private int order = 3; private bool isEven = false; public PolyHarmonicSpline(int order) { this.isEven = order % 2 == 0; this.order = order; } public virtual double F(double x) { if (this.isEven) { return Math.Pow(x, order) * Math.Log(x); } return Math.Pow(x, order); } public virtual double G(double x) { if (order == 1) { return 1; } if (this.isEven) { return order * Math.Pow(x, order - 1) * Math.Log(x) + Math.Pow(x, order - 1); } return order * Math.Pow(x, order - 1); } } public static void Main(String[] args) { int nDim = 2; // Sample, with noise, the function at 100 randomly chosen points int nData = 200; double[,] xData = new double[nData,nDim]; double[] fData = new double[nData]; double[] row = new double[nDim]; Imsl.Stat.Random rand = new Imsl.Stat.Random(123457); rand.Multiplier = 16807; double[] noise = new double[nData * nDim]; for (int k = 0; k < nData; k++) { for (int i = 0; i < nDim; i++) { noise[k * 2 + i] = 1.0d - 2.0d * rand.NextDouble(); xData[k,i] = 3 * noise[k * 2 + i]; } // noisy sample for(int j = 0; j<nDim; j++) row[j]=xData[k,j]; fData[k] = fcn(row) + noise[k * 2] / 10; } // Compute the radial basis approximation using 100 centers int nCenters = 100; RadialBasis rb = new RadialBasis(nDim, nCenters); rb.RadialFunction = new PolyHarmonicSpline(2); rb.Update(xData, fData); // Compute the error at a randomly selected set of points int nTest = 10000; double maxError = 0.0; double aveError = 0.0; double maxMagnitude = 0.0; double[][] x = new double[nTest][]; for (int i2 = 0; i2 < nTest; i2++) { x[i2] = new double[nDim]; } noise = new double[nTest * nDim]; for (int i = 0; i < nTest; i++) { for (int j = 0; j < nDim; j++) { noise[i * 2 + j] = 1.0d - 2.0d * (double) rand.NextDouble(); x[i][j] = 3 * noise[i * 2 + j]; } double error = Math.Abs(fcn(x[i]) - rb.Eval(x[i])); maxMagnitude = Math.Max(Math.Abs(fcn(x[i])), maxMagnitude); aveError += error; maxError = Math.Max(error, maxError); } aveError /= nTest; Console.WriteLine("Average normalized error is " + aveError / maxMagnitude); Console.WriteLine("Maximum normalized error is " + maxError / maxMagnitude); Console.WriteLine("Using even order equation: " + ((PolyHarmonicSpline) rb.RadialFunction).EvenOrder); } }
Average normalized error is 0.0180558727060151 Maximum normalized error is 0.257565310407464 Using even order equation: TrueLink to C# source.