You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
 
 
 
 
 

143 lines
3.3 KiB

  1. /* ---------------------------------------------------------------------------
  2. Polynomial regression, freely adapted from:
  3. NUMERICAL METHODS: C Programs, (c) John H. Mathews 1995
  4. Algorithm translated to C by: Dr. Norman Fahrer
  5. NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
  6. Prentice Hall, International Editions: ISBN 0-13-625047-5
  7. This free software is compliments of the author.
  8. E-mail address: in%"mathews@fullerton.edu"
  9. */
  10. #include <math.h>
  11. #define DMAX 5 /* Maximum degree of polynomial */
  12. #define NMAX 10 /* Maximum number of points */
  13. static void FactPiv(int N, double A[DMAX][DMAX], double B[], double Cf[]);
  14. void polyreg(const int M, const int N, const double X[], const double Y[], double C[]) {
  15. int R, K, J; /* Loop counters */
  16. double A[DMAX][DMAX]; /* A */
  17. double B[DMAX];
  18. double P[2 * DMAX + 1];
  19. double x, y;
  20. double p;
  21. /* Zero the array */
  22. for (R = 0; R < M + 1; R++)
  23. B[R] = 0;
  24. /* Compute the column vector */
  25. for (K = 0; K < N; K++) {
  26. y = Y[K];
  27. x = X[K];
  28. p = 1.0;
  29. for (R = 0; R < M + 1; R++) {
  30. B[R] += y * p;
  31. p = p * x;
  32. }
  33. }
  34. /* Zero the array */
  35. for (J = 1; J <= 2 * M; J++)
  36. P[J] = 0;
  37. P[0] = N;
  38. /* Compute the sum of powers of x_(K-1) */
  39. for (K = 0; K < N; K++) {
  40. x = X[K];
  41. p = X[K];
  42. for (J = 1; J <= 2 * M; J++) {
  43. P[J] += p;
  44. p = p * x;
  45. }
  46. }
  47. /* Determine the matrix entries */
  48. for (R = 0; R < M + 1; R++) {
  49. for (K = 0; K < M + 1; K++)
  50. A[R][K] = P[R + K];
  51. }
  52. /* Solve the linear system of M + 1 equations: A*C = B
  53. for the coefficient vector C = (c_1,c_2,..,c_M,c_(M+1)) */
  54. FactPiv(M + 1, A, B, C);
  55. } /* end main */
  56. /*--------------------------------------------------------*/
  57. static void FactPiv(int N, double A[DMAX][DMAX], double B[], double Cf[]) {
  58. int K, P, C, J; /* Loop counters */
  59. int Row[NMAX]; /* Field with row-number */
  60. double X[DMAX], Y[DMAX];
  61. double SUM, DET = 1.0;
  62. int T;
  63. /* Initialize the pointer vector */
  64. for (J = 0; J < NMAX; J++)
  65. Row[J] = J;
  66. /* Start LU factorization */
  67. for (P = 0; P < N - 1; P++) {
  68. /* Find pivot element */
  69. for (K = P + 1; K < N; K++) {
  70. if (fabs(A[Row[K]][P]) > fabs(A[Row[P]][P])) {
  71. /* Switch the index for the p-1 th pivot row if necessary */
  72. T = Row[P];
  73. Row[P] = Row[K];
  74. Row[K] = T;
  75. DET = -DET;
  76. }
  77. } /* End of simulated row interchange */
  78. if (A[Row[P]][P] == 0) {
  79. /* The matrix is SINGULAR! */
  80. return;
  81. }
  82. /* Multiply the diagonal elements */
  83. DET = DET * A[Row[P]][P];
  84. /* Form multiplier */
  85. for (K = P + 1; K < N; K++) {
  86. A[Row[K]][P] = A[Row[K]][P] / A[Row[P]][P];
  87. /* Eliminate X_(p-1) */
  88. for (C = P + 1; C < N + 1; C++) {
  89. A[Row[K]][C] -= A[Row[K]][P] * A[Row[P]][C];
  90. }
  91. }
  92. }
  93. /* End of L*U factorization routine */
  94. DET = DET * A[Row[N - 1]][N - 1];
  95. /* Start the forward substitution */
  96. for (K = 0; K < N; K++)
  97. Y[K] = B[K];
  98. Y[0] = B[Row[0]];
  99. for (K = 1; K < N; K++) {
  100. SUM = 0;
  101. for (C = 0; C <= K - 1; C++)
  102. SUM += A[Row[K]][C] * Y[C];
  103. Y[K] = B[Row[K]] - SUM;
  104. }
  105. /* Start the back substitution */
  106. X[N - 1] = Y[N - 1] / A[Row[N - 1]][N - 1];
  107. for (K = N - 2; K >= 0; K--) {
  108. SUM = 0;
  109. for (C = K + 1; C < N; C++) {
  110. SUM += A[Row[K]][C] * X[C];
  111. }
  112. X[K] = (Y[K] - SUM) / A[Row[K]][K];
  113. } /* End of back substitution */
  114. /* Output */
  115. for (K = 0; K < N; K++)
  116. Cf[K] = X[K];
  117. }