Поможем написать учебную работу
Если у вас возникли сложности с курсовой, контрольной, дипломной, рефератом, отчетом по практике, научно-исследовательской и любой другой работой - мы готовы помочь.
Если у вас возникли сложности с курсовой, контрольной, дипломной, рефератом, отчетом по практике, научно-исследовательской и любой другой работой - мы готовы помочь.
МИНИСТЕРСТВО ОБРАЗОВАНИЯ РЕСПУБЛИКИ БЕЛАРУСЬ
Полоцкий государственный университет
Кафедра технологий программирования
Отчет
По лабораторной работе №2
по курсу «Численные методы в инженерных расчетах»
«Метод Холецкого»
Выполнила:
студент группы 10-ИТ-1
Бойченко А.В.
Проверила:
Мурашкевич О.Н.
Полоцк, 2012 г.
Цель: изучить метод Холецкого (метод квадратного корня).
Формулы метода квадратного корня несколько сложнее, чем в методе Гаусса. К тому же этот алгоритм может оказаться численно неустойчивым. Однако в двух случаях можно гарантировать устойчивость алгоритма для положительно определенных матриц А и матриц с диагональным преобладанием.
Теорема. Если матрица А симметричная (А=АТ), вещественная и положительно определенная А>0, тогда
А=SТS (3.1)
Для положительной определенности матрицы необходимо и достаточно, чтобы все ее главные миноры были положительны. Это условие на практике проверить сложно. Обычно достаточно требования положительности диагональных элементов и их диагонального преобладания, т.е.
aii > , i=1,2,…,n.
Рассмотрим сначала случай симметричной положительно определенной матрицы А. Подставим (3.1) в уравнение AX=b и введем обозначение SX=Y.
В этом случае решение системы линейных алгебраических уравнений (1) сводится к последовательному решению двух систем с треугольными матрицами
SТY=b (3.2)
SX=Y (3.3)
Т.к. матрица SТ={s*i,j =0, j>i, i=1,2,..,n; j=2,…,n} системы (3.2) является нижней треугольной, т.е. система выглядит следующим образом
, и можно сразу выписать ее решение
i=2,3,…,n (3.4)
Определив, таким образом, вектор Y, можем определить из системы (3.3) искомое решение. Это решение находим обратным ходом метода Гаусса, т.к. матрица S верхняя треугольная.
Имеем:
, i=n-1,n-2,…,1 (3.5)
Вычисления по формулам (3.4)-(3.5) требуют знания коэффициентов матрицы S. Их значения найдем из условия (3.1), т.е. А=STS.
Рассмотрим в качестве примера случай, когда матрица является квадратной матрицей второго порядка.
a11>0 a22>0 и на основании (3.1) имеем:
Для выполнения условия (3.1) необходимо
a11=(s11)2 a12=s11s12 a22= (s12)2 + (s22)2
Из этой системы можем найти конкретное значение коэффициентов матрицы S.
s11= , s12 = s22 = ,
Рассуждая аналогичным образом для матрицы А произвольного порядка получим следующие соотношения.
Таким образом, решение системы (1) методом квадратного корня предусматривает
Листинг программы:
#include <vcl.h>
#pragma hdrstop
#include <iostream.h>
#include <fstream.h>
#include <conio.h>
#include <math.h>
int main()
{
int N,j,k,i;
float temp,del,l,q;
float a[10][10],f[10],x[10],y[10],s[10][10],st[10][10];
fstream f1,f2;
f1.open("input.txt",ios::in);
f2.open("output.txt",ios::out);
f1>>N;
for(i=0;i<N;i++) //ввод матрицы
for(j=0;j<N;j++)
f1>>a[i][j];
for (i=0;i<N;i++) //ввод вектор-столбца
f1>>f[i];
for(i=0;i<N;i++)
for(j=0;j<N;j++)
{
s[i][j]=0;
st[i][j]=0;
}
//подсчёт транформирующеей матрицы S,
//удовлетворяющей равенству A=s*DS
s[0][0]=sqrt(a[0][0]);
for(j=1;j<N;j++)
s[0][j]=a[0][j]/s[0][0];
for(i=1;i<N;i++)
{
q=0;
for(k=0;k<i;k++)
q+=(s[k][i]*s[k][i]);
s[i][i]=sqrt(a[i][i]-q);
for(j=i+1;j<N;j++)
{
l=0;
for(k=0;k<i;k++)
l+=(s[k][i]*s[k][j]);
s[i][j]=(a[i][j]-l)/s[i][i];
}
}
for(i=0;i<N;i++)
{ cout<<endl;
for(j=0;j<N;j++)
cout<<s[i][j]<<' ';}
getch();
cout<<endl;
for(i=0;i<N;i++)
for(j=0;j<N;j++)
st[i][j]=s[j][i];
//решение уранения S*Y=F
y[0]=f[0]/st[0][0];
//обратный ход, нахождение вектора-Y
for(i=1;i<N;i++)
{
del=0;
for(j=0;j<i;j++)
{
del+=(st[i][j]*y[j]);
}
y[i]=(f[i]-del)/st[i][i];
}
for(k=0;k<N;k++)
cout<<'y'<<k<<'='<<x[k]<<endl;
getch();
//решение уранения SX=Y
x[N-1]=y[N-1]/s[N-1][N-1];
//обратный ход, нахождение вектора-решения
for(i=N-2;i>=0;i--)
{
del=0;
for(j=i+1;j<N;j++)
{
del+=(s[i][j]*x[j]);
}
x[i]=(y[i]-del)/s[i][i];
}
//вывод результата
for(k=0;k<N;k++)
f2<<'x'<<k+1<<'='<<x[k]<<endl;
f1.close();
f2.close();
return 0;
}
Результат работы:
Проверка: