|
@@ -0,0 +1,261 @@
|
|
|
+#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);
|
|
|
+}
|