vychmat kursach

Метод Гаусса
Классический метод решения систем линейных алгебраических уравнений (СЛАУ). Это метод последовательного исключения переменных, когда с помощью элементарных преобразований система уравнений приводится к равносильной системе ступенчатого (или треугольного) вида, из которого последовательно, начиная с последних (по номеру) переменных, находятся все остальные переменные.
История
Хотя в настоящее время данный метод повсеместно называется методом Гаусса, он был известен и до К.Ф. Гаусса. Первое известное описание данного метода в китайском трактате "Математика в девяти книгах", составленном между I в. До н.э. и II в. н.э.
Описание метода
Пусть исходная система выглядит следующим образом
Матрица A называется основной матрицей системы, b - столбцом свободных членов.
Тогда, согласно свойству элементарных преобразований над строками основную матрицу этой системы можно привести к ступенчатому виду (эти же преобразования нужно применять к столбцу свободных членов):
При этом будем считать, что базисный минор (ненулевой минор максимального порядка) основной матрицы находится в верхнем левом углу, то есть в него входят только коэффициенты при переменных .
Тогда переменные называются главными переменными. Все остальные называются свободным.
Если хотя бы одно число , где i > r, то рассматриваемая система несовместна.
Пусть для любых i > r.
Перенесем свободные переменные за знаки равенств и поделим каждое из уравнений системы на свой коэффициент при самом левом где i = 1, ...,r; k = i + 1, ..., n.
Если свободным переменным системы (2) придавать все возможные значения и решать новую систему относительно главных неизвестных снизу вверх (то есть от нижнего уравнения к верхнему), то мы получим все решения этой СЛАУ. Так как эта система получена путем элементарных преобразований над исходной системой (1), то по теореме об эквивалентности при элементарных преобразованиях системы (1) и (2) эквивалентны, то есть множества их решений совпадают.
Следствия
1. Если в совместной системе все переменные главные, то такая система является определенной.
2. Если количество переменных в системе превосходит число уравнений, то такая система является либо неопределенной, либо несовместной.
Условие совместности:
Упомянутое выше условие для всех i > r может быть сформулировано в качестве необходимого и достаточного условия совместности:
Напомним, что рангом совместной системы называется ранг её основной матрицы (либо расширенной, так как они равны).
Теорема Кронекера-Капелли:
Система совместна тогда и только тогда, когда ранг её основной матрицы равен рангу её расширенной матрицы.
Следствия:
Количество главных переменных равно рангу системы и не зависит от её решения.
Если ранг совместной системы равен числу переменных данной системы, то она определена.
Алгоритм
Описание:
Алгоритм решения СЛАУ методом Гаусса подразделяется на два этапа.
* На первом этапе осуществляется так называемый прямой ход, когда путём элементарных преобразований над строками систему приводят к ступенчатой или треугольной форме, либо устанавливают, что система несовместна. А именно: среди элементов первого столбца матрицы выбирают ненулевой, перемещают его на крайнее верхнее положение перестановкой строк и вычитают получившуюся после перестановки первую строку из остальных строк, домножив её на величину, равную отношению первого элемента каждой из этих строк к первому элементу первой строки, обнуляя тем самым столбец под ним. После того, как указанные преобразования были совершены, первую строку и первый столбец мысленно вычёркивают и продолжают пока не останется матрица нулевого размера. Если на какой-то из итераций среди элементов первого столбца не нашелся ненулевой, то переходят к следующему столбцу и проделывают аналогичную операцию.
* На втором этапе осуществляется так называемый обратный ход, суть которого заключается в том, чтобы выразить все получившиеся базисные переменные через небазисные и построить фундаментальную систему решений, либо, если все переменные являются базисными, то выразить в численном виде единственное решение системы линейных уравнений. Эта процедура начинается с последнего уравнения, из которого выражают соответствующую базисную переменную (а она там всего одна) и подставляют в предыдущие уравнения, и так далее, поднимаясь по "ступенькам" наверх. Каждой строчке соответствует ровно одна базисная переменная, поэтому на каждом шаге, кроме последнего (самого верхнего), ситуация в точности повторяет случай последней строки.
Метод Гаусса требует порядка действий.
Этот метод опирается на:
Теорема (о приведении матриц к ступенчатому виду)
Любую матрицу путем элементарных преобразований только над строками можно привести к ступенчатому виду.
Достоинства метода Гаусса:
- менее трудоёмкий по сравнению с другими методами;
- позволяет однозначно установить, совместна система или нет, и если совместна, найти её решение;
- позволяет найти максимальное число линейно независимых уравнений - ранг матрицы системы.
Метод Крамера (правило Крамера)
Способ решения квадратных систем линейных алгебраических уравнений с ненулевым определителем основной матрицы (причем для таких уравнений решение существует и единственно). Назван по имени Габриеля Крамера (1704-1752), придумавшего метод.
Описание метода
Для системы n линейных уравнений с n неизвестными (над произвольным полем)
с определителем матрицы системы Δ, отличным от нуля, решение записывается в виде
(i-ый столбец матрицы системы заменяется столбцом свободных членов).
В другой форме правило Крамера формулируется так: для любых коэффициентов c1, c2,..., cn:
В этой форме формуле Крамера справедлива без предположения, что Δ отлично от нуля, не нужно даже, чтобы коэффициенты системы были бы элементами целостного кольца (определитель системы может быть даже делителем нуля в кольце коэффициентов). Можно считать, что либо наборы b1, b2,..., bn и x1, x2,..., xn , либо набор с1, c2,..., cn состоят не из элементов кольца коэффициентов системы, а какого-нибудь модуля над этим кольцом. В этом виде формула Крамера используется, например, при доказательстве формулы для определителя Грамма и леммы Накаямы.
Вычислительная сложность:
Метод Крамера требует вычисления n + 1 определителей размерности n×n. При использовании метода Гаусса для вычисления определителей, метод имеет временную сложность порядка , что хуже, чем если бы метод Гаусса напрямую использовался для решения системы уравнений. Поэтому метод считался непрактичным. Однако в 2010 году было показано, что метод Крамера может быть реализован со сложностью , сравнимой со сложностью метода Гаусса.
Погрешности решения СЛАУ "точными" методами
Число обусловленности матрицы
Числом обусловленности матрицы А называется вещественное число, равное произведению нормы матрицы на норму обратной матрицы. Если Cond A~1 , то говорят, что матрица А хорошо определена. Если Cond A>>1 (значительно превосходит), то матрица плохо определена.
Классификация случаев СЛАУ:
1) Det ≠ 0, Det A = O(1), Cond A = O(1) матрица хорошо определена, имеем единственное решение СЛАУ. Здесь O(1) - символ Ландау;
2) Det A = 0, Cond A → решений СЛАУ не существует;
3) Det A = 0, Cond A → решений СЛАУ бесконечно много;
4) Det A = c1,|c1|<<1, Cond A >> 1. Этот случай самый распространённый на практике так, как определитель меленький и не равен нулю, то решение существует, но найти его сложно. На практике в этом случае имеем бесконечно много приближенных решений СЛАУ.
Виды погрешностей:
1. Неустранимые погрешности обуславливаются методиками измерения экспериментальных данных и округление при измерениях. Бессмысленно решат систему более точно, чем с этой ошибкой.
2. Вычислительная погрешность, которая связана с конечностью числа разрядов машинного слова.
3. Погрешность метода. Реально возмущаются и вектор правой части, и элементы матрицы. Однако для анализа надо предположить, что можно отделить эти два возмущения друг от друга.
Идеальная (невозмущенная) система Ax = b
Теорема 1
Пусть δА = 0. Относительная ошибка решения системы линейных уравнений.
A (x+δx) = b + δb, где функции δb и δx векторы возмущения решений и правой части имеет вид
Оценка справедлива в любой норме.
Теорема 2
Пусть δb = 0. Относительная ошибка решения системы линейных уравнений
(A + δA)(x + δx) = b, где функции δA и δx векторы возмущения решения и элементов матрицы имеет вид
#include "stdafx.h"
#include <math.h>
#include <iostream>
#include <conio.h>
#include <stdlib.h>
#include <complex>
using namespace std;
int w;
int z;
int NaddKr, NmultKr, NaddGauss, NmultGauss;
complex<double> D2,D3,D4,D5;
complex<double>dva(complex<double> a[2][2]);
complex<double>tri(complex<double>a[3][3]);
complex<double>chetire(complex<double>a[4][4]);
complex<double>pyat(complex<double>a[5][5]);
complex<double>normMatr(complex<double>a[5][5]);
complex<double>normVect(complex<double>b[5]);
void KRAMER();
void GAUSS5();
complex<double>Error(complex<double>a[5][5],complex<double>b[5]);
void menu();
complex<double>a[5][5];
complex<double> KR(complex<double>a[5][5], complex<double>x[5], int p);
int i,j,k,m,s,q,n,p;
//=======================================================================================
//функция нахождения определителя матрицы размерности 2
complex<double>dva(complex<double>a[2][2])
{
D2=a[0][0]*a[1][1]-a[0][1]*a[1][0];
NaddKr++;//подсчет аддитивных операций (+,-)
NmultKr=NmultKr+2;//подсчет мультипликативных операций (*,/)
return D2;
}
//функция нахождения определител матрицы размерности 3
complex<double>tri(complex<double>a[3][3])
{
complex<double>b[2][2];
D3=0;
for(s=0;s<3;s++)
{
for(j=0;j<2;j++)
{
m=0;
for(k=0;k<3;k++)
{
if(k!=s)
{
b[j][m]=a[j][k];
if (m<1) m++;
}
}
}
D3+=a[2][s]*dva(b)*pow((double)-1,s);
}
NaddKr++;
NmultKr=NmultKr+2;
return D3;
}
//функция нахождения определителя матрицы размерности 4
complex<double>chetire(complex<double>a[4][4])
{
complex<double>b[3][3];
D4=0;
for(i=0;i<4;i++)
{
for(j=0;j<3;j++)
{
m=0;
for(k=0;k<4;k++)
{
if(k!=i)
{
b[j][m]=a[j][k];
if(m<2) m++;
}
}
}
D4+=tri(b)*a[3][i]*pow((double)-1,(i+1));
NaddKr++;
NmultKr=NmultKr+2;
}
return D4;
}
//функция нахождения определителя матрицы размерности 5
complex<double> pyat(complex<double> a[5][5])
{
complex<double>b[4][4];
D5=0;
for(q=0;q<5;q++)
{
for(j=0;j<4;j++)
{
m=0;
for(k=0;k<5;k++)
{
if(k!=q)
{
b[j][m]=a[j][k];
if(m<3) m++;
}
}
}
D5+=chetire(b)*a[4][q]*pow((double)-1,q);
NaddKr++;
NmultKr=NmultKr+2;
}
return D5;
}
//===========================================================================================
complex<double>KR(complex<double> a[5][5],complex<double>x[5],int p)
{
complex<double>dp=0;
if(z==2) {complex<double>D[2]; complex<double>f[2][2]; for(i=0;i<z;i++) for(j=0;j<z;j++) f[i][j]=a[i][j]; for(n=0;n<z;n++) f[n][p]=x[n];
dp=dva(f);}
//копирование матрицы коэффициентов А в матрицу f и копирование матрицы значений х для матрицы размерности 2
if(z==3) {complex<double>D[3]; complex<double>f[3][3]; for (i=0;i<z;i++) for (j=0;j<z;j++) f[i][j]=a[i][j];for(n=0;n<z;n++) f[n][p]=x[n];
dp=tri(f);}
//копирование матрицы коэффициентов А в матрицу f и копирование матрицы значений х для матрицы размерности 3
if(z==4) {complex<double>D[4]; complex<double>f[4][4]; for (i=0;i<z;i++) for (j=0;j<z;j++) f[i][j]=a[i][j];for(n=0;n<z;n++) f[n][p]=x[n];
dp=chetire(f);}
//копирование матрицы коэффициентов А в матрицу f и копирование матрицы значений х для матрицы размерности 4
if(z==5) {complex<double>D[5]; complex<double>f[5][5]; for (i=0;i<z;i++) for (j=0;j<z;j++) f[i][j]=a[i][j];for(n=0;n<z;n++) f[n][p]=x[n];
dp=pyat(f);}
//копирование матрицы коэффициентов А в матрицу f и копирование матрицы значений х для матрицы размерности 5
return dp;
}
void KRAMER()
{
NaddKr=0;
NmultKr=0;
cout<<"\n Введите размерность от 2 до 5:\n";
cin>>z;
complex<double>a[5][5];
complex<double>y[5];
complex<double>x[5];
complex<double>D;
cout<<"\n Введите коэффициенты при x:\n";
for(i=0;i<5;i++) for(j=0;j<5;j++) a[i][j]=0;
for(i=0;i<z;i++)
for (j=0;j<z;j++) cin>>a[i][j]; //ввод матрицы коэффициентов при х
cout<<"\n Введите матрицу значений:\n";
for(i=0;i<z;i++) cin>>x[i]; //ввод матрицы значений
if(z==2) {complex<double>b[2][2]; for(i=0;i<z;i++) for(j=0;j<z;j++) b[i][j]=a[i][j]; D=dva(b);}//копирование матрицы A в матрицу b
if(z==3) {complex<double>b[3][3]; for(i=0;i<z;i++) for(j=0;j<z;j++) b[i][j]=a[i][j]; D=tri(b);}
if(z==4) {complex<double>b[4][4]; for(i=0;i<z;i++) for(j=0;j<z;j++) b[i][j]=a[i][j]; D=chetire(b);}
if(z==5) D=pyat(a);
if(D==(complex<double>)0) {cout<<"\n ERROR"; getch();} //если определитель равен 0, то ошибка
else for(p=0;p<z;p++)
{
y[p]=KR(a,x,p)/D; //вычисление корней уравнений
if(imag(y[p])!=0) cout<< " x["<<p<<"]="<<y[p];
else cout<< " x["<<p<<"]="<<real(y[p]);
NmultKr++;
}
cout<<"\n К-во мультипликативных операций =" << NmultKr<< "\n К-во аддитивных операций ="<<NaddKr<<"\n Общее к-во операций ="<<NmultKr+NaddKr;
Error(a,x);
getch();
}
//===========================================================================================================================
void GAUSS5(complex<double>a[5][5],complex<double>b[5])
{
complex<double>x[5],a1[1][1],b1[1],a2[2][2],b2[2],a3[3][3],b3[3],a4[4][4],b4[4];
if(z==2) for(i=0;i<2;i++) {for(j=0;j<2;j++) a2[i][j]=a[i][j];b2[i]=b[i];}
if(z==3) for(i=0;i<3;i++) {for(j=0;j<3;j++) a3[i][j]=a[i][j];b3[i]=b[i];}
if(z==4) for(i=0;i<4;i++) {for(j=0;j<4;j++) a4[i][j]=a[i][j];b4[i]=b[i];}
if(z==5)
{
for(i=0;i<4;i++)
{
b4[i]=b[i]-b[4]*a[i][4]/a[4][4];
for (j=0;j<4;j++) a4[i][j]=a[i][j]-a[4][j]*a[i][4]/a[4][4];
NaddGauss=NaddGauss+5;
NmultGauss=NmultGauss+10;
}
}
if(z>=4)
{
for(i=0;i<3;i++)
{
b3[i]=b4[i]-b4[3]*a4[i][3]/a4[3][3];
for(j=0;j<3;j++) a3[i][j]=a4[i][j]-a4[3][j]*a4[i][3]/a4[3][3];
NaddGauss=NaddGauss+4;
NmultGauss=NmultGauss+8;
}
}
if(z>=3)
{
for(i=0;i<2;i++)
{
b2[i]=b3[i]-b3[2]*a3[i][2]/a3[2][2];
for(j=0;j<2;j++) a2[i][j]=a3[i][j]-a3[2][j]*a3[i][2]/a3[2][2];
NaddGauss=NaddGauss+3;
NmultGauss=NmultGauss+6;
}
}
for(i=0;i<1;i++)
{
b1[i]=b2[i]-b2[1]*a2[i][1]/a2[1][1];
for(j=0;j<1;j++) a1[i][j]=a2[i][j]-a2[1][j]*a2[i][1]/a2[1][1];
x[0]=b1[0]/a1[0][0];
cout<<"\ x[0]=";
if(imag(x[0])!=0) cout<<x[0];
else cout<<real(x[0]);
NaddGauss=NaddGauss+2;
NmultGauss=NmultGauss+5;
}
if(z>=2)
{
x[1]=(b2[1]-a2[1][0]*x[0])/a2[1][1];
cout<< " х[1]=";
if(imag(x[1])!=0) cout<<x[1];
else cout<<real (x[1]);
NaddGauss++;
NmultGauss=NmultGauss+2;
}
if(z>=3)
{
x[2]=(b3[2]-a3[2][0]*x[0]-a3[2][1]*x[1])/a3[2][2];
cout<<" х[2]= ";
if(imag(x[2])!=0) cout<<x[2];
else cout<<real(x[2]);
NaddGauss=NaddGauss+2;
NmultGauss=NmultGauss+3;
}
if(z>=4)
{
x[3]=(b4[3]-a4[3][0]*x[0]-a4[3][1]*x[1]-a4[3][2]*x[2])/a4[3][3];
cout<<" х[3]= ";
if(imag(x[3])!=0) cout<<x[3];
else cout<<real(x[3]);
NaddGauss=NaddGauss+3;
NmultGauss=NmultGauss+4;
}
if(z==5)
{
x[4]=(b[4]-a[4][0]*x[0]-a[4][1]*x[1]-a[4][2]*x[2]-a[4][3]*x[3])/a[4][4];
cout<<" х[4]= ";
if(imag(x[4])!=0) cout<<x[4];
else cout<<real(x[4]);
NaddGauss=NaddGauss+4;
NmultGauss=NmultGauss+5;
}
}
void GAUSS()
{
complex<double> a[5][5],b[5],x[5];
for(i=0;i<5;i++) for (j=0;j<5;j++) a[i][j]=0;
for(i=0;i<z;i++) b[i]=0;
NaddGauss=0;
NmultGauss=0;
cout<<"\n Введите размерность от 2 до 5: \n";
cin>>z;
cout<<"\n Введите коэффициенты при x: \n";
for(i=0;i<z;i++) for(j=0;j<z;j++) cin>>a[i][j];
cout<<"\n Введите матрицу значений: \n";
for(i=0;i<z;i++) cin>>b[i];
GAUSS5(a,b); //вызов функции рассчета корней
cout<<"\n К-во аддитивных операций=" <<NaddGauss<<"\n К-во мультипликативных операций ="<<NmultGauss<<"\n Общее к-во операций ="<<NaddGauss+NmultGauss;
Error(a,b);
getch();
}
//=================================================================================================================
//функция вычисления погрешности
complex<double> Error (complex<double> a[5][5], complex<double> b[5])
{
int m, n, q, w;
complex<double> D, osh;
complex<double> dop[5][5], obr[5][5], deltab[5];
for(i=0;i<5;i++) for(j=0;j<5;j++) {dop[i][j]=0; obr[i][j]=0;}
for(i=0;i<5;i++) deltab[i]=0.02;
if(z==2) {complex<double> b[2][2]; for(i=0;i<z;i++) b[i][j]=a[i][j];D=dva(b);}
if(z==3) {complex<double> b[3][3]; for(i=0;i<z;i++) b[i][j]=a[i][j];D=tri(b);}
if(z==4) {complex<double> b[4][4]; for(i=0;i<z;i++) b[i][j]=a[i][j];D=chetire(b);}
if(z==5) D=pyat(a);
for(i=0;i<z;i++)
{
for(j=0;j<z;j++)
{
m=0;
n=0;
if(z==2) for (q=0;q<z;q++) for(w=0;w<z;w++) if(q!=i && w!=j) {dop[i][j]=pow((double)-1,i+j)*a[q][w]/D;obr[j][i]=dop[i][j];}
if(z==3) {complex<double> minor[2][2];for(q=0;q<z;q++) {n=0; for(w=0;w<z;w++) if(q!=i && w!=j) {minor[m][n]=a[q][w];n++;} if(q!=i && w!=j) m++;}
dop[i][j]=pow((double)-1,i+j)*dva(minor)/D; obr[j][i]=dop[i][j];}
if(z==4) {complex<double> minor[3][3];for(q=0;q<z;q++) {n=0; for(w=0;w<z;w++) if(q!=i && w!=j) {minor[m][n]=a[q][w];n++;} if(q!=i && w!=j) m++;}
dop[i][j]=pow((double)-1,i+j)*tri(minor)/D;obr[j][i]=dop[i][j];}
if(z==5) {complex<double> minor[4][4];for(q=0;q<z;q++) {n=0;for(w=0;w<z;w++) if(q!=i && w!=j) {minor[m][n]=a[q][w];n++;} if(q!=i && w!=j)m++;}
dop[i][j]=pow((double)-1,i+j)*chetire(minor)/D;obr[j][i]=dop[i][j];}
}
}
cout<<"\nОбратная :";
cout<<"\n ";
for(i=0;i<z;i++)
{
cout<<"\n ";
for(j=0;j<z;j++)
{
cout<<" ";
if(imag(obr[i][j])!=0) cout<<obr[i][j];
else cout<<real(obr[i][j]); }
}
for(i=0;i<z;i++)
osh=normMatr(a)*normMatr(obr)*normVect(deltab)/normVect(b);
cout<<"\n Погрешность = "<<real(osh);
return osh;
}
//=============================================================================================
complex<double> normMatr(complex<double> a[5][5])
{
double norm=0;
for(i=0;i<z;i++)
for(j=0;j<z;j++) norm+=abs(a[i][j])*abs(a[i][j]);
return sqrt(norm);
}
complex<double> normVect(complex<double> b[5])
{
double norm=0;
for(i=0;i<z;i++) norm+=abs(b[i])*abs(b[i]);
return sqrt(norm);
}
//===============================================================================
void main()
{
setlocale(LC_ALL, "russian");
char t;
cout<<" Выбор пункта: \n Методы: ";
cout<< "\n 1.Крамер\n 2.Гаусс\n";
cin>>t;
switch(t)
{
case '1': KRAMER();break;
case '2': GAUSS();break;
}
getch();
}