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.
 
 
 
 
 

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