Лабораторная работа №8

Постановка задачи

Решить систему линейных алгебраических уравнений методом сопряженных градиентов.

Работу выполнить на задачах лабораторной работы № 7

Теоретические сведения

Метод сопряженных градиентов предназначен для решения систем линейных алгебраических уравнений с симметричной положительно определенной матрицей. Матрица называется положительно определенной, если скалярное произведениедля всех ненулевых векторов. Для того чтобы матрицабыла положительно определенной необходимо и достаточно, чтобы все ее главные миноры были положительны.

Рассмотрим систему линейных алгебраических уравнений (4.1) с положительно определенной симметричной матрицей. Пусть - решение этой системы. Покажем, что в этом случаедает минимум функционала. Имеем

.

Так как матрица положительно определенная, то из последнего равенства следует, что решениесистемы (4.1) дает минимум функционала.

Будем искать этот минимум итерационным методом. Для этого предположим, что заданы два вектора и. Последовательные приближенияк решению системы (4.1) будем строить по формулам

, (4.22)

, (4.23)

где . (4.24)

Параметры ибудем определять из условия минимума функционалав плоскости, проходящей через точкуи натянутой на векторыи, т.е. минимум функционала ищется на множестве. Преобразуем наш функционал на этом множестве

.

Для определения минимума используем необходимое условие экстремума функции. Продифференцируемпоии приравняем к нулю

,

.

Т.о., получаем систему уравнений

, (4.24)

. (4.25)

Так как

, то последнюю систему можно представить в виде

Это система двух уравнений с двумя неизвестными и. Решая ее, получаем

,

Покажем, что метод сопряженных градиентов всегда сходится к решению системы . Учитывая (4.23) - (4.25) и ортогональность векторови, получим

.

Так как векторы иортогональны, то из выражения для коэффициентаимеем

,

где - наибольшее собственное значение матрицы. Окончательно получаем

.

Таким образом, последовательность монотонно убывает и ограничена снизу значением. Так как всякая монотонно убывающая последовательность сходится, т.е., то по признаку Коши имеем. Тогда из последнего неравенства следует, что, а это означает, что последовательность векторовсходится к решениюсистемы (4.1).

Выберем начальные вектора идля метода сопряженных градиентов. Пусть- произвольный вектор. Построим вектор, где векторвыбирается по направлениютак, чтобы векторыибыли ортогональными, т.е.. Обозначим. Тогда

.

Для выполнения условия ортогональности векторов идостаточно потребовать, чтобы

.

Можно показать, что при таком выборе векторов иитерации метода сопряженных градиентов обрываются не позднее, чем на-ом шаге, давая точное решение системы (4.1) (- порядок системы). В сочетании с методом подавления компонент метод сопряженных градиентов позволяет находить решение за число итераций, значительно меньшее.

Основной операцией метода сопряженных градиентов является вычисление произведений и, причем ни в какие другие операции матрицане входит. Поэтому максимальный порядок решаемых на ЭВМ систем уравнений зависит только от объема информации, необходимой для умножения матрицына вектор. Значит, применение метода сопряженных градиентов будет тем эффективнее, чем эффективнее составлена программа умножения матрицы на вектор.

Код программы

Unit1.cpp:

//---------------------------------------------------------------------------

#include <vcl.h>

#pragma hdrstop

#include "Unit1.h"

#include "math.h"

//---------------------------------------------------------------------------

#pragma package(smart_init)

#pragma resource "*.dfm"

TForm1 *Form1;

//---------------------------------------------------------------------------

__fastcall TForm1::TForm1(TComponent* Owner)

: TForm(Owner)

{

}

//---------------------------------------------------------------------------

int n;

//float pr=0,nor=0;

float sk_pr(float* q1, float* q2){

float pr=0;

for(int i=0;i<=n;i++)

pr+=q1[i]*q2[i];

return pr;

}

float norm(float* q){

float nor=0;

nor=sqrt(sk_pr(q,q));

return nor;

}

void __fastcall TForm1::Button1Click(TObject *Sender)

{

int i,j,k=0,n;

float t,t1,t2,e,q,d,z;

float A[10][10];

float y[10],b[10],x[10],r[10],p[10],u[10];

n = StrToInt(Edit1->Text);

for (i=0; i<n; i++)

for (j=0; j<n; j++)

A[i][j] = StrToFloat(SG1->Cells[j][i]) ;

for (j=0; j<n; j++)

b[j] = StrToFloat(SG2->Cells[0][j]) ;

for (i=0; i<n; i++)

for (j=0; j<n; j++)

if (A[i][j]!=A[j][i]) {

ShowMessage("Матрица не симметричная !");

return;

}

x[0]=1;

for (i=1; i<n; i++) x[i]=0;

i=0;

while (i<n) {

t=0;

for (j=0; j<n; j++) {

y[i]=t+A[i][j]*x[j];

t=y[i];

}

r[i]=y[i]-b[i];

i++;

}

t=0;

for (i=0; i<n; i++)

t=t+r[i]*r[i];

i=0;

while (i<n) {

t1=0;

for (j=0; j<n; j++) {

y[i]=t1+A[i][j]*r[j];

t1=y[i];

}

i++;

}

e=0;

for (i=0; i<n; i++)

e=e+y[i]*r[i];

q = -(t/e);

for (i=0; i<n; i++) p[i]=q*r[i];

for (i=0; i<n; i++) x[i]=x[i]+p[i];

k=1;

while (k<n) {

i=0;

while (i<n) {

y[i]=0;

for (j=0; j<n; j++) y[i]=y[i]+A[i][j]*x[j];

r[i]=y[i]-b[i];

i++;

}

i=0;

while (i<n) {

t=0;

for (j=0; j<n; j++) {

y[i]=t+A[i][j]*r[j];

t=y[i];

}

i++;

}

i=0;

while (i<n) {

t=0;

for (j=0; j<n; j++) {

u[i]=t+A[i][j]*p[j];

t=u[i];

}

i++;

}

t=0;

for (i=0; i<n; i++) t=t+r[i]*r[i];

t1=0;

for (i=0; i<n; i++) t1=t1+r[i]*p[i];

t2=0;

for (i=0; i<n; i++) t2=t2+y[i]*p[i];

e=0;

for (i=0; i<n; i++) e=e+y[i]*r[i];

d=0;

for (i=0; i<n; i++) d=d+u[i]*p[i];

z=(t*t2-t1*e)/(d*e-t2*t2);

q=(t1*t2-t*d)/(d*e-t2*t2);

for (i=0; i<n; i++) {

p[i]=z*p[i]+q*r[i];

x[i]=x[i]+p[i];

}

k++;

}

for (int i = 0; i < 1; i++)

for (int j=0; j<=n;j++){

TSG1->Cells[i][j] = FloatToStr(x[j]);

}

float f1[50];

for(int i=0;i<n; i++){

f1[i]=0;

for(int j=0;j<n;j++)

f1[i]+=A[i][j]*x[j];

}

for (int i = 0; i < 1; i++)

for (int j=0; j<n;j++){

TSG2->Cells[i][j] = FloatToStr(f1[j]);

}

}

//---------------------------------------------------------------------------

void __fastcall TForm1::UpDown1Click(TObject *Sender, TUDBtnType Button)

{

SG1->RowCount = StrToInt(Edit1->Text);

SG2->RowCount = StrToInt(Edit1->Text);

SG1->ColCount = StrToInt(Edit1->Text);

TSG1->RowCount = StrToInt(Edit1->Text);

TSG2->RowCount = StrToInt(Edit1->Text);

}

//---------------------------------------------------------------------------

void __fastcall TForm1::Button3Click(TObject *Sender)

{

SG1->Cells[0][0] = "3,6";

SG1->Cells[1][0] = "0,8";

SG1->Cells[2][0] = "-0,7";

SG1->Cells[0][1] = "0,8";

SG1->Cells[1][1] = "3,6";

SG1->Cells[2][1] = "0,9";

SG1->Cells[0][2] = "-0,7";

SG1->Cells[1][2] = "0,9";

SG1->Cells[2][2] = "3,3";

SG2->Cells[0][0] = "3,8";

SG2->Cells[0][1] = "0,4";

SG2->Cells[0][2] = "-1,76";

}

//---------------------------------------------------------------------------

void __fastcall TForm1::Button2Click(TObject *Sender)

{

int n = StrToInt(Edit1->Text);

for(int i = 0; i < n; i++)

for(int j = 0; j <= n; j++)

SG1->Cells[j][i] = "";

for(int i = 0; i < 1; i++)

for(int j = 0; j < n; j++){

TSG2->Cells[i][j] = "";

TSG1->Cells[i][j] = "";

SG2->Cells[i][j] = "";

}

}

//---------------------------------------------------------------------------

Результат работы программы