#include "stdafx.h"
#include <cmath>
#include <iostream>
#include <fstream>
#include <iomanip>
using namespace std;
const double eps = 1e-4;
double vec_norm(double *vec, int dim)
{
double norm = 0;
for (int i = 0; i < dim; i++) norm += pow(vec[i], 2);
return sqrt(norm);
}
double residual(double **A, double *x, double *b, int dim)
{
double *vec = new double[dim];
for (int i = 0; i < dim; i++) vec[i] = -b[i];
for (int i = 0; i < dim; i++) for (int j = 0; j < dim; j++) vec[i] += A[i][j] * x[j];
return vec_norm(vec, dim);
delete[] vec;
}
double eigen_val(double **A, int dim)
{
int iter = 0;
double summ, err, lambda_max, norm_x;
double *x_prev = new double[dim];
double *x_cur = new double[dim];
double *x_cur_l = new double[dim];
double *x_prev_norm = new double[dim];
for (int i = 0; i < dim; i++) //начальное приближение x = (1,0,...,0)т
x_prev[i] = 0;
x_prev[0] = 1;
do
{
summ = 0;
for (int i = 0; i < dim; i++)
summ += x_prev[i] * x_prev[i];
norm_x = sqrt(summ);
for (int i = 0; i < dim; i++)
x_prev_norm[i] = x_prev[i] / norm_x;
for (int i = 0; i < dim; i++)
{//x(k+1) = A*x(k)
x_cur[i] = 0;
for (int j = 0; j < dim; j++)
x_cur[i] += A[i][j] * x_prev_norm[j];
}
summ = 0;
for (int i = 0; i < dim; i++)
summ += x_cur[i] * x_cur[i];
lambda_max = sqrt(summ); //
for (int i = 0; i < dim; i++)
x_prev[i] = x_cur[i];
for (int i = 0; i < dim; i++) x_cur_l[i] = lambda_max*x_cur[i];
err = residual(A, x_cur, x_cur_l, dim);
iter++;
} while (err>eps);
cout << "iter = " << iter << endl;
return lambda_max;
delete[] x_prev;
delete[] x_cur;
delete[] x_prev_norm;
delete[] x_cur_l;
}
int main()
{
cout << setprecision(8);
int dim = 0;
FILE *input_data;
input_data = fopen("C:\\TASK3\\test3.dat", "r");
// input_data = fopen("H:\\task3\\test3.dat", "r"); В УНИВЕРЕ
fscanf(input_data, "%i", &dim);
cout << "for N = "<< dim << ", eps = " << eps << endl;
double ** A = new double *[dim];
for (int i = 0; i < dim; i++)
A[i] = new double[dim];
double *b = new double[dim];
cout << endl;
cout << " A: " << endl;
for (int i = 0; i < dim; i++) for (int j = 0; j < dim; j++) fscanf(input_data, "%Lf", &A[i][j]);
for (int i = 0; i < dim; i++) fscanf(input_data, "%Lf", &b[i]);
for (int i = 0; i < dim; i++)
{
for (int j = 0; j < dim; j++) cout << A[i][j] << " ";
cout << endl;
}
cout << endl;
cout << " B: " << endl;
for (int i = 0; i < dim; i++) cout << b[i] << endl;
//-------------SEIDEL-----------------////
int iter_seidel = 0;
double err_s = 0;
double *x_cur_s = new double[dim];
double *x_prev_s = new double[dim];
for (int i = 0; i < dim; i++) x_cur_s[i] = b[i]; //начальное приближение
do
{
for (int i = 0; i < dim; i++)
x_prev_s[i] = x_cur_s[i];
for (int i = 0; i < dim; i++)
{
double var = 0;
for (int j = 0; j < i; j++)
var += (A[i][j] * x_cur_s[j]);
for (int j = i + 1; j < dim; j++)
var += (A[i][j] * x_prev_s[j]);
x_cur_s[i] = (b[i] - var) / A[i][i];
}
err_s = residual(A, x_cur_s, b, dim);
iter_seidel++;
} while (err_s > eps);
cout << endl;
cout << " SEIDEL method " << endl;
cout << endl;
cout << " Result: " << endl;
for (int i = 0; i < dim; i++)
{
cout << x_cur_s[i] << " ";
}
cout << endl;
cout << " iter = " << iter_seidel << endl;
//-------------RICHARDSON-----------------//
cout << endl;
cout << " RICHARDSON method " << endl;
//Вычисление L_max
cout << endl;
double lambda_max;
lambda_max = eigen_val(A, dim);
cout << " LAMBDA_max = " << lambda_max << endl;
double ** B = new double *[dim];
for (int i = 0; i < dim; i++)
B[i] = new double[dim];
for (int i = 0; i < dim; i++)
{
for (int j = 0; j < dim; j++) B[i][j] = -A[i][j];
B[i][i] += lambda_max;
}
double lambda_0, lambda_min;
lambda_0 = eigen_val(B, dim);
lambda_min = lambda_max - lambda_0;
cout << " LAMBDA_min = " << lambda_min << endl;
double tau = 2 / (lambda_max + lambda_min);
int iter_richard = 0;
double err_r = 0;
double *x_cur_r = new double[dim];
double *x_prev_r = new double[dim];
for (int i = 0; i < dim; i++) x_cur_r[i] = b[i]; //начальное приближение
do
{
for (int i = 0; i < dim; i++) x_prev_r[i] = x_cur_r[i];
for (int i = 0; i < dim; i++)
{
double var = 0;
for (int j = 0; j < dim; j++)
var += (A[i][j] * x_prev_r[j]); //A*x(k)
x_cur_r[i] = x_prev_r[i] + tau*(b[i] - var); //x(k+1) = x(k)+tau*(b - A*x(k))
}
err_r = residual(A, x_cur_r, b, dim);
iter_richard++;
} while (err_r > eps);
cout << endl;
cout << " Result: " << endl;
for (int i = 0; i < dim; i++)
{
cout << x_cur_r[i] << " ";
}
cout << endl;
cout << " iter = " << iter_richard << endl;
cout << endl;
system("pause");
return 0;
}
Be the first to comment
You can use [html][/html], [css][/css], [php][/php] and more to embed the code. Urls are automatically hyperlinked. Line breaks and paragraphs are automatically generated.