3.2. Решение задачи заданной вариантом используя написанную программу:
Решая задачу своего варианта с помощью написанной программы получено следующее приближённое решение:
Также наша программа находит вектор невязки для полученного решения:
Найдём 3 нормы вектора невязки:
Нормированный определитель вектора нормаль равен 1.02477e-008
Проверим систему на устойчивость. Для этого прибавим к случайному элементу матрицы А случайную погрешность в 4-м знаке, получим решение:
При внесении погрешности в 4-й знак элемента матрицы получаем погрешность решения в единицах.
3.3. Решение задачи заданного варианта в математическом пакете maple
Решаем задачу, используя математический пакет maple, получаем вектор решения:
Вычисляем вектор невязки для полученного решения:
>
4. Вывод
Исходя из значения вектора невязки можно сказать, что погрешность вычислений, используя написанную программу использующую метод Гаусса-Жордана составляет порядка 11-го знака после запятой. Математический пакет maple решает приближённо эту же систему с погрешностью порядка более 6-го знака после запятой. Рассматривая решение типовой задачи можно сказать что созданная программа будет решать СЛАУ у которых координаты вектора решений будут целыми, с крайне малой погрешностью, вплоть до того что будет получаться точное решение.
Также, про матрицу А можно сказать следующее: она не является устойчивой, в следствии того что при внесении случайной погрешности порядка 4-го знака после запятой в случайный элемент матрицы получили погрешность в 11-м знаке после запятой. Исходя из того что нормированный определитель лежит в промежутке от 0 до 1 мы можем сказать что матрица плохо обусловлена.
Приложение а
Листинг программы
#include <iostream>
# include <stdio.h>
# include <fstream>
# include <math.h>
# include <time.h>
# include <locale>
using namespace std;
void main()
{
setlocale(LC_ALL, "Russian");// подключение русской раскладки для вывода всех сообщений на русском языке
int i,j,n,k,max;
long double S,p,n1=0,n2=0,n3=0;
//long double MatrixA[4][5]={67.09, -40.45, -16.77, 9.87, 1.05,-40.45, 24.66, 9.86, -5.92, 21.33,-16.77, 9.86, 4.94, -2.97, 7.11, 9.87, -5.92, -2.97, 1.98, -0.23};
long double MatrixA[4][5]={1, -1, 1, 1, 2,2, -1, 2, 1, 4,3, 1, 2, 2, 8,4, -2, 1, 3, 6};
// массив хранит в себе матрицу А в столбцах от 0-го до 3-го и матрицу F в 4-м
long double A[4][5],r[4];
n=4;
srand(time(0));
long double rand1=0;
for (i=0; i<n; i++)// создём копию матрица А и F
{
for (j=0; j<n+1; j++)
{
rand1=(rand()%10+1)/10000.0;
// A[i][j]=MatrixA[i][j]+rand1;// добавляем к случайному элементу матрецы случайную погрешность в 4-м знаке
// MatrixA[i][j]=A[i][j];
A[i][j]=MatrixA[i][j];
}
}
long double det=1;// принимаем начальное значени определителя за единицу, так как будем вычислять определитель используя методгауса жордана
ofstream f("reshenie.txt");
f<<"матрица А:\n";
for (i=0; i<n; i++)
{
for (j=0; j<n+1; j++)
{
cout<< MatrixA[i][j]<<" ";
if(j==4)f<<"\t";
f<<MatrixA[i][j]<<" ";
}
cout<< endl;
f<<endl;
}
// метод Гауса-Жордана
for(k=0;k<4;k++)
{
int m=0;
// находим наибольший элемент столбца на к-м шаге
max=MatrixA[k][k];
m=k;
for(i=k+1;i<4;i++)
{
if(max<MatrixA[i][k])
{
max=MatrixA[i][k];
n=i;
}
}
// меняем строки в матрице
if(m!=k)
{
for(i=k;i<5;i++)
{
p=MatrixA[k][i];
MatrixA[k][i]=MatrixA[m][i];
MatrixA[m][i]=p;
}
det=det*-1;// умножаем определитель на -1 если меняем строки(по свойству определителя)
}
// делим строку на элемент стоящий на главной диагонали при условии что этот элемент не равен 1
if(MatrixA[k][k]!=1)
{
det=det*MatrixA[k][k];//делим определитель на диагональный элемент по свойству определителей
for(i=4;i>k-1;i--)
{
MatrixA[k][i]=MatrixA[k][i]/MatrixA[k][k];
}
}
// вычитаем строку с которой работаем из остальных строк( вычитаем с конца, для того чтобы не обнулять элементы под и над элементом стоящим на главной диагонали)
for (i=0;i<5;i++)
{
if (i!=k)
{
for(j=4;j>k-1;j--)
{
MatrixA[i][j]=MatrixA[i][j]-MatrixA[k][j]*MatrixA[i][k];
}
}
}
}
// выводим полученную после всех преобразований матрицу на экран
cout<< endl;
f<<endl<<"преобразованная матрица А:"<<endl;
for (i=0; i<4; i++)
{
for (j=0; j<5; j++)
{
cout<< MatrixA[i][j]<<" ";
if(j==4)f<<"\t";
f<<MatrixA[i][j]<<" ";
}
cout<< endl;
f<<endl;
}
// считаем вектор невязки
long double var=0;
f<<"вектор невязки:"<<endl;
for (i=0; i<4; i++)
{
var=0;
for (j=0; j<4; j++)
{
var=var+A[i][j]*MatrixA[j][4];
}
r[i]=A[i][4]-var;
cout << r[i]<<endl;
f<<r[i]<<endl;
}
// вычисляем 1-ю норму
f<<"норма 1:\n";// записываем в файл номер нормы, и переводим курсор на следующую строку
var=0;// var хранит значение двойной суммы стоящей под корнем
for (int i=0; i<n; i++)
{
var=var+r[i]*r[i];
}
n1=sqrt(var);
f<< n1<<endl;// записываем в файл 1-ю норму вектора невязки, и переводим курсор на следующую строку
// вычисляем 2-ю норму
f<<"норма 2:\n";// записываем в файл номер нормывектора невязки, и переводим курсор на следующую строку
var=0;// var хранит значение суммы
for (int i=0; i<n; i++)
{
var=var+abs(r[i]);
}
f<<n2<<endl;// записываем в файл 2-ю норму вектора невязки, и переводим курсор на следующую строку
// вычисляем 3-ю норму вектора невязки
f<<"норма 3:\n";// записываем в файл номер нормы вектора невязки, и переводим курсор на следующую строку
for (int i=0; i<n; i++)
{
var=abs(r[i]);
if(var>n3)n3=var;
}
f<<n3<<endl;// записываем в файл 3-ю нормувектора невязки, и переводим курсор на следующую строку
cout<< n1<<endl<<n2<<endl<<n3<<endl;
// вычисляем норму матрицы в пространстве Rn
for (int i=0; i<4; i++)
{
for(j=0;j<4;j++)
{
var=var+A[i][j]*A[i][j];
}
}
long double n4=0;
n4=sqrt(var);
f<<"нормированный определитель:\n"<<det/n4;
cout<<"нормированный определитель:\n"<<det/n4;
f.close();
//cout<<endl<<detA<<endl;
system("pause");
}