VPF::Решение СЛАУ методом минимальных невязок

DUKE009

Дата 12.5.2012, 19:23 (ссылка)

(нет голосов)

Загрузка ...

Новичок

Профиль
Группа: Участник
Сообщений: 3
Регистрация: 12.5.2012

Репутация: нет
Всего: нет

Всем привет!
Помогите разобраться с данным методом, алгоритм я читал здесь:
http://www.physchem.chimfak.rsu.ru/S...ys_lin_eq.html

уже успел реализовать
функцию перемножения матрицы на вектор:

Код


void MatrVekt(int N, double **M, double *V, double *R)
//N- размерность, M- матрица, V- вектор, R- результат
{
for(int i=0; i<N; i++)
        {
        R[i]=0;
        for(int j=0; j<N; j++)
              R[i]+= M[i][j]*V[j];
        }
}

и сам алгоритм:

Код


void MinNev(int N, double **A, double *F, double *X, float eps)
//N- размерность, A- матрица, F- вектор свободных членов, X-вектор результат
{
    int count=0;
double *R = new double [N];
double *Delta = new double [N];
double *TempX = new double[N];
double maxi=0.0, Tau=0.0, TempTau=0.0;
for (int i=0; i<N; i++)
    TempX[i]=0;//первое приближение задаём нулевым
do
{
MatrVekt(N, A, TempX, R);
for(int i=0; i<N; i++)
    {
    Delta[i]=F[i]-R[i];//Вектор невязок
    }
MatrVekt(N, A, Delta, R);
for(int i=0; i<N; i++)
    {
    Tau+=-(R[i]*Delta[i])/(R[i]*R[i]);
    } 
for(int i=0; i<N; i++)
    {
    X[i]=TempX[i]-Tau*Delta[i];
    } 
 
maxi = fabs(X[0] - TempX[0]);
for(int i=0; i<N; i++)
    {
    if(fabs(X[i]-TempX[i])>maxi)
        maxi=fabs(X[i]-TempX[i]);
    }
count++;
}
while (maxi>=eps);
 printf("%d",count);
delete[] R;
delete[] Delta;
delete[] TempX;
}

Если не сложно, помогите найти ошибку. Суть в том, что при запуске функция уходит в бесконечный цикл. Условие выхода вроде реализовано верно. Короче не знаю в чём дело..

W4FhLF

Дата 12.5.2012, 19:31 (ссылка)

(нет голосов)

Загрузка ...

found myself

Профиль
Группа: Участник Клуба
Сообщений: 2831
Регистрация: 2.12.2006

Репутация: 2
Всего: 121

Не знаю насчёт правильности алгоритма, но с чего вы взяли, что он уходит в бесконечный цикл? Метод может просто не сходится на тех матрице и векторе, что вы пытаетесь решить.

--------------------

"Бог умер" © Ницше
"Ницше умер" © Бог

DUKE009

Дата 12.5.2012, 22:06 (ссылка)

(нет голосов)

Загрузка ...

Новичок

Профиль
Группа: Участник
Сообщений: 3
Регистрация: 12.5.2012

Репутация: нет
Всего: нет

W4FhLF, на матрице на которой проверяю, метод сходится(у неё есть диагональное преобладание)
Во всём виноваты мои кривые руки smile

Все ошибки алгоритма исправил, вроде считает верно.
В нете никак не мог нарыть готовой реализации метода, может кому пригодится smile

Код


int MinNev(int N, double **A, double *F, double *X, double eps)
//N- размерность, A- матрица, F- вектор свободных членов, X-вектор результат eps- точность
{
int count=0;//  количество итераций
double *R = new double [N];
double *Delta = new double [N];
double *TempX = new double[N];
double maxi=0.0, Tau=0.0, TempTau=0.0;
for (int i=0; i<N; i++)
    TempX[i]=0;//первое приближение задаём нулевым
do
{
MatrVekt(N, A, TempX, R);
for(int i=0; i<N; i++)
    {
    Delta[i]=R[i]-F[i];//Вектор невязок
    }

MatrVekt(N, A, Delta, R);

Tau=0.0; 
TempTau=0.0;

for(int i=0; i<N; i++)
    {
    Tau+=R[i]*Delta[i];
    TempTau+=R[i]*R[i];
    } 
Tau=Tau/TempTau; 

for(int i=0; i<N; i++)
    X[i]=TempX[i]-Tau*Delta[i];

maxi = fabs(X[0] - TempX[0]);
for(int i=0; i<N; i++)
    {
    if(fabs(X[i]-TempX[i])>maxi)
        maxi=fabs(X[i]-TempX[i]);
    TempX[i]=X[i];
    }
count++;
}
while (maxi>=eps);
delete[] R;
delete[] Delta;
delete[] TempX;
return count;
}

teesync

Дата 27.7.2022, 00:07 (ссылка)

(нет голосов)

Загрузка ...

Новичок

Профиль
Группа: Участник
Сообщений: 36
Регистрация: 26.7.2022

Репутация: нет
Всего: нет

Модератор: Сообщение скрыто.

squilia

Дата 15.8.2022, 06:14 (ссылка)

(нет голосов)

Загрузка ...

Новичок

Профиль
Группа: Участник
Сообщений: 13
Регистрация: 15.8.2022

Репутация: нет
Всего: нет

Модератор: Сообщение скрыто.

1 Пользователей читают эту тему (1 Гостей и 0 Скрытых Пользователей)
0 Пользователей:
« Предыдущая тема \| C/C++: Для новичков \| Следующая тема »