#define _USE_MATH_DEFINES #include #include #include #include // Необязательный флаг для дебага (можно удалить) #define DEBUG 0 // Фиксированная размерность матрицы (обрабатывается препроцессором, просто alias) #define N 10 // struct N_eps { double eps; int I; }; // Тип вектора typedef double vec[N]; // Тип матрицы typedef double mat[N][N]; // Скалярное произведение 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, vec result) { for(int i = 0; i < N; i++) { result[i] = v[i] + u[i]; } } // Вычитание вектора void sub(vec v, vec u, vec result) { for(int i = 0; i < N; i++) { result[i] = v[i] - u[i]; } } // Умножение на скаляр void mul(vec v, double a) { for(int i = 0; i < N; i++) { v[i] *= a; } } // Умножение матрицы на вектор (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]; } } } // Копирование матрицы (copy = m) void copy_vec(vec copy, vec x) { for(int i = 0; i < N; i++) { copy[i] = x[i]; } } //Ошибка между двумя векторами double error(vec a, vec b){ double error = 0; for (int i = 0; i < N; i++) { if (error < fabs(a[i] - b[i])) { error = fabs(a[i] - b[i]); } } return error; } // Чтение матрицы из файла 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 min_max_from_file(FILE* f, vec v) { double d; for (int i = 0; i < 2; 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"); } void mx_b(mat m, vec x, vec b, vec result){ apply(m,x,result); sub(b,result,result); } // Решение СЛАУ методом Ричардсона struct N_eps solve(mat m, double eps, int M, double min, double max) { struct N_eps result; result.I = 0; vec x = { 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9}; vec x1; vec b; vec var; vec x_k = { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}; apply(m,x_k,b); print_vec(b); for (int k = 0; k < M; k++) { double t = cos(M_PI * (2 * (k + 1) - 1) / (2 * M)); printf("\n t = %lf\n", t); double alpha = 2.0 / ((min + max) + ((max - min)*t)); mx_b(m, x, b, var); //print_vec(var); mul(var,alpha); //print_vec(var); add(x,var,x1); print_vec(x1); if (error(x,x1) < eps){ result.eps = error(x,x1); result.I += 1; break; } result.eps = error(x,x1); result.I += 1; copy_vec(x,x1); print_vec(x); } return result; } //взять из мат лаба int main(){ FILE* f; FILE* f2; if ((f = fopen("lab4_1.txt", "r")) == NULL) return 1; if ((f2 = fopen("min_max.txt", "r")) == NULL) return 1; mat m; vec min_max; FILE* file; FILE* file2; if ((file = fopen("lab_4_cond.csv", "w+")) == NULL) return 1; if ((file2 = fopen("lab_4_time.csv", "w+")) == NULL) return 1; for (int k = 0; k < 8; k++){ mat_from_file(f, m); min_max_from_file(f2, min_max); clock_t begin = clock(); struct N_eps result = solve(m, pow(10, -5), 1000, min_max[0], min_max[1]); clock_t end = clock(); printf("\n\n\n"); fprintf(file, "%.5lf, %e\n", pow(10, k), result.eps); fprintf(file2, "%.5lf, %e\n", pow(10, k), ((float)(end - begin) / CLOCKS_PER_SEC)); } fclose(f); fclose(f2); fclose(file); fclose(file2); FILE* fbad = fopen("lab_4_error_bad.txt", "r"); FILE* fminbad = fopen("min_max2.txt", "r"); FILE* fbade = fopen("lab_4_error_bad.csv", "w+"); FILE* fbadm = fopen("lab_4_M_bad.csv", "w+"); mat_from_file(fbad,m); min_max_from_file(fminbad, min_max); for (int i = 1; i < 6; i++) { struct N_eps result = solve(m, pow(10, -i), 1000, min_max[0], min_max[1]); fprintf(fbade, "%e, %e\n", pow(10, -i), result.eps); fprintf(fbadm, "%e, %d\n", pow(10, -i), result.I); } fclose(fbad); fclose(fminbad); fclose(fbade); fclose(fbadm); FILE* fgood = fopen("lab_4_error_good.txt", "r");\ FILE* fmingood = fopen("min_max3.txt", "r"); FILE* fgoode = fopen("lab_4_error_good.csv", "w+"); FILE* fgoodm = fopen("lab_4_M_good.csv", "w+"); mat_from_file(fgood,m); min_max_from_file(fmingood, min_max); for (int i = 1; i < 16; i++) { struct N_eps result = solve(m, pow(10, -i), 1000, min_max[0], min_max[1]); fprintf(fgoodm, "%e, %d\n", pow(10, -i), result.I); fprintf(fgoode, "%e, %e\n", pow(10, -i), result.eps); } fclose(fgood); fclose(fmingood); fclose(fgoode); fclose(fgoodm); FILE* fbad2 = fopen("lab_4_error_bad.txt", "r"); FILE* fminbad2 = fopen("min_max2.txt", "r"); FILE* fbadmt = fopen("lab_4_M_bad_tol.csv", "w+"); mat_from_file(fbad2,m); min_max_from_file(fminbad2, min_max); for (int i = 0; i < 41; i++) { struct N_eps result = solve(m, pow(10, -15), i, min_max[0], min_max[1]); fprintf(fbadmt, "%e, %d\n", result.eps, result.I); } fclose(fbad2); fclose(fminbad2); fclose(fbadmt); FILE* fgood2 = fopen("lab_4_error_good.txt", "r"); FILE* fmindoof2 = fopen("min_max3.txt", "r"); FILE* fgoodmt2 = fopen("lab_4_M_good_tol.csv", "w+"); mat_from_file(fgood2,m); min_max_from_file(fmindoof2, min_max); for (int i = 0; i < 41; i++) { struct N_eps result = solve(m, pow(10, -15), i, min_max[0], min_max[1]); fprintf(fgoodmt2, "%e, %d\n", result.eps, i); } fclose(fgood2); fclose(fmindoof2); fclose(fgoodmt2); }