|
@@ -0,0 +1,181 @@
|
|
|
+#include<stdio.h>
|
|
|
+#include<stdlib.h>
|
|
|
+#include<math.h>
|
|
|
+
|
|
|
+// Необязательный флаг для дебага (можно удалить)
|
|
|
+#define DEBUG 0
|
|
|
+
|
|
|
+// Фиксированная размерность матрицы (обрабатывается препроцессором, просто alias)
|
|
|
+#define N 10
|
|
|
+
|
|
|
+// Тип вектора
|
|
|
+typedef double vec[N];
|
|
|
+
|
|
|
+// Тип матрицы
|
|
|
+typedef double mat[N][N];
|
|
|
+
|
|
|
+// Структура для удобного сбора данных
|
|
|
+struct result {
|
|
|
+ double lambda;
|
|
|
+ int M;
|
|
|
+};
|
|
|
+
|
|
|
+// Скалярное произведение
|
|
|
+double dot(vec v, vec u) {
|
|
|
+ double res = 0;
|
|
|
+ for(int i = 0; i < N; i++) res += v[i]*u[i];
|
|
|
+
|
|
|
+ return res;
|
|
|
+}
|
|
|
+
|
|
|
+// Прибавление вектора
|
|
|
+void add(vec v, vec u) {
|
|
|
+ for(int i = 0; i < N; i++) {
|
|
|
+ v[i] += u[i];
|
|
|
+ }
|
|
|
+}
|
|
|
+
|
|
|
+// Функция для вычисления сдвига
|
|
|
+void shift(mat m, double s) {
|
|
|
+ for(int i = 0; i < N; i++) {
|
|
|
+ m[i][i] = m[i][i] - s;
|
|
|
+ }
|
|
|
+}
|
|
|
+
|
|
|
+// Умножение матрицы на вектор (res = m*v)
|
|
|
+void apply(mat m, vec v, vec res) {
|
|
|
+ for(int i = 0; i < N; i++) {
|
|
|
+ res[i] = 0;
|
|
|
+ for(int j = 0; j < N; j++) {
|
|
|
+ res[i] += m[i][j] * v[j];
|
|
|
+ }
|
|
|
+ }
|
|
|
+}
|
|
|
+
|
|
|
+// Функция нормализации
|
|
|
+void normalize(vec v){
|
|
|
+ double norm = sqrt(dot(v, v));
|
|
|
+ for(int i = 0; i < N; i++) {
|
|
|
+ v[i] /= norm;
|
|
|
+ }
|
|
|
+}
|
|
|
+
|
|
|
+// Функция для вычисления ошибки между векторами
|
|
|
+double error_vector(vec x, vec x1) {
|
|
|
+ double error = 0;
|
|
|
+ for (int i = 0; i < N; i++) {
|
|
|
+ error = fmax(error, fabs(fabs(x[i]) - fabs(x1[i])));
|
|
|
+ }
|
|
|
+ return error;
|
|
|
+}
|
|
|
+
|
|
|
+// Копирование вектора
|
|
|
+void copy_vec(vec v, vec copy) {
|
|
|
+ for(int i = 0; i < N; i++) {
|
|
|
+ copy[i] = v[i];
|
|
|
+ }
|
|
|
+}
|
|
|
+
|
|
|
+// Чтение матрицы из файла
|
|
|
+void mat_from_file(FILE* f, mat m) {
|
|
|
+ double d;
|
|
|
+ for (int i = 0; i < N; i++) {
|
|
|
+ for (int j = 0; j < N; j++) {
|
|
|
+ fscanf(f, "%lf,", &d);
|
|
|
+ m[i][j] = d;
|
|
|
+ }
|
|
|
+ }
|
|
|
+}
|
|
|
+
|
|
|
+// Чтение вектора из файла
|
|
|
+void vec_from_file(FILE* f, vec v) {
|
|
|
+ double d;
|
|
|
+ for (int i = 0; i < N; i++) {
|
|
|
+ fscanf(f, "%lf,", &d);
|
|
|
+ v[i] = d;
|
|
|
+ }
|
|
|
+}
|
|
|
+
|
|
|
+// Красивая печать матрицы
|
|
|
+void print_mat(mat m) {
|
|
|
+ for(int j = 0; j < N; j++) {
|
|
|
+ printf("-----------|");
|
|
|
+ }
|
|
|
+ printf("\n");
|
|
|
+ for(int i = 0; i < N; i++) {
|
|
|
+ for(int j = 0; j < N; j++) {
|
|
|
+ printf("%10.4lf |", m[i][j]);
|
|
|
+ }
|
|
|
+ printf("\n");
|
|
|
+ for(int j = 0; j < N; j++) {
|
|
|
+ printf("-----------|");
|
|
|
+ }
|
|
|
+ printf("\n");
|
|
|
+ }
|
|
|
+}
|
|
|
+
|
|
|
+// Красивая печать вектора
|
|
|
+void print_vec(vec v) {
|
|
|
+ for(int j = 0; j < N; j++) {
|
|
|
+ printf("%10.4lf |", v[j]);
|
|
|
+ }
|
|
|
+ printf("\n");
|
|
|
+}
|
|
|
+
|
|
|
+struct result solve(mat m, vec x, double shift_v, double eps, vec eigen_vector) {
|
|
|
+ struct result res;
|
|
|
+ shift(m, shift_v);
|
|
|
+ res.M = 0;
|
|
|
+ vec y;
|
|
|
+ double lambda2;
|
|
|
+ normalize(x);
|
|
|
+ for (int i = 0; i < 100; i++) {
|
|
|
+ apply(m, x, y);
|
|
|
+ res.lambda = shift_v + dot(x, y)/dot(x,x);
|
|
|
+ res.M++;
|
|
|
+ normalize(y);
|
|
|
+ copy_vec(y, eigen_vector);
|
|
|
+ copy_vec(y, x);
|
|
|
+ if (fabs((res.lambda-lambda2)/lambda2) < eps){
|
|
|
+ break;
|
|
|
+ }
|
|
|
+ lambda2 = res.lambda;
|
|
|
+ }
|
|
|
+ return res;
|
|
|
+}
|
|
|
+
|
|
|
+
|
|
|
+int main(){
|
|
|
+ vec x;
|
|
|
+ FILE* f;
|
|
|
+ if ((f = fopen("matrix.txt", "r")) == NULL) return 1;
|
|
|
+ FILE* fs;
|
|
|
+ if ((fs = fopen("value_shift.txt", "r")) == NULL) return 1;
|
|
|
+ FILE* vf;
|
|
|
+ if ((vf = fopen("vector_shift.txt", "r")) == NULL) return 1;
|
|
|
+
|
|
|
+ FILE* file;
|
|
|
+ if ((file = fopen("lab_6_error.csv", "w+")) == NULL) return 1;
|
|
|
+ mat m;
|
|
|
+ vec eigen_vector;
|
|
|
+ for (int i = 0; i < 9; i++) {
|
|
|
+ for(int j = 0; j < N; j++) {
|
|
|
+ x[j] = 1;
|
|
|
+ }
|
|
|
+ mat_from_file(f, m);
|
|
|
+ vec_from_file(vf, eigen_vector);
|
|
|
+
|
|
|
+ double shift_v;
|
|
|
+ fscanf(fs, "%lf\n", &shift_v);
|
|
|
+
|
|
|
+ vec eigen_approximate;
|
|
|
+ struct result solu = solve(m, x, shift_v, pow(10,-10), eigen_approximate);
|
|
|
+ print_vec(eigen_vector);
|
|
|
+ print_vec(eigen_approximate);
|
|
|
+ printf("\n");
|
|
|
+ fprintf(file, "%e, %e, %e\n", pow(10,(i-2)), fabs(pow(10,i-1) - solu.lambda)/pow(10,i-1), error_vector(eigen_vector,eigen_approximate));
|
|
|
+ }
|
|
|
+ fclose(f);
|
|
|
+ fclose(fs);
|
|
|
+ fclose(file);
|
|
|
+}
|