OPOLY
Generates orthogonal polynomials with respect to x‑values and specified weights.
Required Arguments
X — Vector of length N containing the x‑values. (Input)
NDEG — Degree of highest degree orthogonal polynomial to be generated. (Input)
SMULTC — Multiplicative constant used to compute a scaled version of x on the interval ‑2 to 2, inclusive. (Output)
SADDC — Additive constant used to compute a scaled version of x on the interval ‑2 to 2, inclusive. (Output)
SX — Vector of length N containing the scaled version of x on the interval ‑2 to 2, inclusive, computed as follows: SX(i) = SMULTC * X(i) + SADDC where i = 1, 2, …, N. (Output)
If X is not needed, SX and X can occupy the same storage locations.
A — Vector of length NDEG containing constants used to generate orthogonal polynomials. (Output)
B — Vector of length NDEG containing constants used to generate orthogonal polynomials. (Output)
POLY — Matrix, N by NDEG, containing the orthogonal polynomials evaluated at SX(i) for i = 1, 2, …, N. (Output)
Optional Arguments
N — Number of x‑values. (Input)
Default: N = size (POLY,1).
IWT — Weighting option. (Input)
IWT = 0 means that all weights are 1.0. For IWT = 1, WT contains the weights.
Default: IWT = 0.
WT — Vector of length N containing the weights. (Input, if IWT = 1)
If IWT = 0, WT is not referenced and can be a vector of length one.
LDPOLY — Leading dimension of POLY exactly as specified in the dimension statement in the calling program. (Input)
Default: LDPOLY = size (POLY,1).
FORTRAN 90 Interface
Generic: CALL OPOLY (X, NDEG, SMULTC, SADDC, SX, A, B, POLY [, …])
Specific: The specific interface names are S_OPOLY and D_OPOLY.
FORTRAN 77 Interface
Single: CALL OPOLY (N, X, IWT, WT, NDEG, SMULTC, SADDC, SX, A, B, POLY, LDPOLY)
Double: The double precision name is DOPOLY.
Description
Routine OPOLY generates orthogonal polynomials over a set of xi’s and with respect to weights wi. The routine OPOLY is based on the algorithm of Forsythe (1957). (See also Kennedy and Gentle 1980.) A modification to Forsythe’s algorithm is made for the inclusion of weights (Kelly 1967, page 68).
Let xi be a value of the independent variable. The xi’s are scaled to the interval [‑2, 2] for computational accuracy. The scaled version of the independent variable is computed by the formula zi = mxi + c. The multiplicative scaling constant m (stored in SMULTC) is
The additive constant c (stored in SADDC) is
Orthogonal polynomials are generated using the three-term recurrence relationship
beginning with the initial polynomials
The aj’s and bj’s (stored in A and B) are computed to make the pj(z)’s orthogonal, with respect to the the set of weights wi, and over the set zi.
Comments
1. Informational error
|
Type |
Code |
Description |
|
3 |
8 |
N must be greater than NDEG in order for higher order polynomials to be nonzero. Columns N + 1 through NDEG of POLY are set to zero. |
2. The orthogonal polynomials evaluated at each scaled X value are computed from A and B as follows:
POLY(I, 1) = SX(I) ‑ A(1)
POLY(I, 2) = (SX(I) ‑ A(2)) * POLY(I, 1) ‑ B(2)
POLY(I, J) = (SX(I) ‑ A(J)) * POLY(I, J ‑ 1) ‑ B(J) * POLY(I, J ‑ 2) for J = 3 through NDEG.
3. If NDEG is greater than 10, the accuracy of the results may be questionable.
Example
First-degree and second-degree orthogonal polynomials are generated using equally spaced x values 1, 2, …, 12. (Equally spaced x values are not required by OPOLY.)
USE OPOLY_INT
USE UMACH_INT
USE WRRRN_INT
IMPLICIT NONE
INTEGER LDPOLY, N, NDEG
PARAMETER (N=12, NDEG=2, LDPOLY=N)
!
INTEGER NOUT
REAL A(NDEG), B(NDEG), POLY(LDPOLY,NDEG), SADDC, SMULTC, &
SX(N), WT(1), X(N)
!
DATA X/1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, &
12.0/
!
CALL OPOLY (X, NDEG, SMULTC, SADDC, SX, A, B, POLY)
CALL UMACH (2, NOUT)
WRITE (NOUT,99999) SMULTC, SADDC
99999 FORMAT (' SMULTC = ', F7.3, ' SADDC = ', F7.3)
CALL WRRRN ('A', A, 1, NDEG, 1)
CALL WRRRN ('B', B, 1, NDEG, 1)
CALL WRRRN ('POLY', POLY)
END
Output
SMULTC = 0.364 SADDC = -2.364
A
1 2
-5.960E-08 -1.009E-07
B
1 2
0.000 1.576
POLY
1 2
1 -2.000 2.424
2 -1.636 1.102
3 -1.273 0.044
4 -0.909 -0.749
5 -0.545 -1.278
6 -0.182 -1.543
7 0.182 -1.543
8 0.545 -1.278
9 0.909 -0.749
10 1.273 0.044
11 1.636 1.102
12 2.000 2.424
