Решение СЛАУ методом Якоби

@klimkin3teers · Регистрация: 23.03.2016

Author24 — интернет-сервис помощи студентам

Пытаюсь реализовать метод решения системы линейных уравнений методом Якоби.

Кликните здесь для просмотра всего текста

C++

#include "stdafx.h"
#include "conio.h"
#include "iostream"
using namespace std;
double eps = 1;
const int n=3;
double a[n][n], c[n], xx[n];
double E[n][n], D[n][n];
void enterA() , enterX0(), enterC();
void inversion(double **A, int N);
void out();void returnB();
double **D1 = new double *[n];
double returnDA(int I,int J);
double B[n][n];double returnG(int I);
double g[n];
double returnBx(int I);
double xXx[n];
 
int main()
{
    double norm;
    setlocale(LC_ALL, "Russian");
    for (int i = 0; i < n; i++)
    D1[i] = new double[n];
    enterA();
    enterC();
    enterX0();
    out();
    returnB();
    
 
    
    for (int i = 0;i < n;i++)
    {
        xXx[i] = xx[i];
    }
    double Bx[n];
    do
    {
        for (int i = 0;i < n;i++)
        {
            Bx[i] = returnBx(i);
        }
        for (int i = 0;i < n;i++)
        {
            xXx[i] = Bx[i] - g[i];
        }
        norm = xXx[0] - xx[0];
        for (int i = 0;i < n;i++)
        {
            xx[i] = xXx[i];
        }
    } while (norm > eps);
    cout << "Корни:" << endl;
    for (int i = 0;i < n;i++)
    {
        cout <<xXx[i]<<" "<<endl;
    }
        
    getch();
 
}
 
void enterA()
{
 
    cout << "Введите коэффициенты: "<<endl;
    for (int i = 0;i <= (n - 1);i++)
    {
        for (int j = 0;j <= (n - 1);j++)
        {
            cout << "Введите a[" << i << "][" << j << "]: ";
            cin >> a[i][j];
            cout << endl;
        }   
    }
}
 
void enterC()
{
    cout << "Введите свободные члены: " << endl;
    for (int i = 0;i <= (n - 1);i++)
    {
        cout << "Введите c[" << i << "]: ";
        cin >> c[i];
    }cout << endl;
}
 
void enterX0()
{
    cout << "Введите начальные приближения: " << endl;
    for (int i = 0;i <= (n - 1);i++)
    {
        cout << "Введите хх[" << i << "]: ";
        cin >> xx[i];
    }cout << endl;
}
 
void out()
{
    cout << "Введенные коэффициенты представляют следующую матрицу:" << endl;
    for (int i = 0;i <= (n - 1);i++)
    {
        for (int j = 0;j <= (n - 1);j++)
        {
            cout << a[i][j] << " ";
        }
        cout << endl;
    }
 
    cout << "Столбец свободных членов:" << endl;
    for (int i = 0;i < n;i++)
    {
        cout << c[i] << endl;
    }
 
    cout << "Столбец начальных приближений:" << endl;
    for (int i = 0;i < n;i++)
    {
        cout << xx[i] << endl;
    }
 
 
    for (int i = 0;i < n;i++)
    {
        for (int j = 0;j < n;j++)
        {
            E[i][j] = 1;
        }
    }
    for (int i = 0;i < n;i++)
    {
        for (int j = 0;j < n;j++)
        {
            if (j == i) D[i][j] = a[i][j];
            else D[i][j] = 0;
        }
    }
 
    cout << "Матрица единиц:" << endl;
 
    for (int i = 0;i < n;i++)
    {
        for (int j = 0;j < n;j++)
        {
            cout << E[i][j] << " ";
        }cout << endl;
    }
 
    cout << "Матрица D:" << endl;
 
    for (int i = 0;i < n;i++)
    {
        for (int j = 0;j < n;j++)
        {
            cout << D[i][j] << " ";
        }cout << endl;
    }
 
    
 
    for (int i = 0; i < n; i++)
        for (int j = 0; j < n; j++)
        {
            D1[i][j] = D[i][j];
        }
 
    inversion(D1, n);
 
    cout << "Матрица D^-1:" << endl;
 
    for (int i = 0;i < n;i++)
    {
        for (int j = 0;j < n;j++)
        {
            cout << D1[i][j] << " ";
        }cout << endl;
    }
}
 
void inversion(double **A, int N)
{
    double temp;
 
    double **E = new double *[N];
 
    for (int i = 0; i < N; i++)
        E[i] = new double[N];
 
    for (int i = 0; i < N; i++)
        for (int j = 0; j < N; j++)
        {
            E[i][j] = 0.0;
 
            if (i == j)
                E[i][j] = 1.0;
        }
 
    for (int k = 0; k < N; k++)
    {
        temp = A[k][k];
 
        for (int j = 0; j < N; j++)
        {
            A[k][j] /= temp;
            E[k][j] /= temp;
        }
 
        for (int i = k + 1; i < N; i++)
        {
            temp = A[i][k];
 
            for (int j = 0; j < N; j++)
            {
                A[i][j] -= A[k][j] * temp;
                E[i][j] -= E[k][j] * temp;
            }
        }
    }
 
    for (int k = N - 1; k > 0; k--)
    {
        for (int i = k - 1; i >= 0; i--)
        {
            temp = A[i][k];
 
            for (int j = 0; j < N; j++)
            {
                A[i][j] -= A[k][j] * temp;
                E[i][j] -= E[k][j] * temp;
            }
        }
    }
 
    for (int i = 0; i < N; i++)
        for (int j = 0; j < N; j++)
            A[i][j] = E[i][j];
 
    for (int i = 0; i < N; i++)
        delete[] E[i];
 
    delete[] E;
}
 
void returnB()
{
    double D1A[n][n];
        for (int i = 0;i < n;i++)
        {
            for (int j = 0;j < n;j++)
            {
                D1A[i][j] = returnDA(i,j);
            }
        }
        cout << "D1A:" << endl;
        for (int i = 0;i < n;i++)
        {
            for (int j = 0;j < n;j++)
            {
                cout << D1A[i][j] << " ";
            }cout << endl;
        }
        
        for (int i = 0;i < n;i++)
        {
            for (int j = 0;j < n;j++)
            {
                B[i][j] = E[i][j] - D1A[i][j];
            }
        }
        cout << "Матрица B:" << endl;
        for (int i = 0;i < n;i++)
        {
            for (int j = 0;j < n;j++)
            {
                cout << B[i][j] << " ";
            }cout << endl;
        }
 
        for (int i = 0;i < n;i++)
        {
            g[i] = returnG(i);
        }
        cout << "Матрица G" << endl;
        for (int i = 0;i < n;i++)
        {
        
                cout << g[i] << " ";
            
        }
 
}
 
double returnDA(int I,int J)
{
    double da=0;
    
    for (int i = 0;i < n;i++)
    {
        da += D1[I][i] * a[i][J];
    }
 
    return da;
}
 
double returnG(int I)
{
    double G=0;
 
    for (int i = 0;i < n;i++)
    {
        G += D1[I][i] * c[i];
    }
    return G;
}
 
double returnBx(int I)
{
    double Bx = 0;
 
    for (int i = 0;i < n;i++)
    {
        Bx += B[I][i] * xx[i];
    }
    return Bx;
}

Может кому-то этот код не будет резать глаза, и получится посмотреть что делаю не так. Все делал по порядку на основе формул из википедии, последовательно получал все необходимые матрицы-множители, все вроде правильно находилось. И вот я сделал на основе всех полученных множителей подсчет корней, и полученные как-то совсем не сходятся с проверкой по онлайн калькулятору СЛАУ.

@Antikl · 05.06.2017, 09:24

C++

#include <math.h>
const double eps = 0.001; ///< желаемая точность 
 
..........................
 
/// N - размерность матрицы; A[N][N] - матрица коэффициентов, F[N] - столбец свободных членов,
/// X[N] - начальное приближение, ответ записывается также в X[N];
void Jacobi (int N, double** A, double* F, double* X)
{
    double* TempX = new double[N];
    double norm; // норма, определяемая как наибольшая разность компонент столбца иксов соседних итераций.
 
    do {
        for (int i = 0; i < N; i++) {
            TempX[i] = F[i];
            for (int g = 0; g < N; g++) {
                if (i != g)
                    TempX[i] -= A[i][g] * X[g];
            }
            TempX[i] /= A[i][i];
        }
        norm = fabs(X[0] - TempX[0]);
        for (int h = 0; h < N; h++) {
            if (fabs(X[h] - TempX[h]) > norm)
                norm = fabs(X[h] - TempX[h]);
            X[h] = TempX[h];
        }
    } while (norm > eps);
    delete[] TempX;
}

@klimkin3teers · 05.06.2017, 09:29 **[ТС]**

я пробовал брать этот код, но как-то не получается им пользоваться, я немного не догоняю объявление массивов через new double*[N]

@Antikl · 05.06.2017, 09:39

C++

const unsigned int DIM1 = 3;
const unsigned int DIM2 = 5;
 
//создание двумерного динамического масива
int **ary= new int * [DIM1];
 
 for (int i = 0; i < DIM1; i++) { 
        ary[i] = new int [DIM2]; 
    }
 
// работа с массивом
    for (int i = 0; i < DIM1; i++) {
        for (int j = 0; j < DIM2; j++) {
            ary[i][j] =  i +10;
        }
    }
 
вывод
for (int i = 0; i < DIM1; i++) {
        for (int j = 0; j < DIM2; j++) {
            cout << setw(4) << ary[i][j];
        }
        cout << endl;
    }
 
// уничтожение
    for (int i = 0; i < DIM1; i++) {
        delete [] ary[i];
    }
    delete [] ary;

Добавлено через 2 минуты
так же можно с double и другими числовыми типами

@klimkin3teers · 05.06.2017, 12:12 **[ТС]**

спасибо, дошло)

Добавлено через 16 секунд
а вот все таки не правильно работает код из вкипедии, ну видимо я что-то не так делаю, while зацикливается, а все корни в каждом цикле неправильные, мой изначальный код тоже зацикливался, но корни находил несколько больше похожие на действительные.

Кликните здесь для просмотра всего текста

C++

#include "stdafx.h"
#include <math.h>
#include <iomanip>
#include "iostream"
using namespace std;
const double eps = 0.1; ///< желаемая точность 
void Jacobi(int N, double** A, double* F, double* X);
 
int main()
{
    setlocale(LC_ALL, "Russian");
    const unsigned int DIM1 = 3;
    const unsigned int DIM2 = 3;
 
    //создание двумерного динамического масива
    double **a = new double *[DIM1];
 
    for (int i = 0; i < DIM1; i++) {
        a[i] = new double[DIM2];
    }
 
    double * c = new double[DIM1];
    double * x = new double[DIM1];
 
    // работа с массивом
    for (int i = 0; i < DIM1; i++) {
        for (int j = 0; j < DIM2; j++) {
            cout << "Введите a[" << i << "][" << j << "]: ";
            cin >> a[i][j];cout << endl;
        }
    }
 
    for (int j = 0; j < DIM2; j++) {
        cout << "Введите c[" << j << "][" << j << "]: ";
        cin >> c[j];cout << endl;
    }
 
    for (int j = 0; j < DIM2; j++) {
        cout << "Введите x[" << j << "][" << j << "]: ";
        cin >> x[j];cout << endl;
    }
    Jacobi(DIM1, a, c, x);
    //вывод
    //for (int i = 0; i < DIM1; i++) {
    for (int j = 0; j < DIM2; j++) {
        cout << setw(4) << x[j];
    }
    cout << endl;
    //}
    system("pause");
 
    // уничтожение
    for (int i = 0; i < DIM1; i++) {
        delete[] a[i];
    }
    delete[] a;
    return 0;
}
 
/// N - размерность матрицы; A[N][N] - матрица коэффициентов, F[N] - столбец свободных членов,
/// X[N] - начальное приближение, ответ записывается также в X[N];
void Jacobi(int N, double** A, double* F, double* X)
{
    double* TempX = new double[N];
    double norm; // норма, определяемая как наибольшая разность компонент столбца иксов соседних итераций.
 
    do {
        for (int i = 0; i < N; i++) {
            TempX[i] = F[i];
            for (int g = 0; g < N; g++) {
                if (i != g)
                    TempX[i] -= A[i][g] * X[g];
            }
            TempX[i] /= A[i][i];
        }
        norm = fabs(X[0] - TempX[0]);
        for (int h = 0; h < N; h++) {
            if (fabs(X[h] - TempX[h]) > norm)
                norm = fabs(X[h] - TempX[h]);
            X[h] = TempX[h];
        }for (int i = 0;i < 3;i++)
    {
        cout << X[i]<<" ";
        }cout << endl;
    } while (norm > eps);
 
    
 
    delete[] TempX;
}

@Antikl · 05.06.2017, 12:34

пример в интернете нашол может он подойдет

C++

#include <iostream>
#include <cmath>
using namespace std;
// приведение матрицы коэффициентов к виду с ненулевой диагональю и соответствующее изменение вектора правых частей
// возвращает true - если приведение - успешно
bool getDiagonal( double **coefficients, double *rightPart, int currColumn, int numberOfEquation ) {
    bool result = false;
    int i, row = currColumn;
    double tempItem;
 
    // для матрицы 1x1 ответом является ненулевость ее единственного элемента
    if ( currColumn == numberOfEquation - 1 ) {
        result = coefficients[currColumn][currColumn] != 0;
    }
    else {
        // пока не найдена перестановка приводящая матрицу к виду с ненулевой диагональю и пока мы можем еще просматривать строки ищем перестановку
        while ( !result && row < numberOfEquation ) {
            // если текущий элемент на диагонали нулевой ищем в столбце дальше ненулевой
            if ( coefficients[row][currColumn] != 0 ) {
                // если элемент ненулевой и не лежит на диагонали меняем местами сотвествующие строки в матрице и элементы в векторе прваых частей
                if ( row != currColumn ) {
                    for ( i = 0; i < numberOfEquation; i++ ) {
                        tempItem = coefficients[currColumn][i];
                        coefficients[currColumn][i] = coefficients[row][i];
                        coefficients[row][i] = tempItem;
                    }               
                    tempItem = rightPart[currColumn];
                    rightPart[currColumn] = rightPart[row];
                    rightPart[row] = tempItem;                  
                }
                // рекурсивный вызов фактически для подматрицы с размерностью на 1 меньше
                result = getDiagonal( coefficients, rightPart, currColumn + 1, numberOfEquation );
                if ( result ) {
                    break;
                }
            }
            row ++;
        }
    }
 
    return result;
}
 
// было ли найдено решение, если да - итог в параметре solution
int iteration( double **coefficients, double *rightPart, int numberOfEquation, double *solution, double precision ) {
    bool result;
    int i, j, k, step = 1;
    double temp;
    double* tempSolution;
 
    tempSolution = new double[numberOfEquation];
    
    // приведение матрицы коэффициентов к виду с ненулевой диагональю и соответствующее изменение вектора правых частей
    result = getDiagonal( coefficients, rightPart, 0, numberOfEquation );   
 
    // если приведение успешно - работаем дальше
    if ( result ) {
        double fault = precision + 1;
 
        // преобразуем матрицу и правую часть для дальнейшего решения
        for ( i = 0; i < numberOfEquation; i ++ ) {
            for ( j = 0; j < numberOfEquation; j ++ ) {
                if ( i != j ) {
                    coefficients[i][j] = - coefficients[i][j] / coefficients[i][i];
                }
            }
            rightPart[i] = rightPart[i] / coefficients[i][i];
            coefficients[i][i] = 0;
        }
 
        // первое приближение решения - преобразованный вектор правых частей
        for ( i = 0; i < numberOfEquation; i ++ ) {
            solution[i] = rightPart[i];
        }
 
        // пока не найдено решение с допустимй погрешнстью или пока не исчерпан лимит шагов... если все расходится например
        while ( fault > precision && step <= 1000 ) {
 
            // поиск новой итерации с ее "самоиспользованием" при расчетах            
            for ( j = 0; j < numberOfEquation; j ++ ) {
                tempSolution[j] = 0;
            }
            for ( i = 0; i < numberOfEquation; i ++ ) {
                for ( j = 0; j < numberOfEquation; j ++ ) {
                    tempSolution[i] = tempSolution[i] + coefficients[i][j] * solution[j]; 
                }
                tempSolution[i] = tempSolution[i] + rightPart[i];
            }
 
            // расчет погрешности
            fault = 0.0;
            for ( j = 0; j < numberOfEquation; j ++ ) {
                fault = fault + (solution[j] - tempSolution[j])*(solution[j] - tempSolution[j]);
            }
            fault = sqrt( fault );
 
            // сохранение полученной новой итерации
            for ( j = 0; j < numberOfEquation; j ++ ) {
                solution[j] = tempSolution[j];
            }
            step++;
        }
    }
    else {
        result = 1001;
    }
    
 
    return step;
}
 
 
int main() {
    int i, j;
    int size;
    double **coefficients, *rightPart, *solution, precision;
 
    cout << "Simple iteration method.\nEnter system dimension: ";
    cin >> size;
 
    coefficients = new double*[size];
    for ( i = 0; i < size; i++ ) {
        coefficients[i] = new double[size];
    }
    rightPart = new double[size];
    solution = new double[size];
 
    for ( i = 0; i < size; i ++ ){
        cout << "Enter " << i + 1 << " row: ";
        for ( j = 0; j < size; j ++ ){
            cin >> coefficients[i][j];
        }
    }
 
    cout << "Enter right part: ";
    for ( j = 0; j < size; j ++ ){
        cin >> rightPart[j];
    }
 
    cout << "Enter precision: ";
        cin >> precision;
 
    int steps = iteration( coefficients, rightPart, size, solution, precision );
    if ( steps > 1000 ) {
        cout << "Solution for this matrix of coefficients not exist or not found";
    }
    else {
        cout << "Solution is:\n";
        for ( j = 0; j < size; j ++ ){
            cout << solution[j] << "\n";
        }
        cout << "Number of step: " << steps;
    }
 
    cout << "\nPress "Enter" to continue..." << endl; 
return 0;
}

@klimkin3teers 0 / 0 / 0 Регистрация: 23.03.2016 Сообщений: 15
	05.06.2017, 09:29 [ТС]	3
	я пробовал брать этот код, но как-то не получается им пользоваться, я немного не догоняю объявление массивов через new double*[N] 0

	05.06.2017, 12:34