4. Реализация алгоритма Гаусса в пакете Mathcad

Для реализации алгоритма Гаусса в пакете Mathcad необходимо знать некоторые функции определения матриц и операции с блоками матриц:

matrix(m,n,f) – создает и заполняет матрицу размерности , элемент которой, расположенный в i-ой строке, j-м столбце, равен значению функции ;

diag(v) – создает диагональную матрицу, элементы главной диагонали которой хранятся в векторе v;

identity(n) – создает единичную матрицу порядка n;

augment(A,B) – формирует матрицу, в первых столбцах которой содержится матрица A, а в последних – матрица B (матрица A и Bдолжны иметь одинаковое число строк);

stak(A,B) – формирует матрицу, в первых строках которой содержится матрица A, а в последних – матрица B(матрица A и Bдолжны иметь одинаковое число столбцов);

submatrix(A,ir,jr,ic,jc) – формирует матрицу, которая является блоком матрицы А, расположенным в блоках с ir по jr и в столбцах с ic по jc, ir ≤ jr, ic ≤ jc.

Номер первой строки (столбца) матрицы или первой компоненты вектора хранится в Mathcad в переменной ORIGIN. По умолчанию в Mathcad координаты векторов, столбцы и строки матрицы нумеруются начиная с 0 (ORIGIN:=0). Поскольку в математической записи чаще используется нумерация с 1, здесь и в дальнейшем перед началом работы с матрицами будем определять значение переменной ORIGIN:=1.

Приведем пример решения системы методом Гаусса в Mathcad:

5. Реализация алгоритма Гаусса на языке Turbo Pascal

Введем двумерный массив

const n_max=20;

type dim2=array[1..n_max, 1..n_max] of real;

и запишем процедуру решения СЛАУ методом исключения Гаусса:

procedure Gauss_01( var a:dim2; n:integer);

var i,j,k,np,kp :integer; s:real;

begin

np:=n+1;

for k:=1 to n do { прямой ход}

begin kp:=k+1;

for j:=kp to np do { нормировка на единицу диагонального элемента }

a[k,j]:=a[k,j]/a[k,k];

for i:=kp to np do { последовательное исключение всех элементов в к-ом столбце}

for j:=kp to np do a[i,j]:=a[i,j]-a[i,k]*a[k,j];

end;

for i:=n-1 downto 1 do {обратный ход}

begin s:=a[i,n+1];

for j:=i+1 to n do s:=s+a[i,j]*a[j,np];

a[i,np];=s

end

end;

Обратите внимание, где размещены результаты вычислений. Прежде всего, необходимо отладить программу и произвести тестовые расчеты. Для тестирования можно рассмотреть различные способы формирования матрицы . Приведем, простую процедуру:

procedure Data (var a:dim2; var n:integer);

var i, j:integer;

begin n:=4;

for i:=1 to n do begin

for j:=1 to n do

if i=j

then a[i,j]:=4.0

else a[i,j]:=1.0;

a[i, n+1]:=i*1.0

end

end;

Таким образом, формируется исходная система уравнений:

, (1.7)

результатом вычислений которой является вектор

Прямой ход метода Гаусса строится на трех вложенных циклах (см. листинг подпрограммы), следовательно общее число арифметических операций, требуемых для приведения матрицы к верхнетреугольному виду, будет пропорционально Эту оценку принято записывать как Обратный ход метода исключения Гаусса требует на порядок меньше операций Для сравнения отметим, что в методе Крамера на вычисление только одного определителя требуется порядка n! операций. Так, при n=100 по формуле Стирлинга получим Эта вычислительная проблема, связанная с использованием метода Крамера, названа математиками «проклятием размерности». Алгоритм Гаусса требует всего порядка операций. Он прост, нагляден и является одним из лучших алгоритмов решения линейных систем.