2.3. Программы для решения систем линейных алгебраических уравнений

Файл GAUSS.CPP.

// Подключение математики:

#include <math.h>

// Математика подключена.

gauss(a,b,x,n)

float a[][nmax]; // Исходная матрица коэффицентов :a

float b[]; // Правые части

float x[]; // Результат (будет получен при решении)

int n; // Количество элементов в векторах и размерность матрицы

{

int i,j,k,imax;

int nm1=n-1; // Используется для удобства вместо n, так как

// элементы массивов отсчитываются с 0, а не с 1

int err=0; // Начальное значение кода ошибки =0 (т.е. ошибки нет),

// перед окончанием функции переменная err передается

// в качестве значения функции

float amax,t,mi,sum;

// Проверка корректности переменной n и константы nmax

// в случае противоречий выход из функции:

// endofproc -метка конца функции

if (n<2) { err=2; goto endofproc; } // В матрице <2 элементов

if (n>nmax) { err=3; goto endofproc; }

// n>реальной размерности =nmax

// Теперь все корректно.

// Далее следует ...

// Цикл из n-1 шагов, на каждом шаге находится главный элемент и

// преобразовывается матрица.

// Главный элемент определяется только из очередного столбца,

// а не из всей матрицы.

// n-й шаг цикла не нужен так как остается один элемент a[n-1][n-1],

// который сам по себе главный (рассматривается после цикла).

for (k=0; k<=nm1-1; k++)

{

// Поиск главного элемента в текущем столбце (k):

for (i=k,amax=0.0; i<=nm1; i++)

if ( fabs(a[i][k]) > fabs(amax) ) amax=a[(imax=i)][k];

// Главный элемент найден.

// Если он =0, то ошибка, так как матрица вырождена и выход из функции

// (метка endofproc стоит в конце функции) :

if (amax==0.0) { err=1; goto endofproc; }

// Замена строки с найденным главным элементом на строку с номером k

for (j=k;j<=nm1;j++)

{

t=a[imax][j];

a[imax][j]=a[k][j];

a[k][j]=t;

}

// Строки заменены.

// Замена элементов правой части для тех же строк:

t=b[imax];

b[imax]=b[k];

b[k]=t;

// Заменены.

// Определение mi для каждой строки (ниже главной)

// и ее преобразование для найденного mi:

for (i=k+1;i<=nm1;i++)

{

mi=-a[i][k]/amax;

a[i][k]=0.0;

for (j=k+1;j<=nm1;j++)

a[i][j]+=a[k][j]*mi;

b[i]+=b[k]*mi;

}

// На текущем шаге k матрица и правые части преобразованы.

} // Конец основного цикла (k увеличивается на единицу и все с начала)

// Проверка единственного оставшегося (и поэтому главного) элемента

// если он=0, то ошибка, так как матрица вырождена и выход из функции

if (a[nm1][nm1]==0.0) { err=1; goto endofproc; }

// Теперь матрица a и вектор b преобразованы (подготовлены к поиску x)

// Решение полученной системы (обратный ход):

for (i=nm1;i>=0;i--)

{

sum=0.0;

for (j=i+1;j<=nm1;j++)

sum+=a[i][j]*x[j];

X[i]=(b[i]-sum)/a[i][i];

}

// Решение получено.

endofproc: // Конец процедуры: метка для аварийного выхода

return(err); // Функция будет иметь значение ошибки err

}

Файл MITER.CPP.

// Подключение математики:

#include <math.h>

// Математика подключена.

// Функция получения мантиссы и порядка float числа:

mantpor(x,mant,por)

float *x,*mant;

double *por;

{

int st;

st=(int)floor(log10(*x));

*por=pow10(st);

*mant=*x/(*por);

}

// Конец функции mantpor.

// Функция, округляющая float число до n разрядов:

float round(x,n)

float *x;

int n;

{

float xabs,mant,dpown,xr;

double por;

int signx;

if (*x!=0.0)

{

xabs=fabs(*x);

signx=(*x>=0)? 1:-1;

mantpor(&xabs,&mant,&por);

dpown=pow10(n);

xr=signx*(floor(mant*dpown+0.5)/dpown)*por;

}

else xr=0.0;

return (xr);

}

// Конец функции round.

// Функция, решающая приведенную систему методом итераций:

miter(a,b,x,n,ntch,it)

float a[][nmax]; // Приведенная матрица коэффициентов.

float b[]; // Вектор свободных членов.

float x[]; // Вектор неизвестных.

int n, // Количество элементов в векторах и размерность матрицы.

ntch, // Количество знаков после запятой в мантиссе результата,

// которые должны остаться после округления.

*it; // Количество итераций, потребовавшихся для получения x.

{

float sum,xst[nmax];

int err=0; // Значение ошибки (сейчас=0, т.е. ее нет)

// в конце функции в качестве ее значения возвращается

// эта переменная.

int nm1=n-1; // Используется вместо n, так как

// массивы отсчитываются с 0, а не с 1

int ntchp1=ntch+1;

int i,j,tochnok;

// Проверка корректности переменной n и константы nmax

// в случае противоречий ошибка и выход из функции:

// endofproc - метка в конце функции

if (n<2) {err=2; goto endofproc;} // В матрице <2 элементов

if (n>nmax) {err=3; goto endofproc;} // n>реальной размерности =nmax

// Теперь все корректно.

// Проверка на корректность количества разрядов мантиссы для округления:

// Если n выходит за ограничения, то ошибка и выход из функции:

if ( ! ((0<=ntch)&&(ntch<=6)) ) {err=1; goto endofproc;}

// Проверка закончена.

// Проверка соответствия матрицы a условиям сходимости процесса:

// По строкам:

for (i=0,sum=0.0; (i<=nm1)&&(sum<1.0); i++)

for (j=0,sum=0.0; j<=nm1; j++) sum+=fabs(a[i][j]);

// По строкам проверено.

// Если по строкам не подходит, то по столбцам:

if (sum>=1)

{

// По столбцам:

for (j=0,sum=0.0; (j<=nm1)&&(sum<1.0); j++)

for (i=0,sum=0.0; i<=nm1 ;i++) sum+=fabs(a[i][j]);

// По столбцам проверено.

// Если и по столбцам не подходит, то ошибка и выход из функции:

if (sum>=1) {err=1; goto endofproc;}

}

// Проверка на сходимость закончена.

// Начальная установка массива x (в принципе может быть любая):

for (j=0;j<=nm1;j++) x[j]=b[j];

// Закончена.

*it=0; // Обнуляется счетчик итераций.

// Основной цикл: на каждом шаге находятся новые значения x,

// выход при достижении необходимой точности x:

do

{

(*it)++; // Увеличение счетчика итераций.

for (j=0;j<=nm1;j++) xst[j]=x[j]; // Запоминает старое x в xst

// Получает новые x:

for (i=0;i<=nm1;i++)

{

for (j=0,sum=0.0; j<=nm1; j++) sum+=a[i][j]*xst[j];

X[i]=b[i]+sum;

}

// Новые значения x получены.

// Проверяет отличие нового x от старого,

// т. е. точность до ntch+1 знаков мантиссы:

for (j=0,tochnok=1; (j<=nm1)&&tochnok; j++)

tochnok=( round(&x[j],ntchp1)==round(&xst[j],ntchp1) );

// Точность проверена, если все нормально, то tochnok=1 иначе 0

} // Конец тела основного цикла.

while ( !(tochnok) ); // Выход, если достигнута необходимая точность.

// Округление найденного x до ntch знаков мантиссы:

for (j=0;j<=nm1;j++) x[j]=round(&x[j],ntch);

// Округлены.

endofproc: // Метка для аварийного выхода в случае ошибок:

return (err); // Код ошибки возвращается в качестве значения функции.

}