Seidel

#include "stdafx.h" #include <cmath> #include <iostream> #include <fstream> #include <iomanip> using namespace std; const int N = 2; const double eps = 1e-5; void Id_matrix(double **A, int k) { for (int i = 0; i < k; i++) { for (int j = 0; j < N + 1; j++) { A[i][j] = 0; A[j][i] = 0; } } for (int i = 0; i < k; i++) { for (int j = 0; j < k; j++) { if (i == j) A[i][j] = 1; } } } void Id_vector(double *vec, int zeros, int j) //зануляет первые zeros элементов, на j-е место ставит 1 { for (int i = 0; i < zeros; i++) vec[i] = 0; vec[j] = 1; } void Mul_matmat(double **A, double **B, double **Res) { for (int i = 0; i < N + 1; i++) { for (int j = 0; j < N + 1; j++) { for (int k = 0; k < N + 1; ++k) { Res[i][j] += A[i][k] * B[k][j]; } } } } void Mul_matvec(double **A, double *b, double *res) { for (int i = 0; i < N + 1; i++) { for (int j = 0; j < N + 1; j++) { res[i] += A[i][j] * b[j]; } } } double vTu(double *vec1, double *vec2) { double res = 0; for (int i = 0; i < N + 1; i++) { res += vec1[i] * vec2[i]; } return res; } void vuT(double *vec1, double *vec2, double **res) { for (int i = 0; i < N + 1; i++) { for (int j = 0; j < N + 1; j++) { res[i][j] = vec1[i] * vec2[j]; } } } void Mul_numvec(double *vec, double m) { for (int j = 0; j < N + 1; j++) { vec[j] *= m; } } void Mul_nummat(double **A, double m) { for (int i = 0; i < N + 1; i++) { for (int j = 0; j < N + 1; j++) { A[i][j] *= m; } } } double vec_norm(double *vec) { double norm = 0; for (int i = 0; i < N + 1; i++) norm += pow(vec[i], 2); return sqrt(norm); } bool STOP(double **A, double *x_k, double *b, double eps) //критерий остановки - ||Ax_k-b||<eps { bool stop = false; double *Ax = new double[N + 1]; double *err = new double[N + 1]; Mul_matvec(A, x_k, Ax); for (int i = 0; i < N + 1; i++) err[i] = Ax[i] - b[i]; if (vec_norm(err) < eps) stop = true; return stop; } void seidel(double **A, double *b, double *res) { double *x_prev = new double[N + 1]; double *x_cur = new double[N + 1]; for (int i = 0; i < N + 1; i++) { x_cur[i] = b[i]; } do { for (int i = 0; i < N + 1; i++) x_prev[i] = x_cur[i]; for (int i = 0; i < N + 1; i++) { double var = 0; for (int j = 0; j < i; j++) var += (A[i][j] * x_cur[j]); for (int j = i + 1; j < N + 1; j++) var += (A[i][j] * x_prev[j]); x_cur[i] = (b[i] - var) / A[i][i]; } } while (!STOP(A, x_cur, b, eps)); for (int i = 0; i < N + 1; i++) res[i] = x_cur[i]; } int main() { double ** A = new double *[N + 1];//объявление матрицы, которую можно подать как аргумент функции for (int i = 0; i < N+1; i++) A[i] = new double[N+1]; Id_matrix(A, N + 1); Mul_nummat(A, 4); double *b = new double[N + 1];//объявление вектора, который можно подать как аргумент функции for (int i = 0; i < N + 1; i++) b[i] = i + 1; double *res = new double[N + 1]; //seidel(A, b, res); //for (int i = 0; i < N + 1; i++) cout << res[i] << " "; cout << endl; system("pause"); return 0; }