Алгоритм решения СЛАУ методом Гаусса C++
Автор: SysaninФев 16
Попросили как то найти алгоритм решения СЛАУ (системы линейных алгебраических уравнений) методом Гаусса на pascal'e или хотя бы на каком другом языке. Но все алгоритмы были оформлены: а ля решения по лабораторным работам, естественно за звонкую монету. В итоге недавно мне самому понадобилось написать лабу на C++ по этому алгоритму. Вот и решил выложить.
Сам алгоритм математического решения приводить не буду (там всякие прямые прогонки, обратные), это легко гуглится при желании. Тока скажу саму суть алгоритма: составляем квадратную матрицу n*n из коэфицентов уравнений. Далее приводим матрицу к треугольному типу. Я делал это путем деления коэфицентов 1,2...(n-1) -го из первого уравнения на соответствующие коэфиценты 2,3..n-го уравнений. Полученное частное я умножал на первое уравнение и отнимал от текущего (2,3..n-го). Если на диагонали лежал нулевой элемент, то я менял строку местами со следующий у которой был не нулевой элемент на диагонали.
После такой прогонки мы получаем треугольную матрицу. Единственное, что не всякую матрицу можно привести к треугольному виду. В частном случае возможны нулевые элементы на главной диагонали. Тогда решать методом Гаусса нельзя.
Вот собственно икод, сдобренный комментариями. Комментарии /**/ содержат руководство к тому что вам надо дописать.
// меняем местами строки в массиве с коэфицентами уравнения
void exchange(double **mas, int razmer,int from, int to){
for(int i=0; i < razmer; i++){
mas[from][i] = mas[from][i] + mas[to][i];
mas[to][i] = mas[from][i] - mas[to][i];
mas[from][i] = mas[from][i] - mas[to][i];
}
}
// определяем ближайшую строку с не 0 диоганалью вниз
int get_n_no_empty(double **a, int razmer, int begin){
for(int i=begin+1; i < razmer; i++){
if(a[i][i] != 0){
return i;
}
}
cout << "Данная система не решается методом Гауса" << endl;
getch();
exit(1);
}
// функция ждет пока юзер прочтет и нажмет на кнопку, не обязательная функция
void pause(){
cout << "Нажмите любую кнопку для продолжения" << endl;
getch();
}
Непосредственный код для выполнения:
double **a, // матрица с коэфицентами
*x, // х
*b; // свободные члены
/* Здесь переменной razmer необходимо присвоить размер матрицы */
// создание динамических переменных
b = new double[razmer];
x = new double[razmer];
a = new double*[razmer];
for(int i=0; i < razmer; i++){
a[i] = new double[razmer];
}
/* Здесь нам надо присвоить значения матрице а, и заполнить свободные члены b */
// начинаем решать)
// частный случай, если размер матрицы == 1
if(razmer == 1){
if(a[0][0] != 0){
cout << "x=" << b[0]/a[0][0] << endl;
}else{
cout << "Х любое" << endl;
}
getch();
exit(0);
}
// прямой проход
// идем слево на право по коэфицентам х
for(int k=0; k < razmer-1; k++){
// вычитаем по строкам
for(int m=k+1; m < razmer; m++){
// если на диагонали элемент = 0, то поменяем местами строки
if(a[m][m] == 0) exchange(a, razmer, m, get_n_no_empty(a, razmer, m));
double koeficent=a[m][k]/a[0][k];
//вычисление новых коэфицентов уравнения
b[m] = b[m] - b[0] * koeficent;
for(int z=0; z < razmer; z++){
a[m][z] = a[m][k] - a[0][k] * koeficent;
}
}
}
// ищем решения
for(int m=razmer-1; m >= 0; m--){
double sum=0;
// идем по строке спаво налево, считая сумму корень*коэфицент, до текущего корня
for(int i=razmer-1; i > m; i--){
sum += x[i] * a[m][i];
}
x[m] = (b[m] - sum)/a[m][m];
}
// вывод решений
cout << "###Решения:" << endl;
for(int i=0; i < razmer; i++){
cout << "x[" << i << "]=" << x[i] << endl;
}
pause();
Для проверки полученных значений можно использовать этот код
for(int m=0; m < razmer; m++){
double sum=0;
for(int i=0; i < razmer; i++){
sum += a[m][i] * x[i];
}
cout << m << ":" << (sum - b[m]) << endl;
}
pause();
Для корректной работы подключите следующие модули:
#include <conio.h>
#include <tchar.h>
#include <iostream>
#include <stdlib.h>
#include <time.h>
7 коммент.
Пишет Loopy | дата: 25 февраля 2009 в 10:25
Лучше для такких штук юзать классы( эт на будущее ), а вообще чтобы не произошло того что диагональный элемент равен 0 и ошибка округления была меньше лучше юзать алгоритм с выбором гланого элемента хотябы по строке, тогда решения не будет существовать только в случае невырожденности матрицы. Ну или строку
if(a[m][m] == 0) exchange(a, razmer, m, get_n_no_empty(a, razmer, m));
заменить на
if(abs(a[m][m]) <E) exchange(a, razmer, m, get_n_no_empty(a, razmer, m));
где E выбирается в зависимости от порядка входных данных.;)
а вообще переходи на std классы.
Пишет Sysanin | дата: 25 февраля 2009 в 12:16
ну классы да:) тока я еще плохо вкурил модель проектирования на них
а на счет алгоритма – ты прав, но как бы мне надо придерживаться методички:) хотя замену надо будет сделать.
Пишет Sysanin | дата: 25 февраля 2009 в 14:46
тока, наверно, вместо abs(a[m][m]), надо fabs, abs же тока для целых чисел…
Пишет Loopy | дата: 26 февраля 2009 в 21:38
посмотри math.h там она перегружена для разных типов;)
Пишет Sysanin | дата: 26 февраля 2009 в 22:05
хм… значит в старой версии билдера она была не перегруженной, мб вобщем
Пишет bazz | дата: 26 апреля 2010 в 15:18
// прямой проход
// идем слево на право по коэфицентам х
for(int k=0; k < razmer-1; k++){
// вычитаем по строкам
for(int m=k+1; m < razmer; m++){
// если на диагонали элемент = 0, то поменяем местами строки
if(a[m][m] == 0) exchange(a, razmer, m, get_n_no_empty(a, razmer, m));
double koeficent=a[m][k]/a[k][k];
//вычисление новых коэфицентов уравнения
b[m] = b[m] – b[k] * koeficent;
for(int z=0; z < razmer; z++){
a[m][z] = a[m][k] – a[k][z] * koeficent;
}
}
}
Пишет Андрей | дата: 13 октября 2010 в 20:43
проверьте плизз правильность написания прямого хода метода
#include
#include
#include
#define eps 0.0000000001
class CMatrix
{private: int m;
int n;
double*p;
public: CMatrix (){};
CMatrix(int,int);
friend CMatrix operator+(CMatrix&,CMatrix&);
friend CMatrix operator*(CMatrix&,CMatrix&);
friend CMatrix gauss(CMatrix&,CMatrix&);
void input();
void output(int,int);
void output1(int,int);
void input1();
double det();
};
void gauss();
void adddouble ();
void multdouble ();
void det();
void main ()
{int k;
clrscr();
cout<<» add 1\n multiply 2\n schutat’ opredelitel 3\n metod gaussa 4\n»<>k;
if(k==1) adddouble();
else if (k==2) multdouble();
else if (k==3) det();
else if (k==4) gauss();
}
CMatrix::CMatrix(int a,int b)
{m=a;
n=b;
p=new double [m*n];
}
void CMatrix::output(int x, int y)
{int i;
int j;
for(i=0;i<m;i++)
{
for(j=0;j<n;j++)
{
gotoxy(5*j+x,2*i+y);
cout<<p[i*n+j];
}
}
}
void CMatrix::output1(int x, int y)
{int i;
int j;
for(i=0;i<m;i++)
{
for(j=1;j<n;j++)
{
gotoxy(5*j+x,2*i+y);
cout<<p[i*n+j];
}
}
}
void det()
{int n;
double d;
cout<<»vvedite razmer matrici n» <>n;
CMatrix A(n,n);
clrscr();
A.input();
A.output(31,9);
d=A.det();
gotoxy(28,15);
cout<<»opredel raven «<<d;
getch();
}
void adddouble ()
{int m,n;
clrscr();
cout<<»vvedite razmeri matriz m i n»<>m>>n;
CMatrix A(m,n);
CMatrix B(m,n);
CMatrix C(m,n);
A.input();
B.input();
clrscr();
C=A+B;
getch();
}
void multdouble ()
{int m,n,l;
cout<<»vvedite m, n, l»<>m>>n>>l;
CMatrix A(m,n);
CMatrix B(n,l);
CMatrix C(m,l);
A.input();
B.input();
clrscr();
C=A*B;
A.output(1,1);
B.output(41,1);
C.output(21,9);
getch();
}
void CMatrix::input()
{int i, j;
cout<<»vvedite info»<<endl;
for(i=0;i<m;i++)
for(j=0;j>p[i*n+j];
}
void CMatrix::input1()
{int i, j;
cout<<»vvedite info»<<endl;
for(i=1;i<m;i++)
for(j=1;j>p[i*n+j];
}
CMatrix operator+(CMatrix&obj1,CMatrix&obj2)
{int i,j;
CMatrix temp(obj1.m,obj1.n);
for(i=0;i<obj1.m;i++)
{
for(j=0;j<obj1.n;j++)
{
temp.p[i*obj1.n+j]=obj1.p[i*obj1.n+j]+obj2.p[i*obj1.n+j];
gotoxy(5*j+1,2*i+1);
cout<<obj1.p[i*obj1.n+j];
gotoxy(5*j+41,2*i+1);
cout<<obj2.p[i*obj1.n+j];
gotoxy(5*j+21,2*i+8);
cout<<temp.p[i*obj1.n+j];
}
}
return temp;
}
CMatrix operator*(CMatrix&obj1,CMatrix&obj2)
{
int i,j,k;
double s;
CMatrix temp(obj1.m,obj2.n);
for(i=0;i<obj1.m;i++)
{
for(j=0;j<obj2.n;j++)
{
s=0;
for(k=0;k2)
{
CMatrix OBR(n-1,n-1);
for(i=0;i<n;i++)
{ for(j=1;j<n;j++)
{for (k=0;k<n;k++)
if(ki)
OBR.p[(j-1)*(n-1)+(k-1)]=p[j*n+k];
}
znak*=-1;
d+=(p[i]*OBR.det()*znak);
}
}
return d;
}
void gauss()
{
int n;
cout<<»vvedite n «<>n;
CMatrix a (n,n);
CMatrix b (n,1);
CMatrix x (1,n);
a.input1();
b.input1();
x=gauss(a,b);
gotoxy (28,15);
cout<<»otvet = «;
x.output1(38,15);
getch;
}
CMatrix gauss(CMatrix&a,CMatrix&b)
{
CMatrix x(1,a.n);
int i;
int j;
int k;
int n=a.n;
for (k=1;k<n;k++)//nomer obnul stolbcha
{
/*if(fabs(a.p+k*(a.n+1)+k)<eps )
{
}*/
for(i=k;i<n;i++)//i nomer obr str
{ clrscr();
a.output1 (1,1);
b.output1 (37,1);
for(j=k;j<n;j++)// nomer stolbcha v obr stroke
{
a.p[i*n+j]=a.p[i*n+j]-a.p[i*n+k]/a.p [k*n+k] * a.p[k*n+j];
}
b.p[i]=b.p[i]-a.p[i*n+k]/a.p[k*n+k]*b.p[k];
a.output1(40,1);
b.output1(70,1);
getch;
}
}
return x;
}