123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261 |
- #define _USE_MATH_DEFINES
- #include<stdio.h>
- #include<stdlib.h>
- #include<math.h>
- #include<time.h>
- // Необязательный флаг для дебага (можно удалить)
- #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);
- }
|