|
- /* ---------------------------------------------------------------------------
-
- Polynomial regression, freely adapted from:
-
- NUMERICAL METHODS: C Programs, (c) John H. Mathews 1995
- Algorithm translated to C by: Dr. Norman Fahrer
- NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
- Prentice Hall, International Editions: ISBN 0-13-625047-5
- This free software is compliments of the author.
- E-mail address: in%"mathews@fullerton.edu"
- */
-
- #include <math.h>
-
- #define DMAX 5 /* Maximum degree of polynomial */
- #define NMAX 10 /* Maximum number of points */
-
- static void FactPiv(int N, double A[DMAX][DMAX], double B[], double Cf[]);
-
- void polyreg(const int M, const int N, const double X[], const double Y[], double C[]) {
- int R, K, J; /* Loop counters */
- double A[DMAX][DMAX]; /* A */
- double B[DMAX];
- double P[2 * DMAX + 1];
- double x, y;
- double p;
-
- /* Zero the array */
- for (R = 0; R < M + 1; R++)
- B[R] = 0;
-
- /* Compute the column vector */
- for (K = 0; K < N; K++) {
- y = Y[K];
- x = X[K];
- p = 1.0;
- for (R = 0; R < M + 1; R++) {
- B[R] += y * p;
- p = p * x;
- }
- }
-
- /* Zero the array */
- for (J = 1; J <= 2 * M; J++)
- P[J] = 0;
-
- P[0] = N;
-
- /* Compute the sum of powers of x_(K-1) */
- for (K = 0; K < N; K++) {
- x = X[K];
- p = X[K];
- for (J = 1; J <= 2 * M; J++) {
- P[J] += p;
- p = p * x;
- }
- }
-
- /* Determine the matrix entries */
- for (R = 0; R < M + 1; R++) {
- for (K = 0; K < M + 1; K++)
- A[R][K] = P[R + K];
- }
-
- /* Solve the linear system of M + 1 equations: A*C = B
- for the coefficient vector C = (c_1,c_2,..,c_M,c_(M+1)) */
- FactPiv(M + 1, A, B, C);
- } /* end main */
-
-
- /*--------------------------------------------------------*/
- static void FactPiv(int N, double A[DMAX][DMAX], double B[], double Cf[]) {
- int K, P, C, J; /* Loop counters */
- int Row[NMAX]; /* Field with row-number */
- double X[DMAX], Y[DMAX];
- double SUM, DET = 1.0;
- int T;
-
- /* Initialize the pointer vector */
- for (J = 0; J < NMAX; J++)
- Row[J] = J;
-
- /* Start LU factorization */
- for (P = 0; P < N - 1; P++) {
- /* Find pivot element */
- for (K = P + 1; K < N; K++) {
- if (fabs(A[Row[K]][P]) > fabs(A[Row[P]][P])) {
- /* Switch the index for the p-1 th pivot row if necessary */
- T = Row[P];
- Row[P] = Row[K];
- Row[K] = T;
- DET = -DET;
- }
- } /* End of simulated row interchange */
- if (A[Row[P]][P] == 0) {
- /* The matrix is SINGULAR! */
- return;
- }
-
- /* Multiply the diagonal elements */
- DET = DET * A[Row[P]][P];
-
- /* Form multiplier */
- for (K = P + 1; K < N; K++) {
- A[Row[K]][P] = A[Row[K]][P] / A[Row[P]][P];
-
- /* Eliminate X_(p-1) */
- for (C = P + 1; C < N + 1; C++) {
- A[Row[K]][C] -= A[Row[K]][P] * A[Row[P]][C];
- }
- }
- }
- /* End of L*U factorization routine */
- DET = DET * A[Row[N - 1]][N - 1];
-
- /* Start the forward substitution */
- for (K = 0; K < N; K++)
- Y[K] = B[K];
-
- Y[0] = B[Row[0]];
-
- for (K = 1; K < N; K++) {
- SUM = 0;
- for (C = 0; C <= K - 1; C++)
- SUM += A[Row[K]][C] * Y[C];
- Y[K] = B[Row[K]] - SUM;
- }
-
- /* Start the back substitution */
- X[N - 1] = Y[N - 1] / A[Row[N - 1]][N - 1];
- for (K = N - 2; K >= 0; K--) {
- SUM = 0;
- for (C = K + 1; C < N; C++) {
- SUM += A[Row[K]][C] * X[C];
- }
- X[K] = (Y[K] - SUM) / A[Row[K]][K];
- } /* End of back substitution */
-
- /* Output */
- for (K = 0; K < N; K++)
- Cf[K] = X[K];
- }
|