lab3.c 5.0 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208
  1. #include<stdio.h>
  2. #include<stdlib.h>
  3. #include<math.h>
  4. // Необязательный флаг для дебага (можно удалить)
  5. #define DEBUG 0
  6. // Фиксированная размерность матрицы (обрабатывается препроцессором, просто alias)
  7. #define N 10
  8. // Тип вектора
  9. typedef double vec[N];
  10. // Тип матрицы
  11. typedef double mat[N][N];
  12. // Структура для удобного сбора данных
  13. struct result {
  14. double lambda;
  15. int M;
  16. };
  17. // Скалярное произведение
  18. double dot(vec v, vec u) {
  19. double res = 0;
  20. for(int i = 0; i < N; i++) res += v[i]*u[i];
  21. return res;
  22. }
  23. // Прибавление вектора
  24. void add(vec v, vec u) {
  25. for(int i = 0; i < N; i++) {
  26. v[i] += u[i];
  27. }
  28. }
  29. // Функция для вычисления сдвига
  30. void shift(mat m, double s) {
  31. for(int i = 0; i < N; i++) {
  32. m[i][i] -= s;
  33. }
  34. }
  35. // Умножение матрицы на вектор (res = m*v)
  36. void apply(mat m, vec v, vec res) {
  37. for(int i = 0; i < N; i++) {
  38. res[i] = 0;
  39. for(int j = 0; j < N; j++) {
  40. res[i] += m[i][j] * v[j];
  41. }
  42. }
  43. }
  44. // Функция нормализаци
  45. void normalize(vec v){
  46. double norm = sqrt(dot(v, v));
  47. for(int i = 0; i < N; i++) {
  48. v[i] /= norm;
  49. }
  50. }
  51. // Функция для вычисления ошибки между векторами
  52. double error_vector(vec x, vec x1) {
  53. double error = 0;
  54. for (int i = 0; i < N; i++) {
  55. error = fmax(error, fabs(fabs(x[i]) - fabs(x1[i])));
  56. }
  57. return error;
  58. }
  59. // Копия матрицы
  60. void copy_mat(mat m, mat copy) {
  61. for(int i = 0; i < N; i++) {
  62. for(int j = 0; j < N; j++) {
  63. copy[i][j] = m[i][j];
  64. }
  65. }
  66. }
  67. // Копирование вектора
  68. void copy_vec(vec v, vec copy) {
  69. for(int i = 0; i < N; i++) {
  70. copy[i] = v[i];
  71. }
  72. }
  73. // Чтение матрицы из файла
  74. void mat_from_file(FILE* f, mat m) {
  75. double d;
  76. for (int i = 0; i < N; i++) {
  77. for (int j = 0; j < N; j++) {
  78. fscanf(f, "%lf,", &d);
  79. m[i][j] = d;
  80. }
  81. }
  82. }
  83. // Чтение вектора из файла
  84. void vec_from_file(FILE* f, vec v) {
  85. double d;
  86. for (int i = 0; i < N; i++) {
  87. fscanf(f, "%lf,", &d);
  88. v[i] = d;
  89. }
  90. }
  91. // Красивая печать матрицы
  92. void print_mat(mat m) {
  93. for(int j = 0; j < N; j++) {
  94. printf("-----------|");
  95. }
  96. printf("\n");
  97. for(int i = 0; i < N; i++) {
  98. for(int j = 0; j < N; j++) {
  99. printf("%10.4lf |", m[i][j]);
  100. }
  101. printf("\n");
  102. for(int j = 0; j < N; j++) {
  103. printf("-----------|");
  104. }
  105. printf("\n");
  106. }
  107. }
  108. // Красивая печать вектора
  109. void print_vec(vec v) {
  110. for(int j = 0; j < N; j++) {
  111. printf("%10.4lf |", v[j]);
  112. }
  113. printf("\n");
  114. }
  115. // Реализация метода
  116. struct result solve(mat m, vec x, double shift_v, double eps, vec eigen_vector) {
  117. struct result res;
  118. shift(m, shift_v);
  119. res.M = 0;
  120. vec y;
  121. double lambda2;
  122. normalize(x);
  123. for (int i = 0; i < 10000; i++) {
  124. apply(m, x, y);
  125. res.lambda = shift_v + dot(x, y)/dot(x,x);
  126. res.M++;
  127. normalize(y);
  128. copy_vec(y, eigen_vector);
  129. copy_vec(y, x);
  130. if (fabs((res.lambda-lambda2)/lambda2) < eps){
  131. break;
  132. }
  133. lambda2 = res.lambda;
  134. }
  135. return res;
  136. }
  137. int main(){
  138. vec x;
  139. FILE* f;
  140. if ((f = fopen("matrix2.txt", "r")) == NULL) return 1;
  141. FILE* fs;
  142. if ((fs = fopen("value_shift2.txt", "r")) == NULL) return 1;
  143. FILE* vf;
  144. if ((vf = fopen("vector_shift2.txt", "r")) == NULL) return 1;
  145. FILE* file;
  146. if ((file = fopen("lab_6_M_low.csv", "w+")) == NULL) return 1;
  147. FILE* file2;
  148. if ((file2 = fopen("lab_6_M_high.csv", "w+")) == NULL) return 1;
  149. mat m;
  150. vec eigen_vector;
  151. double shift_v;
  152. fscanf(fs, "%lf\n", &shift_v);
  153. mat_from_file(f, m);
  154. vec_from_file(vf, eigen_vector);
  155. print_mat(m);
  156. vec eigen_approximate;
  157. for (int i = 2; i < 16; i++) {
  158. mat m_copy;
  159. copy_mat(m,m_copy);
  160. for(int j = 0; j < N; j++) {
  161. x[j] = 1;
  162. }
  163. struct result solu = solve(m_copy, x, shift_v, pow(10,(-i)), eigen_approximate);
  164. print_vec(eigen_vector);
  165. print_vec(eigen_approximate);
  166. printf("\n");
  167. fprintf(file, "%e, %i\n", pow(10,(-i)), solu.M);
  168. }
  169. fscanf(fs, "%lf\n", &shift_v);
  170. mat_from_file(f, m);
  171. print_mat(m);
  172. vec_from_file(vf, eigen_vector);
  173. for (int i = 2; i < 16; i++) {
  174. mat m_copy;
  175. copy_mat(m,m_copy);
  176. for(int j = 0; j < N; j++) {
  177. x[j] = j;
  178. }
  179. struct result solu = solve(m_copy, x, shift_v, pow(10,(-i)), eigen_approximate);
  180. print_vec(eigen_vector);
  181. print_vec(eigen_approximate);
  182. printf("\n");
  183. fprintf(file2, "%e, %i\n", pow(10,(-i)), solu.M);
  184. }
  185. fclose(f);
  186. fclose(fs);
  187. fclose(file);
  188. }