Матрица зависит от параметра α и формируется следующим образом: pii= α, pij=1 для i≠j.

Задание

Вариант 2

Матрица зависит от параметра α и формируется следующим образом: p_ii = α, p_ij =1 для i ≠ j.

Решить систему y= b, где b =(1, 1, …, 1) ^T, используя процедуры DECOMP и SOLVE. Порядок матрицы равен 6. Вычислительный эксперимент состоит в исследовании априорной и апостериорной оценок качества матриц при использовании различных форматов представления данных. Априорной оценкой является стандартное число обусловленности программы DECOMP, а апостериорной относительная погрешность сравниваемых решений

α =0,9, 0,99, 0,999 и т. д. вплоть до нештатной реакции компьютера.

1. Математическая постановка задачи

Процедура DECOMP находит априорную оценку, а SOLVE решает систему y= b. Результатом процедурs DECOMP является априорные оценки матриц а SOLVE – решения x ₁, x ₂ систем y= b и y= b. Апостериорная относительная погрешность вычисляется по следующей формуле:

2. Схема алгоритма главной программы

3. Текст программы

program Lab3;

{ Use modules}

uses FMM, printer, crt;

procedure GenerateMatrixP(n: Integer; var P: floatmatrix; param: float);

var i, j: Integer;

begin

for i:= 1 to n do

begin

for j:= 1 to n do P[ i, j]:= 1;

P[ i, i]:= param;

end;

end;

procedure GenerateVectorb(n: Integer; var b: floatvector);

var i: Integer;

begin

for i:= 1 to n do b[ i]:= 1;

end;

function GenerateParam(i: Integer): float;

begin

GenerateParam:= 1.0 - powri(0.1, i);

end;

function GetAposteriorEstimation(n: Integer; X1, X2: floatvector): float;

var

i: Integer;

t1, t2: float;

begin

t1:= 0; t2:= 0;

for i:= 1 to n do

begin

t1:= t1 + powri(X1[ i] - X2[ i], 2);

t2:= t2 + powri(X1[ i], 2);

end;

GetAposteriorEstimation:= sqrt(t1 / t2);

end;

procedure MultMatrixMatrix(n: Integer; P, b: floatmatrix; var Res: floatmatrix);

var

i, j, k: Integer;

begin

for i:= 1 to n do

begin

for j:= 1 to n do

begin

Res[ i, j]:= 0.0;

for k:= 1 to n do Res[ i, j]:= Res[ i, j] + P[ i, k] * b[ k, j];

end;

end;

end;

procedure MultMatrixVector(n: Integer; P: floatmatrix; b: floatvector; var Res: floatvector);

var

i, j: Integer;

begin

for i:= 1 to n do

begin

Res[ i]:= 0.0;

for j:= 1 to n do Res[ i]:= Res[ i] + P[ i, j] * b[ j];

end;

end;

procedure PrintMatrix(n: Integer; P: floatmatrix; var pFile: TEXT);

var

i, j: Integer;

begin

for i:= 1 to n do

begin

for j:= 1 to n do write(pFile, P[ i, j]:20:7);

writeln(pFile);

end;

end;

procedure PrintVector(n: Integer; P: floatvector; var pFile: TEXT);

var

i: Integer;

begin

for i:= 1 to n do writeln(pFile, P[ i]:20:7);

end;

procedure Transposition(n: Integer; P: floatmatrix; var Pt: floatmatrix);

var

i, j: Integer;

begin

for i:= 1 to n do

for j:= 1 to n do Pt[ i, j]:= P[ j, i];

end;

{ Constants}

const n = 6; { Dimension of matrix P and b}

nP = 7; { Number of matrix P}

nPt = 5;

Var

P, P1, Pt: floatmatrix;

b1, W, b: floatvector;

cond: float;

Ip: ivector;

param: float;

i, j, nRes: Integer;

pFile: TEXT;

ResP: array [1..nA, 1..n] of float;

ResPt: array [1..nAt, 1..n] of float;

Begin

{ Clear screen}

clrscr;

{ Open file to write result data for P*y=b equation}

assign (pFile, 'resultPXB.txt');

rewrite(pFile);

for i:= 1 to nP do

begin

{ Get param for matrix P}

param:= GenerateParam(i);

{ Generate initial data}

GenerateMatrixP(n, P, param);

GenerateVectorb(n, b);

{ Print initial data "param, P, b"}

write(pFile, '-----------------------------------------------'); writeln(pFile);

write(pFile, 'Parameter alpha = '); writeln(pFile, param);

writeln(pFile, 'Matrix P:'); PrintMatrix(n, P, pFile);

writeln(pFile, 'Matrix b:'); PrintVector(n, b, pFile);

{ Calculate solution of equation P*y=b}

decomp(n, P, cond, Ip, W); { Calculate gaussian matrix}

solve(n, P, b, Ip); { Calculate solution, result in b}

{ Print result of solution P*y=b}

writeln(pFile, 'Result P*y=b:'); PrintVector(n, b, pFile);

for j:= 1 to n do ResP[ i, j]:= b[ j];

{ Aprior estimation of quality}

write(pFile, 'Aprior estimation (cond) = '); writeln(pFile, cond);

end;

close(pFile);

{Open file to write result data for Pt*P*y=Pt*b equation}

assign (pFile, 'resultPtPXPtb.txt');

rewrite(pFile);

for i:= 1 to nPt do

begin

{ Get param for matrix P}

param:= GenerateParam(i);

{ Generate initial data}

GenerateMatrixA(n, P1, param);

GenerateVectorB(n, b1);

Transposition(n, P1, Pt);

MultMatrixMatrix(n, Pt, P1, P); { Get Pt*P}

MultMatrixVector(n, Pt, b1, b); { Get Pt*P }

{ Print initial data "param, P, b"}

write(pFile, '-----------------------------------------------'); writeln(pFile);

write(pFile, 'Parameter an = '); writeln(pFile, param);

writeln(pFile, 'Matrix Pt*b:'); PrintMatrix(n, P, pFile);

writeln(pFile, 'Matrix Pt*b:'); PrintVector(n, P, pFile);

{ Calculate solution of equation Pt*P*y=Pt*b}

decomp(n, P, cond, Ip, W); { Calculate gaussian matrix}

solve(n, P, b, Ip); { Calculate solution, result in }

{ Print result of solution Pt*P*y=Pt*b }

writeln(pFile, 'Result Pt*P*y=Pt*b:'); PrintVector(n, b, pFile);

for j:= 1 to n do ResPt[ i, j]:= b[ j];

{ Aprior estimation of quality}

write(pFile, 'Aprior estimation (cond) = '); writeln(pFile, cond);

end;

close(pFile);

assign (pFile, 'est.txt');

rewrite(pFile);

writeln(pFile, 'Aposterior estimation:');

if nA > nPt then nRes:= nPt

else nRes:= nP;

for i:= 1 to nRes do

begin

for j:= 1 to n do b [ j]:= ResP[i, j];

for j:= 1 to n do b1[ j]:= ResPt[i, j];

writeln(pFile, GetAposteriorEstimation(n, b, b1));

end;

close(pFile);

end.

4. Результаты в виде таблиц и графиков

5. Интерпретация результата

С увеличение точности параметра α, увеличиваются априорные оценки и апостериорная относительная погрешность, сильнее друг от друга отличаются решения x ₁ и x ₂.

Дата добавления: 2015-08-29