|
||||||||
| PREV NEXT | FRAMES NO FRAMES | |||||||
See:
Description
| Packages | |
| ru.sscc.matrix | Collection of classes describing algebraic matrices. |
| ru.sscc.matrix.solve | Collection of Solvers intended for solving of Systems of Linear Algebraic Equations (SLAE). |
| ru.sscc.spline | Basic spline approximation classes. |
| ru.sscc.spline.analytic | Multi-dimensional Duchon's splines (not prepared yet). |
| ru.sscc.spline.base | The collection of auxiliary classes used in spline approximation. |
| ru.sscc.spline.polynomial | One-dimensional polynomial splines. |
| ru.sscc.util | Collection of miscellaneous utility classes. |
| ru.sscc.util.data | Collection of data access classes implementing such concepts as a data container, data vector, data pointer, and array of vectors. |
Version 1.22.
The first and main aim of the development of this Library was an implementation of various splines to be constructed by function measurements on 1D or multi-dimensional mesh of scattered points. This job was originated on the spline software library LIDA-3 developed in 1985 at the Novosibirsk Computing Center under the leadership of Prof. Vladimir A. Vasilenko. The current job began at the middle of 1999. While developing, the initial scope of the Library was enlarged. It now contains 3 main sections:
The scope of the Library will be continuously enlarged to cover other fields of the Numerical Mathematics.
The ru.sscc.util.data package provides unified access to user's data vectors. The enumeration of vector's entries always starts from zero. All package classes support four primitive operations upon data entries:
Four conceptual classes are introduced here:
RealContainerrelativeAccuracy()
method returns the relative accuracy of data type used in the
container, e.g. the smallest value e such that 1+e
isn't equal to 1 for the data type used. This constant is widely useful in
numeric calculations, but it is absent in Java.
RealVector
RealPointerSimplePointer,
and to point to entries of the RealVectors instance, the
RealPointers.
RealVectorssquaredDistance(index,pointer)
method is used.
All above mentioned classes are tied together: the RealContainer class allows to create a vector or pointer upon it; the RealVector class allows to create a subvector or pointer to a subvector in it; the RealVectors class provides getting of i-th vector from the set of vectors, creation of a vector consisting of j-th component of all vectors, creation of a pointer to i-th vector. All (sub)vector creation operations don't create a copy of entries. They create an instance that refers to the same underlying data container as the instance with the help of which the operation was performed.
A container can move or shift the starting index of a vector created by it. This operation is useful in matrix calculations when a matrix is considered as a set of rows or columns, and the same operation must be produced on every element of the set.
To simplify conversions, a number of static by methods is introduced.
For example, to create a RealContainer from a double[] or
float[] array, you can write down the following:
double[] array = new double[count];
//... Fill the array in
RealContainer container = RealContainer.by(array);
Finally, we notice the RealMath class. It
contains many useful static methods that perform some BLAS operations upon
double[], float[], and RealPointer instances. The
calculation of factorial of n, integer power of x, and a polynomial value
by the Gorner's algorithm is also implemented here.
The ru.sscc.matrix package and its subpackages present the Numerical Linear Algebra Section of the Library. The root package describes matrices. Only real matrices are implemented now.
The Matrix interface describes the most common properties
of any matrix: a matrix is a rectangular table consisting of a constant number
of rows and columns and a matrix is serializable object.
The RealAlgebraicMatrix interface extends
the Matrix by adding two algebraic methods: multiplication of the matrix
by a vector and multiplication of the transposed matrix by a vector. In other
words, the algebraic matrix is the operator that linearly maps a source vector
to the target vector and it has the transposed operator acting analogously. In
addition, the method relativeAccuracy() must return the relative
accuracy of data type used for matrix entries.
When a matrix is used for solving of a linear algebraic system, the decomposition coefficients are stored in it. In this case, the matrix loses the algebraicity and the multiplication by a vector has no sense more. We say that the matrix becomes non-algebraic and the try to multiply it by a vector throws the IllegalStateException.
The RealMatrix class is the basic class for all real
matrices. It is the abstract class. It partially implements the
RealAlgebraicMatrix interface and adds four primitive operations for the
access to matrix values,
get(i,j),
set(i,j,value),
add(i,j,value), and
mul(i,j,value).
The indices of matrix entries start from zero. So, if the matrix has m rows and n columns then i must belong to the set {0,1,...,m-1} and j must belong to the set {0,1,...,n-1}.
All real matrices store their coefficients in the body, an instance of
the
RealContainer type, e.g. the 1D array is used as
the storage for matrix entries. This solution was selected by the efficiency
reason. You can get the matrix body by the
getContainer()
method. The access to the matrix dimensions may be done in two way:
rowsNumber() and
columnsNumber() methods inherited
from the Matrix interface and
nRows and
nColumns.
The RealMatrix class supports cloning. The clone is the matrix having the separate body copied from the body of the original matrix. The clone is the algebraic matrix, independently to the state of the original matrix. If the original matrix is a submatrix of a dense matrix, the clone will contain the submatrix entries only.
Three types of matrices are implemented in the current version of the Library:
DenseMatrix
SymBandedMatrix
RectBandedMatrixThe switching between algebraic and non-algebraic state of matrices is implemented by the lock and reuse methods: the lock method changes the matrix state to be non-algebraic and the reuse method allows the reusing of the matrix as the algebraic matrix.
The package supports matrix-by-matrix multiplication operations as static
methods of the RealMatrix class.
The ru.sscc.matrix.solve package is devoted to the solving of a System of Linear Algebraic Equations (SLAE)
|
Every real solver implements the
RealSolver
interface, containing the
solve(b,x)
method that receives the right hand side vector b and calculates the solution
vector x. Two auxiliary methods are also described in the interface:
sourceSize() and
targetSize().
The RealSolver interfaces supports serialization.
We distinct two main types of real solvers: direct solvers and iterative
solvers. Only direct solvers are implemented yet. The basic abstract class for
all direct solvers is the
RealDirectSolver.
It implements the basic interface for all direct solvers and adds the abstract
method
factorize()
that provides the matrix decomposition by factors (e.g. factorization). The
matrix to be factorized is attached to the solver. A direct solver can have two
states: before the factorization and after it. If the matrix is nonfactorized
yet, the solve method is inaccessible. When the matrix is factorized, it
becomes non-algebraic, and all algebraic operations with it are blocked. The
reuse method implemented in all subclasses changes the state of the
solver and the attached matrix to the initial state: the factorization tag is
cleared and the matrix becomes algebraic again. So, after that, you can fill in
the attached matrix again.
The basic principal for all direct solvers is: "Once factorize, many times solve". So, if you need to solve many SLAEs that differ in the right hand sides only, you should factorize the matrix only once, and then can find the solution many times for different right hand sides.
To improve the solution accuracy, use the iterative refinement algorithm implemented in the RealDirectSolver. It is most efficient if the attached matrix is constructed on the base of the float[] type vector. The construction of the inverse matrix is also implemented in this class.
The hierarchy of the classes based on the RealDirectSolver is the following:
CholeskyBandSolver
that solves the SLAE with a positive definite symmetric banded matrix by
Cholesky Method also known as Square Root Method. The solver having the
specialized algorithm, the
GreenSolver,
will be added soon. It will be used in spline approximation algorithm based on
analytic representation of spline.
RealSquareSolver
is the basic class for solvers oriented to the systems with dense square
matrices. We implemented the following subclasses: (1) the
CholeskySolver
solving the system with a positive definite symmetric matrix; (2) the
GaussSolver
solving the system with a positive or negative definite non-symmetric matrix
(in this case the Gauss eliminations may be produced directly and pivoting
isn't needed); and (3) the
CrautSolver
solving the system with a common matrix by the Craut method with partial
pivoting (the Craut method is the modification of Gauss eliminations).
RealCommonSolver
is the basic class for solvers oriented to the systems with common dense
rectangular matrices. The system with square matrix may be solved by a solver
of this class also. This class supports the balancing of matrices that may
improve the accuracy of solution. Three subclasses are implemented now: (1) the
RotationSolver
based on Givens rotations; (2) the
ReflectionSolver
based on Hausholder reflections (it 30% more efficient for speed, but loses in
accuracy; to provide the same accuracy as the previous solver, the calculation
of inner products with quadruple precision is needed, but this is impossible in
pure Java); and (3) the
EliminationSolver
based on Gauss eliminations with partial pivoting (it is twice faster that the
previous solver, but its accuracy is worst of all).
The ru.sscc.spline package contains basic spline approximation classes. We use the term spline for a function
|
|
A spline instance consists of 3 parts:
SplineBody interface.
SplineWorkspace is basic
class for a spline workspace.
The BaseSpline class is the basic class for all
splines. It implements common operations for all splines: the
dimension() operation returns the
dimension m of the space of points Rm for which the spline was
constructed; the
clone() operation creates a clone of the spline
(the clone shares the body and coefficients with the parent instance, but has
its own workspace); and the
flush() operation removes the workspace from
a spline. You can also use the
getBody() and
getContainer() operations for the
access to the spline body and coefficients.
All splines are serializable objects.
The BaseSpline instance cannot be created directly. But it has two subclasses having public constructors: the Spline and Splines. The first class performs operations with a unique spline, and the last class performs operations with many splines constructed on the same spline body.
The Spline class has a number of value()
operations. For 1D spline, you can use any of them, but, for multi-D spline,
you must use only the multi-D variant of the operation. If 1D spline allows to
calculate a spline derivative, the
value(point,derivative) operation
may be used. For the multi-D spline, you can calculate a spline value only. If
the calculation cannot be done for the parameters specified, the
IllegalArgumentException is thrown.
The Splines class has more value() operations: as
individual spline value calculation so as vector operations for all splines of
the instance.
The ru.sscc.spline.polynomial package is intended for the construction of 1D piecewise polynomial splines. The polynomial splines differ by the spline degree and spline representation. The construction of polynomial splines of odd degree with polynomial representation in mesh nodes is available now. We plan to add the construction of even degree splines and also implement the B-spline representation for all splines in future releases.
The POddSplineCreator class allows
to construct polynomial splines of odd degree. To describe its features, we
introduce the notation on the example of cubic splines.
Let x0 < x1 < ¼ < xn-1 be mesh nodes on the real axis and fi = f(xi) be the values of the function to be interpolated. The well-known cubic spline is the function constructed with the cubic polynomials defined at mesh cells (xi,xi+1), i = 0,¼,n-2, in such a way to provide its second derivative to be continuous function at mesh nodes. The interpolating cubic spline satisfies the interpolation conditions
|
|
To construct and use the natural interpolating cubic spline, write the
following code:
// Create a mesh array. It may be a double[], float[], or RealContainer instance.
// For example:
double[] mesh = new double[n];
// ... fill the mesh
// Create an array of interpolation values at mesh nodes. For example:
double[] values = new double[n];
// ... fill the array
// Create the natural interpolating cubic spline
Spline spl = POddSplineCreator.createSpline(2, mesh, values);
// Use the spline in calculations: x - a point on the real axis, f - the spline value
f = spl.value(x);
The first parameter of all createSpline methods in the
POddSplineCreator class
is the difference order that is tied with the spline degree by the
following rule: degree=2*order-1. So, to create the cubic spline, the
order must be equal to 2. If order=1, the piecewise linear spline will be
constructed. For order=3,4,..., the quintic spline, 7th degree spline, etc. will be constructed. The more order is used, the higher spline derivative is
continuous. Exactly, the (2*order-2)-th spline derivative is the
continuous function. The extrapolation (e.g. calculation of the spline
values for x < x0 and x > xn-1) is provided by the polynomials of
(order-1)-th degree. So, the cubic spline is extended by linear
functions, the quintic spline is extended by quadratic functions and so on.
The number of mesh nodes must satisfy the inequality n>=order. The mesh
nodes must be ordered in increasing or decreasing order and must be different.
If the mesh is incorrect, or some mesh nodes are very closed, or the difference
order is too large, the
CalculatingException will be thrown while
the spline construction. This means that the spline cannot be constructed. You
must catch this exception anywhere.
If the distance between mesh nodes is constant, you need not fill in the mesh
array. Use the following method instead:
Spline spl = POddSplineCreator.createSpline(order, x0, step, n, values);
Here x0 is the first mesh node, step is the distance between
neighbour mesh nodes, n is the number of mesh nodes, and values
is the array of interpolating values at mesh nodes.
When the function measurements are non-strict, the constructing of smoothing spline is more preferable than the interpolation. The smoothing spline goes "near" the function measurements. It is more smooth and may suppress unexpected overshoots in the measurements. The degree of smoothing is controlled by the smoothing parameter a > 0. When a tends to zero, the smoothing spline converges to the interpolation spline. When a tends to the positive infinity, the smoothing spline converges to the polynomial of the degree order-1 (this is the linear function for cubic spline). The problem is to select proper a providing the "best" smoothing results.
We select the smoothing parameter to satisfy the residual criterion
|
It isn't necessary to solve this equation with high accuracy. It is usually
enough to construct the spline, such that the residual belongs to the segment
[e(1-d),e(1+delta)]. Here d is the accuracy
level. The default accuracy level is set to 0.1. It may be changed by the
SmoothingPreparator.setAccuracy(delta)
method.
To construct a smoothing spline satisfying the residual criterion, the
following methods are used:
// For arbitrary mesh:
Spline spl = POddSplineCreator.createSpline(order, mesh, values, epsilon);
// For uniform mesh:
Spline spl = POddSplineCreator.createSpline(order, x0, step, n, values, epsilon);
If epsilon isn't less than emax, the spline will
coincide with the polynomial of the (order-1)-th degree-the solution
to the smoothing problem for a = ¥.
If you know the function measurement error for every mesh node, you'll want to construct the smoothing spline that may strongly deviate from the function measurements in "bad" nodes, and weekly deviate in the "good" nodes. To provide this, use the weighted residual criterion
|
To construct a smoothing spline satisfying the weighted residual criterion, the
following methods are used:
// For arbitrary mesh:
Spline spl = POddSplineCreator.createSpline(order, mesh, values, epsilon, weights);
// For uniform mesh:
Spline spl = POddSplineCreator.createSpline(order, x0, step, n, values, epsilon, weights);
Here the weights array contains positive weights assigned to mesh nodes.
If you want to construct many splines distinct in function measurements only, you must do the following:
POddSplineCreator class.
The Step 3 may be repeated many times.
Example 1 To construct a number of interpolating splines on an
arbitrary mesh, write down the following:
double[] mesh = new double[n];
// ... fill the mesh
POddSplineCreator creator = new POddSplineCreator(order, RealVector.by(mesh));
creator.prepareSolver(0); // Prepare solver for alpha=0
// Create array for the measurements
double[] values = new double[n];
// ... prepare function measurements in the array
// Construct the first spline:
Spline spl1 = creator.constructSpline(RealVector.by(values));
// ... prepare function measurements in the same array
// Construct the second spline:
Spline spl2 = creator.constructSpline(RealVector.by(values));
// And so on ...
Example 2 To construct the weighted smoothing splines with the same
smoothing parameter a (1 in this example), write down the following:
POddSplineCreator creator = new POddSplineCreator(order, RealVector.by(mesh));
double[] weights = new double[n];
// ... fill the weight vector
creator.prepareWeights(weights);
creator.prepareSolver(1); // Prepare solver for alpha=1
// ... construct splines by the same manner as in the previous example
Example 3 To construct the smoothing spline with selection of
a from the residual criterion, you need not prepare the solver. Instead
of the preparation of the solver you can immediately use the constructing
method:
Spline spl1 = creator.constructSpline(RealVector.by(values), epsilon);
After that, the solver will be prepared. So, to construct other splines with
the same parameter a, you can call:
Spline spl2 = creator.constructSpline(RealVector.by(values));
// ...
Example 4 When many splines must be constructed in one Splines
instance, you should prepare the RealVectors instance containing all
functions measurements and use one of the following methods:
RealVectors measurements;
// ... prepare the measurements
// Construct splines by the rows of the measurements array:
Splines spls = creator.constructSplinesByRows(measurements);
// or construct splines by the columns of the measurements array:
Splines spls = creator.constructSplinesByColumns(measurements);
Example 5 If you want to create the Splines instance with
sequential calculation the splines' coefficients, do the following. Manually
construct the Splines instance:
Splines spls = new Splines(creator.getBody(), count);
Here count is a number of splines to be calculated. After that, call the
calculate method in a cycle:
for(int i=0; i<count; i++) {
// ... prepare i-th measurements in the values array ...
// Calculate the coefficients of i-th spline
creator.calculate(RealVector.by(values), spls.getVector(i));
}
|
||||||||
| PREV NEXT | FRAMES NO FRAMES | |||||||