123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208 |
- #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] -= 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_mat(mat m, mat copy) {
- for(int i = 0; i < N; i++) {
- for(int j = 0; j < N; j++) {
- copy[i][j] = m[i][j];
- }
- }
- }
- // Копирование вектора
- 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 < 10000; 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("matrix2.txt", "r")) == NULL) return 1;
- FILE* fs;
- if ((fs = fopen("value_shift2.txt", "r")) == NULL) return 1;
- FILE* vf;
- if ((vf = fopen("vector_shift2.txt", "r")) == NULL) return 1;
- FILE* file;
- if ((file = fopen("lab_6_M_low.csv", "w+")) == NULL) return 1;
- FILE* file2;
- if ((file2 = fopen("lab_6_M_high.csv", "w+")) == NULL) return 1;
- mat m;
- vec eigen_vector;
- double shift_v;
- fscanf(fs, "%lf\n", &shift_v);
- mat_from_file(f, m);
- vec_from_file(vf, eigen_vector);
- print_mat(m);
- vec eigen_approximate;
- for (int i = 2; i < 16; i++) {
- mat m_copy;
- copy_mat(m,m_copy);
- for(int j = 0; j < N; j++) {
- x[j] = 1;
- }
- struct result solu = solve(m_copy, x, shift_v, pow(10,(-i)), eigen_approximate);
- print_vec(eigen_vector);
- print_vec(eigen_approximate);
- printf("\n");
- fprintf(file, "%e, %i\n", pow(10,(-i)), solu.M);
- }
- fscanf(fs, "%lf\n", &shift_v);
- mat_from_file(f, m);
- print_mat(m);
- vec_from_file(vf, eigen_vector);
- for (int i = 2; i < 16; i++) {
- mat m_copy;
- copy_mat(m,m_copy);
- for(int j = 0; j < N; j++) {
- x[j] = j;
- }
- struct result solu = solve(m_copy, x, shift_v, pow(10,(-i)), eigen_approximate);
- print_vec(eigen_vector);
- print_vec(eigen_approximate);
- printf("\n");
- fprintf(file2, "%e, %i\n", pow(10,(-i)), solu.M);
- }
- fclose(f);
- fclose(fs);
- fclose(file);
- }
|