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

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

Решить систему линейных алгебраических уравнений методом ортогонализации

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

Этот метод используется для решения систем с произвольной невырожденной матрицей, но по скорости уступает многим прямым методам примерно в 1,5-3 раза. Метод основан на построении вспомогательной системы векторов, связанных с матрицей исходной системы уравнений и ортогональные в некоторой метрике.

Рассмотрим систему (4.1) и запишем ее в виде

, (4.19)

где ,. Обозначим,,. Тогда система (4.19) перепишется в виде

Решение системы уравнений с невырожденной матрицей сводится к нахождению такого вектора, который имеет последнюю координату, равную единице, и ортогонален к линейно независимым векторам. Ортогональность векторак векторамвлечет за собой ортогональность ко всему подпространству, натянутому на них, и, следовательно, к любому его базису. И наоборот, ортогональность векторак некоторому базису подпространствавлечет ортогональность ко всем векторам. Поэтому для решения системы достаточно построить ненулевой вектор, ортогональный к какому-нибудь базису подпространства. Если- последняя координата вектора, то.

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

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

,. (4.20)

На первом шаге метода положим ,. Пусть для некоторого шагауже построен ортогональный базиси ортонормированный базисподпространства. Векторбудем искать как линейную комбинацию векторов

Условие ортогональности вектора к ортогональным векторамили, что то же самое, к векторамдает. Поэтому

. (4.21)

Таким образом, с помощью этого итерационного процесса и соотношения (4.20) строится ортонормированная система векторов . Векторыявляются базисом подпространства. Следовательно, решение исходной системы имеет вид,, гдеи- соответственно-я ия

координаты вектора .

Метод ортогонализации легко реализуется на ЭВМ и для решения системы уравнений требуетопераций умножения и деления иизвлечений квадратного корня. Однако удовлетворительные по точности результаты этот метод дает не для всех матриц. Это связано с неустойчивостью рекуррентного процесса (2.21), нарушающей ортогональность векторов. Чтобы избежать этого недостатка используется алгоритм Уилкинсона. В соответствии с этим алгоритмом векторывычисляются из соотношения

,

где ,.

При метод ортогонализации принимает обычный вид и может быть неустойчивым. Значительно лучшие результаты дает. Обычно же для достижения хорошей точности достаточно сделать две-три итерации.

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

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)

{

const float EPS=0.00000001;

float a[10][10],ac[10][10],f[10],u[10],u1[10],x[10],c[10];

n = StrToInt(Edit1->Text);

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

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

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

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

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

ac[i][j]=a[i][j];

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

c[i]=0;

u[i]=0;

u1[i]=0;

}

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

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

f[j] = StrToFloat(SG1->Cells[n][j]);

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

a[i][n]=-f[i];

a[n][i]=0;

}

a[n][n]=1;

float s;

float norm0=norm(*a);

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

a[0][k]=a[0][k]/norm0;

for(int k=0;k<=n-1;k++){

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

s=0;

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

s+=a[k+1][j]*a[i][j];

c[i]=-s;

}

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

s=0;

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

s+=c[i]*a[i][j];

u[j]=a[k+1][j]+s;

}

do{

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

c[j]=sk_pr(u,*(a+j));

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

s=0;

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

s+=c[i]*a[i][j];

u1[j]=u[j]-s;

}

s=0;

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

s+=pow(u1[i]-u[i],2);

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

u[j]=u1[j];

}

while(s>EPS);

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

norm0=norm(u);

a[k+1][j]=u[j]/norm0;

}

}

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

x[j]=a[n][j]/a[n][n];

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-1; i++){

f1[i]=0;

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

f1[i]+=ac[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);

//StringGrid1->RowCount = StrToInt(Edit1->Text);

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

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";

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

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

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

}

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