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.
 
 
 
 
 

152 lines
3.9 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
  15. polyreg (const int M, const int N, const double X[], const double Y[],
  16. double C[])
  17. {
  18. int R, K, J; /* Loop counters */
  19. double A[DMAX][DMAX]; /* A */
  20. double B[DMAX];
  21. double P[2 * DMAX + 1];
  22. double x, y;
  23. double p;
  24. /* Zero the array */
  25. for (R = 0; R < M + 1; R++)
  26. B[R] = 0;
  27. /* Compute the column vector */
  28. for (K = 0; K < N; K++)
  29. {
  30. y = Y[K];
  31. x = X[K];
  32. p = 1.0;
  33. for (R = 0; R < M + 1; R++)
  34. {
  35. B[R] += y * p;
  36. p = p * x;
  37. }
  38. }
  39. /* Zero the array */
  40. for (J = 1; J <= 2 * M; J++)
  41. P[J] = 0;
  42. P[0] = N;
  43. /* Compute the sum of powers of x_(K-1) */
  44. for (K = 0; K < N; K++)
  45. {
  46. x = X[K];
  47. p = X[K];
  48. for (J = 1; J <= 2 * M; J++)
  49. {
  50. P[J] += p;
  51. p = p * x;
  52. }
  53. }
  54. /* Determine the matrix entries */
  55. for (R = 0; R < M + 1; R++)
  56. {
  57. for (K = 0; K < M + 1; K++)
  58. A[R][K] = P[R + K];
  59. }
  60. /* Solve the linear system of M + 1 equations : A*C = B
  61. for the coefficient vector C = (c_1,c_2,..,c_M,c_(M+1)) */
  62. FactPiv (M + 1, A, B, C);
  63. } /* end main */
  64. /*--------------------------------------------------------*/
  65. static void
  66. FactPiv (int N, double A[DMAX][DMAX], double B[], double Cf[])
  67. {
  68. int K, P, C, J; /* Loop counters */
  69. int Row[NMAX]; /* Field with row-number */
  70. double X[DMAX], Y[DMAX];
  71. double SUM, DET = 1.0;
  72. int T;
  73. /* Initialize the pointer vector */
  74. for (J = 0; J < N; J++)
  75. Row[J] = J;
  76. /* Start LU factorization */
  77. for (P = 0; P < N - 1; P++) {
  78. /* Find pivot element */
  79. for (K = P + 1; K < N; K++) {
  80. if (fabs (A[Row[K]][P]) > fabs (A[Row[P]][P])) {
  81. /* Switch the index for the p-1 th pivot row if necessary */
  82. T = Row[P];
  83. Row[P] = Row[K];
  84. Row[K] = T;
  85. DET = -DET;
  86. }
  87. } /* End of simulated row interchange */
  88. if (A[Row[P]][P] == 0) {
  89. printf ("The matrix is SINGULAR !\n");
  90. printf ("Cannot use algorithm --> exit\n");
  91. exit (1);
  92. }
  93. /* Multiply the diagonal elements */
  94. DET = DET * A[Row[P]][P];
  95. /* Form multiplier */
  96. for (K = P + 1; K < N; K++) {
  97. A[Row[K]][P] = A[Row[K]][P] / A[Row[P]][P];
  98. /* Eliminate X_(p-1) */
  99. for (C = P + 1; C < N + 1; C++) {
  100. A[Row[K]][C] -= A[Row[K]][P] * A[Row[P]][C];
  101. }
  102. }
  103. } /* End of L*U factorization routine */
  104. DET = DET * A[Row[N - 1]][N - 1];
  105. /* Start the forward substitution */
  106. for (K = 0; K < N; K++)
  107. Y[K] = B[K];
  108. Y[0] = B[Row[0]];
  109. for (K = 1; K < N; K++) {
  110. SUM = 0;
  111. for (C = 0; C <= K - 1; C++)
  112. SUM += A[Row[K]][C] * Y[C];
  113. Y[K] = B[Row[K]] - SUM;
  114. }
  115. /* Start the back substitution */
  116. X[N - 1] = Y[N - 1] / A[Row[N - 1]][N - 1];
  117. for (K = N - 2; K >= 0; K--) {
  118. SUM = 0;
  119. for (C = K + 1; C < N; C++) {
  120. SUM += A[Row[K]][C] * X[C];
  121. }
  122. X[K] = (Y[K] - SUM) / A[Row[K]][K];
  123. } /* End of back substitution */
  124. /* Output */
  125. for (K = 0; K < N; K++)
  126. Cf[K] = X[K];
  127. }