OPEM2 - orthonormal polynomial expansion method in two dimensions and its realization in fortran package

Author: Nina B.Bogdanova
Institute of Nuclear Research and Nuclear Energy, Sofia, Bulgaria
You are
visitor here.

Language: Fortran 77, Windows

A package of Fortran 77 programs is presented. This package has developed since 1981 for the calibration problems of measurement devices in High Energy Physics in collaboration with LIT JINR, Dubna.

Last version of the package is developed in part under the contract of Bulgarian National Council for Sientific Research - Phi 1001.

1.Description of the Method

Orthonormal Polynomial Expansion Method (OPEM2) is the method of orthogonalization of polynomials based on Householder-Forsythe 3-term relation [1,2] (one-dimensional case) and the Weisfeld recurrence [3] (many-dimensional case). The recurrence family of polynomials {PL} is given by:

f1 (1.1)

where Pk is a mother term. K and I have to be determined. {PL} are defined in a finite discrete subset D of rn, D={q1,q2,...,qM}. Every point qj=qj(x1j1,x2j2,...,xnjn) presents a n-coordinate vector (in this case n=2). The weights, associated with them, wj=1/S2j, depend on the standard deviation in every point. The points, the weights and the experimental function values {qj,wj,fjexp}j=1M are given. The recurrence coefficients {alfaL} and {betaL-1} and normalizing factor {cL} are evaluated by scalar products of the given values. The description here follows the descriptions [3,4,7]. In OPEM2 the approximated function is given by the expansion of the orthonormal polynomials using special orthonormal coefficients ak. Optimal value of polynomial degree - Lr is evaluated using two criteria.

f2 (1.2)

Then orthonormal coefficients {ak} in (1.2) are easily computed by given values:

f3 (1.3)

A pyramid of two-dimensional polynomials is given - the hierarchy order of polynomials in two variables:


In the cell on the top the number of polynomial is given, and the degrees j1 and j2 of the x1 and x2 variables also - on the bottom. One example for evaluating K and I in (1.1) gives for the polynomial P13 (k=2, j1=2, j2=2):

f4 (1.4)

(P8 is a mother term: K=8, I=4).

New theoretical investigations are presented in two-dimensional case in [7] for the case of differentiation and integration - the fourth member in the equation (1.1) is added. More detailed description of the mathematical method with many-dimensional polynomials including differentiation and integration will be published later in suitable journal according the description in JINR library with some implementations.

2.Description of the package OPEM2 (fortran 77)

It consists of main program CALIB and five subroutines: the principle TSTWO, ORTTWO, PRETWO and 2 auxiliary DEGREE, NUMBER. The main program CALIB gives:

M - the number of points (nodes);
X Y - the arrays of coordinates of points;
F - the values of function to be approximated in nodes;
W - the values of weights.

The principal subroutines contain a named common /LINKS/ where the recurrence factors and auxiliary variables are recorded. The present size of /LINKS/ allows for overall degrees as high as 16 and for Lmax=(16+1)(16+2)/2=153. The array dimensions in /LINKS/ may be changed.

The subroutine TSTWO(X,Y,W,M,F,A,NDEG) is the subroutine, which calls subprograms. First, it calls PRETWO for preparing the recurrence factors with overall degree NDEG. Afterwards it calculates the orthonormal coefficients A, calling the subroutine ORTTWO. After second call of ORTTWO it calculates the approximating values and deviations in the points, selects optimal degree Lr, selects min (xi)2, calculates the errors of series coefficients after third call of ORTTWO and then prints the results.

X Y - two given one-dimensional arrays of nodes;
W - given weights in every node;
M - given number of points;
F - given function in two dimensions;
A - orthonormal coefficients to be calculated;
NDEG - the max overall degree of polynomials.

The subroutine ORTTWO(NUMBMX,X,Y,POLY) returns the specific values of all polynomials in an orthonormal series from 1 to NUMBMX at given values of x and y in POLY-array.

NUMBMX - the value of Lr - optimal degree of polynomials;
X Y - the coordinates of points, where polynomials are calculated;
POLY - one-dimensional array of results.

The subroutine PRETWO(M,MAXD,X,Y,W,POLY) prepares the recurrence factors of two-dimensional polynomials, orthonormal over a discrete point set in a single call, which must precede any working call of ORTTWO.

M - the number of points;
MAXD - the max overall degree for which recurrence factors are computed (the same as NDEG in TSTWO);
X Y - two one-dimensional arrays containing the coordinates of points;
W - weights;
POLY - the same as in ORTTWO, and here when computing betaL-1,alfarL and CL.

PRETWO makes successive recursive calls to ORTTWO with NUMBMX = L-1, then L is incremented, etc, until Lmax is reached.

The subroutine DEGREE(N,J,JX,JY) calculates at given ordinary number of polynomials N the overall degree of the leading term J, the degree in x-variable JX, and that in y-variable JY.

The integer function NUMBER(JX,JY) at given x degree JX and y degree JY returns the ordinary number of polynomials in a series according to Weisfeld [3] with particular stress on the leading term.


Some applications for receiving calibration coefficients between ideal and measured grid in the optical devices in HEP are given in [5,6,7]. The last application of the package is for H1 experiment in DESY, Hamburg - to look for the coordinates of the maximum of the surface formed by gaussians. The analytic surface over the 400 point's grid is our approximation.

It is done in collaboration with Dubna (LIT).


  1. A.Householder, Principles of numerical analysis (Mc Graw-Hill, New York,1953), p.221.
  2. G.Forsythe, J.Soc.Ind.Math. 5 (1957) 74.
  3. W.Weisfeld, Numer.Math. 1 (1959) 38.
  4. V.Gadjokov, Commun.JINR E11-80-782, Dubna (1980).
  5. N.Bogdanova, V.Gadjokov, Comp.Phys.Commun.,24 (1981) p.225-229.
  6. N.Bogdanova, J.Comp.and Appl.Mathematics, 14 (1986) p.345-351.
  7. N.Bogdanova, T.Kupenova, Comp.Phys.Commun.,73 (1992) p.170-178.

Sources of OPEM2 are submitted.

home up e-mail