Метод Гаусса решения СЛАУ
Доброго времени суток . Пишу программу , которая должна решать системы линейных уравнений методом Гаусса(методом ведущего элемента в столбце) . Код написан , но почему-то возникают проблемы "на ровном месте" , например , неправильно находит ведущий элемент , хотя алгоритм вроде бы правильный . Буду благодарен за любую помощь .
Код:
#include<iostream.h>
#include<math.h>
void swap_rows(double **, int, int, int);
void swap(double *, int, int);
void div_row(double **, int, int, double);
void Gauss_method(double **, double *, double *, int);
int main(){
int i, j, n;
cout<<"System : Ax=b (A - n x n, b - n x 1)\n";
cout<<"Enter n : ";
cin>>n;
double ** A = new double * [n];
for(i = 0; i < n; ++i) {
A[i] = new double [n];
}
double * b = new double[n];
cout<<"Enter the matrix A : \n";
for(i = 0; i < n; ++i) {
for(j = 0; j < n; ++j) {
cin>>A[i][j];
}
}
cout<<"Enter the vector b : \n";
for(i = 0; i < n; ++i) {
cin>>b[i];
}
// Решение системы
double * x = new double [n];
Gauss_method(A, b, x, n);
cout<<"The solution : \n";
for(i = 0; i < n; ++i) {
cout<<"x_"<<i+1<<" = "<<x[i]<<"\n";
}
cout<<endl;
// Очистка памяти
for(i = 0; i < n; ++i) {
delete A[i];
}
delete [] A;
delete [] b;
delete [] x;
return 0;
}
void swap_rows(double ** A, int n, int i, int j) {
double tmp = 0;
for(int k = 0; k < n; ++k) {
tmp = A[i][k];
A[i][k] = A[j][k];
A[j][k] = tmp;
}
}
void swap(double * b, int i, int j) {
double tmp = b[i];
b[i] = b[j];
b[j] = tmp;
}
void div_row(double ** A, int n, int i, double d) {
for(int k = 0; k < n; ++k){
A[i][k] /= d;
}
}
void Gauss_method(double ** A, double * b, double * x, int n) {
int i, j, k;
int row;
double d = 0, tmp = 0;
for(k = 0; k < n; ++k) {
for(i = k; i < n; ++i) {
if(fabs(A[i][k]) > d){
d = fabs(A[i][k]);
row = i;
}
}
if(row != k) {
swap(b, row, k);
swap_rows(A, n, row, k);
}
b[k] /= A[k][k];
div_row(A, n, k, A[k][k]);
for(i = k + 1; i < n; ++i){
tmp = A[i][k];
b[i] -= tmp * b[k];
for(j = k; j < n; ++j) {
A[i][j] -= tmp * A[k][j];
}
}
}
for(i = n - 1; i >= 0; --i) {
x[i] = b[i];
for(j = i+1; j < n; ++j){
x[i] -= x[j] * A[i][j];
}
}
}
#include<math.h>
void swap_rows(double **, int, int, int);
void swap(double *, int, int);
void div_row(double **, int, int, double);
void Gauss_method(double **, double *, double *, int);
int main(){
int i, j, n;
cout<<"System : Ax=b (A - n x n, b - n x 1)\n";
cout<<"Enter n : ";
cin>>n;
double ** A = new double * [n];
for(i = 0; i < n; ++i) {
A[i] = new double [n];
}
double * b = new double[n];
cout<<"Enter the matrix A : \n";
for(i = 0; i < n; ++i) {
for(j = 0; j < n; ++j) {
cin>>A[i][j];
}
}
cout<<"Enter the vector b : \n";
for(i = 0; i < n; ++i) {
cin>>b[i];
}
// Решение системы
double * x = new double [n];
Gauss_method(A, b, x, n);
cout<<"The solution : \n";
for(i = 0; i < n; ++i) {
cout<<"x_"<<i+1<<" = "<<x[i]<<"\n";
}
cout<<endl;
// Очистка памяти
for(i = 0; i < n; ++i) {
delete A[i];
}
delete [] A;
delete [] b;
delete [] x;
return 0;
}
void swap_rows(double ** A, int n, int i, int j) {
double tmp = 0;
for(int k = 0; k < n; ++k) {
tmp = A[i][k];
A[i][k] = A[j][k];
A[j][k] = tmp;
}
}
void swap(double * b, int i, int j) {
double tmp = b[i];
b[i] = b[j];
b[j] = tmp;
}
void div_row(double ** A, int n, int i, double d) {
for(int k = 0; k < n; ++k){
A[i][k] /= d;
}
}
void Gauss_method(double ** A, double * b, double * x, int n) {
int i, j, k;
int row;
double d = 0, tmp = 0;
for(k = 0; k < n; ++k) {
for(i = k; i < n; ++i) {
if(fabs(A[i][k]) > d){
d = fabs(A[i][k]);
row = i;
}
}
if(row != k) {
swap(b, row, k);
swap_rows(A, n, row, k);
}
b[k] /= A[k][k];
div_row(A, n, k, A[k][k]);
for(i = k + 1; i < n; ++i){
tmp = A[i][k];
b[i] -= tmp * b[k];
for(j = k; j < n; ++j) {
A[i][j] -= tmp * A[k][j];
}
}
}
for(i = n - 1; i >= 0; --i) {
x[i] = b[i];
for(j = i+1; j < n; ++j){
x[i] -= x[j] * A[i][j];
}
}
}
Вот вариант, написанный на C#. Если кому-нибудь вдруг понадобится.
Самая простая форма реализации данного метода - без выбора ведущего элемента.
Не совсем уверен, что правильно написал.
// Метод Гаусса
private double[] Gauss(double[,] B, double[] RightPart)
{
double[] x = new double[N];
double R;
try
{
// Прямой ход
for (int q = 0; q < N; q++)
{
R = 1 / B[q, q];
B[q, q] = 1;
for (int j = q + 1; j < N; j++)
B[q, j] *= R;
RightPart[q] *= R;
for (int k = q + 1; k < N; k++)
{
R = B[k, q];
B[k, q] = 0;
for (int j = q + 1; j < N; j++)
B[k, j] = B[k, j] - B[q, j] * R;
RightPart[k] = RightPart[k] - RightPart[q] * R;
}
}
}
catch (DivideByZeroException)
{
IER = q;
return x;
}
// Обратный ход
for (int q = N - 1; q >= 0; q--)
{
R = RightPart[q];
for (int j = q + 1; j < N; j++)
R -= B[q, j] * x[j];
x[q] = R;
}
return x;
}
Самая простая форма реализации данного метода - без выбора ведущего элемента.
Не совсем уверен, что правильно написал.
// Метод Гаусса
private double[] Gauss(double[,] B, double[] RightPart)
{
double[] x = new double[N];
double R;
try
{
// Прямой ход
for (int q = 0; q < N; q++)
{
R = 1 / B[q, q];
B[q, q] = 1;
for (int j = q + 1; j < N; j++)
B[q, j] *= R;
RightPart[q] *= R;
for (int k = q + 1; k < N; k++)
{
R = B[k, q];
B[k, q] = 0;
for (int j = q + 1; j < N; j++)
B[k, j] = B[k, j] - B[q, j] * R;
RightPart[k] = RightPart[k] - RightPart[q] * R;
}
}
}
catch (DivideByZeroException)
{
IER = q;
return x;
}
// Обратный ход
for (int q = N - 1; q >= 0; q--)
{
R = RightPart[q];
for (int j = q + 1; j < N; j++)
R -= B[q, j] * x[j];
x[q] = R;
}
return x;
}