2.2. Ðешение систем линейных алгебраических уравнений методом простой итерации
Рассматривается система уравнений вида
,
(2.3)
ãäå
- заданная числовая квадратная матрицаn-го
порядка, аb- заданный вектор (свободный член). Метод
простой итерации состоит в следующем.
Выбирается произвольный векторx
(начальное приближение),и строится итерационная последовательность
векторов по формуле
,
(2.4)
где k =1, 2, …
Доказана
теорема, что если норма
,
то система уравнений (2.3) имеет единственное
решение
и итерации (2.4) сходятся к решению со
скоростью геометрической прогрессии[13].
Для оценки погрешностиk-го
приближения широко применяется
неравенство
,
(2.5)
которое может быть использовано для принятия решения об останове итерационного процесса при выполнении условия
,
ãäå
- некоторая заданная погрешность
вычислений.
Лабораторная работа № 12
В
работе студенты должны найти решение
системы линейных уравнений с nнеизвестными, заданной матрицей
коэффициентов
и вектором свободных членовb,
методом простых итераций. Выполнение
работы включает следующие этапы:
1) с помощью преподавателя определить систему уравнений, которую нужно решить. Привести исходную систему к виду (2.3), пригодному для использования метода простых итераций;
2) задать необходимую точность получения результата (количество знаков мантиссы числа);
3) разработать программу решения задачи на языке С с использованием подпрограммы-функции MITERиз файлаMITER.CPP. ФункцияMITERимеет следующие параметры:
-
- матрица коэффициентов преобразованной
к виду (2.3) системы уравнений размера
,
òèï
;
-
- вектор свободных членов преобразованной
системы размера
,
òèï
;
-
- полученный в результате
проведения итераций вектор решения
размера
,
òèï
;
-
- размер системы уравнений, тип
;
-
- количество знаков после запятой
в мантиссе результата, остающихся после
округления, тип
,
;
- it - выходной параметр, равный количеству произведенных итераций, тип int.
В
качестве значения функции типа
возвращается одно из следующих значений:
а)
0 - все нормально, получено решение
;
б) 1 - не выполняются условия сходимости итерационного процесса;
в)
2 - размер
;
г)
3 - значение
.
Константа
должна быть задана при разработке
головной программы аналогично тому,
как это делается при выполнении
лабораторной работы ¹ 11;
4) произвести вычисления с использованием разработанной программы и построить график зависимости числа итераций от задаваемой точности.
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); // Код ошибки возвращается в качестве значения функции.
}