The nsplines package contains a set of *.m files for use with GNU Octave and a *.sci file for use with SciLab to perform natural spline interpolation for curves specified as a set of discrete points.
Polynomial fitting attempts to create one polynomial to represent the
entire curve. The summary of the quadratic distances between the polynomial
and the points is minimized. Depending on the polynomial's degree the
distances between polynomial and given points differs, for higher degree
polynomials you typically have to expect oscillations.
Natural spline interpolation constructs a sequence of cubic splines, hitting all the specified points. The distance between the interpolated curve and the given points is 0. In addition to "the curve must hit each point" the following constraints are used to calculate the polynomial coefficients:
For modified natural splines we are allowed to specify the value of the first derivative in some or all points. In this case the constraints are changed, the above constraints are replaced by
For a comparison here are two images showing a natural spline interpolation for some measurement points (first picture) and several polynomial fittings for the same points (second image). None of the polynomials exactly hits all the specified points, none of the polynomials matches a common sense free-hand drawn curve for the points.
See the install.html page.
The curve points are specified as a matrix. The matrix has two columns:
the first column is for x-values, the second column is for the y-value.
Each point is specified in one matrix line.
You have to use a 4-column matrix. Each matrix line describes one point.
The x-value is placed in the first column, the y-value is placed in the
If the first derivative for a point is given, the third column must have a value > 0 assigned, the fourth column contains the first derivative of the function in the point.
Otherwise the third column must have a value <=0, the contents of the fourth column doesn't matter.
Invoke the function ns_create() on your input data matrix. This function returns a matrix containing interval borders and polynomial coefficients.
In GNU Octave all variables are local by default (I don't know exactly
about SciLab). If a function wants to access a global variable, the
variable must be declared to be global twice: once on the file level before
a value is assigned and again in the function to access the variable.
The purpose of nsplines is to have functions available we can use for calculation and iteration.
This is an internal function used by ns_value(). You should not need to deal directly with this function.
Once you have run ns_create() and saved the result to variable M you can run ns_value(M, x) to calculate the function result for a given x.
The ns_gpfile() function only writes the data points
and the polynomial coefficients and functions to the GnuPlot file.
I recommend to customize this file before running GnuPlot on it, you should add labels for x- and y-axis, specify output file type and file name...
An attempt to write a complete GnuPlot file would require to have a pendant for each possible GnuPlot setting in the nsplines package. Users would have to learn the GnuPlot syntax, the GNU Octave/SciLab command for the setting and the nsplines command for the same setting...
I avoided this overhead, to obtain good results from GnuPlot users should only need to know GnuPlot commands.
See the example page.
Sometimes this problem occurs if x- and y-value ranges differ too much
(several powers of 10).
You can attempt to do the following:
For an example, change the GNU Octave file
mw = [1.8850e+03 5.2600e+01 3.7700e+03 5.2600e+01 5.6550e+03 5.2600e+01 8.1680e+03 5.2600e+01 9.4250e+03 5.2600e+01 1.2566e+04 5.2700e+01 1.5708e+04 5.2800e+01 1.8850e+04 5.3000e+01 3.7699e+04 5.3200e+01 9.4248e+04 5.4100e+01 1.5708e+05 5.6400e+01 2.1991e+05 5.8900e+01 2.8274e+05 6.3000e+01 3.4558e+05 6.7000e+01 4.0841e+05 7.2000e+01 4.7124e+05 7.6000e+01 5.3407e+05 8.1000e+01 5.9690e+05 8.5000e+01 6.5973e+05 9.0000e+01 7.2257e+05 9.4500e+01]; global M; M = nssp_create(mw); nssp_gpfile(mw,M,'test.gp');
mw = [1.8850e+03 5.2600e+01 3.7700e+03 5.2600e+01 5.6550e+03 5.2600e+01 8.1680e+03 5.2600e+01 9.4250e+03 5.2600e+01 1.2566e+04 5.2700e+01 1.5708e+04 5.2800e+01 1.8850e+04 5.3000e+01 3.7699e+04 5.3200e+01 9.4248e+04 5.4100e+01 1.5708e+05 5.6400e+01 2.1991e+05 5.8900e+01 2.8274e+05 6.3000e+01 3.4558e+05 6.7000e+01 4.0841e+05 7.2000e+01 4.7124e+05 7.6000e+01 5.3407e+05 8.1000e+01 5.9690e+05 8.5000e+01 6.5973e+05 9.0000e+01 7.2257e+05 9.4500e+01]; factor = 1000.0; [r,c]=size(mw); for i = 1:r mw(i,1) = mw(i,1) / factor; endfor global M; M = nssp_create(mw); for i = 1:r mw(i,1) = mw(i,1) * factor; endfor [r,c]=size(M); for i = 1:r M(i,1) = M(i,1) * factor; M(i,2) = M(i,2) * factor; M(i,3) = M(i,3) / (factor * factor * factor); M(i,4) = M(i,4) / (factor * factor); M(i,5) = M(i,5) / factor; endfor nssp_gpfile(mw,M,'test.gp');