Você não pode selecionar mais de 25 tópicos Os tópicos devem começar com uma letra ou um número, podem incluir traços ('-') e podem ter até 35 caracteres.
 
 
 
 
 

145 linhas
3.5 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. y = Y[K];
  30. x = X[K];
  31. p = 1.0;
  32. for (R = 0; R < M + 1; R++) {
  33. B[R] += y * p;
  34. p = p * x;
  35. }
  36. }
  37. /* Zero the array */
  38. for (J = 1; J <= 2 * M; J++)
  39. P[J] = 0;
  40. P[0] = N;
  41. /* Compute the sum of powers of x_(K-1) */
  42. for (K = 0; K < N; K++) {
  43. x = X[K];
  44. p = X[K];
  45. for (J = 1; J <= 2 * M; J++) {
  46. P[J] += p;
  47. p = p * x;
  48. }
  49. }
  50. /* Determine the matrix entries */
  51. for (R = 0; R < M + 1; R++) {
  52. for (K = 0; K < M + 1; K++)
  53. 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++)
  69. Row[J] = J;
  70. /* Start LU factorization */
  71. for (P = 0; P < N - 1; P++) {
  72. /* Find pivot element */
  73. for (K = P + 1; K < N; K++) {
  74. if (fabs(A[Row[K]][P]) > fabs(A[Row[P]][P])) {
  75. /* Switch the index for the p-1 th pivot row if necessary */
  76. T = Row[P];
  77. Row[P] = Row[K];
  78. Row[K] = T;
  79. DET = -DET;
  80. }
  81. } /* End of simulated row interchange */
  82. if (A[Row[P]][P] == 0) {
  83. /* The matrix is SINGULAR ! */
  84. return;
  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++)
  100. Y[K] = B[K];
  101. Y[0] = B[Row[0]];
  102. for (K = 1; K < N; K++) {
  103. SUM = 0;
  104. for (C = 0; C <= K - 1; C++)
  105. SUM += A[Row[K]][C] * Y[C];
  106. Y[K] = B[Row[K]] - SUM;
  107. }
  108. /* Start the back substitution */
  109. X[N - 1] = Y[N - 1] / A[Row[N - 1]][N - 1];
  110. for (K = N - 2; K >= 0; K--) {
  111. SUM = 0;
  112. for (C = K + 1; C < N; C++) {
  113. SUM += A[Row[K]][C] * X[C];
  114. }
  115. X[K] = (Y[K] - SUM) / A[Row[K]][K];
  116. } /* End of back substitution */
  117. /* Output */
  118. for (K = 0; K < N; K++)
  119. Cf[K] = X[K];
  120. }