Первое  неизвестное оказалось исключенным из всех уравнений, кроме первого. Далее предполагаем, что a¹₂₂ 0, делим второе уравнение на a¹₂₂ и исключаем неизвестное x₂ из всех уравнений, начиная со второго, и т.д. В результате последовательного исключения неизвестных система уравнений преобразуется в систему уравнений с треугольной матрицей:

                                            (7)

     Совокупность  проведенных действий называется прямым ходом метода Гаусса.

     Из n-го уравнения системы (6) определяем x_n, из (n-1)-го - x_n-1 и т.д. до x₁. Совокупность таких действий называется обратным ходом метода Гаусса.

     Реализация  прямого хода требует арифметических операций, а обратного -  арифметических операций.

 

Блок-схема алгоритма  программы на языке Pascal  

 

Практическая  часть

Текст программы на языке Pascal, перечень использованных в программе идентификаторов и полученные результаты

     Листинг программы:

program MultiChain;

const

  xy: array [1..40,1..2] of real = ((0.15,    7.26),      (-0.01,   7.95),

      (0.46,    8.77),      (1.71,    9.72),      (3.94,    10.78),      (7.40,    12.45),      (12.41,   14.74),      (19.47,   18.18),      (29.21,   23.36),      (42.49,   30.99),      (60.43,   41.97),      (84.42,   57.38),      (116.18,  78.45),      (157.80,  106.65),      (211.78,  143.64),      (281.04,  191.28),      (369.00,  251.71),      (479.59,  327.27),      (617.30,  420.58),      (787.21,  534.52),      (995.04,  672.24),      1247.17, 837.2),      (1550.71, 1033.14),      (1913.5,  1264.13),     (2344.18, 1534.54),      (2852.21, 1849.11),      (3447.91, 2212.91),    (4142.52, 2631.36),      (4948.19, 3110.26),      (5878.09, 3665.8),     (6946.36, 4274.55),      (8768.24, 4973.5),      (9560.03, 5760.63),     (11139.19,6641.98),      (12924.34,7627.6),      (14935.29,8725.62),     (17193.14,9945.22),      (19720.24,11296.05),      (22540.29,12788.24),      (25678.32,14432.44));

var

  matr: array [1..20,1..20] of real; { матрица решения }

  sum: array [1..21] of real;      { коэффициенты матрицы }

  B: array [1..20] of real;       { временный вектор }

  A: array [1..20] of real;

  n: integer;

  m: integer; 

procedure FillMatr;

var

  i,j,z: integer;

  ss: real;

begin

  for i:=1 to m+1 do sum[i]:=0;   { инициализируем переменные }

  for i:=1 to 20 do B[i]:=0; 

   

  for z:=1 to 40 do begin

    ss:=1;

    for i:=1 to m+1 do begin      { подсчитываем коэффициенты - степени переменных }

      sum[i] := sum[i] + ss;

      ss := ss * xy[z,1];

    end;

    ss:=xy[z,2];

    for i:=1 to n+1 do begin

      b[i] := b[i] + ss;

      ss := ss * xy[z,1];

    end;

  end;

  for i:=1 to n+1 do matr[n+2,i] := b[i];

 

  for i:=1 to n+1 do

    for j:=1 to n+1 do

    begin

      matr[i,j] := sum[i + (j-1)];

    end;

 

end; 

procedure TryFind;

var

  i,j,z:integer;

  Q: real;

begin

  {1 step}

 

  for z:=1 to n do begin

    Q := matr[z,z];

    for j:=1 to n+2 do

      matr[j,z]:= matr[j,z] / Q;

    for i:=z+1 to n+1 do begin

      Q := (-1) / matr[z,i];

      for j:=1 to n+2 do

        matr[j,i]:= Q * matr[j,i] + matr[j,z];

    end;

  

  end;

  Q := matr[n+1,n+1];

  for j:=n+1 to n+2 do

    matr[j,n+1]:= matr[j,n+1] / Q;

   

  {2 step}

  for z:=n downto 1 do begin

    Q := matr[n+2,z+1];

    for i:=z downto 1 do begin

      matr[n+2,i]:=matr[n+2,i] - matr[z+1,i] * Q;

      matr[z+1,i]:=0;

    end;

      end;

 

  Writeln;

  Write('F(x)=');

  for i:=1 to n+1 do begin

    A[i]:=matr[n+2,i];

    if i=1

      then Write('(',A[i]:0:4,')')

      else Write('(',A[i]:0:4,')*x^',i-1);

    if i<>n+1

      then Write('+')

      else Writeln

  end;

 

end; 

function func(x:real):real;

{ функция для  вычисления значения искомой  функции по вычисленному полиному }

var

  i: integer;

  c,ff: real;

begin

  c:=1;

  ff:=0;

  for i:=1 to n+1 do begin

    ff := ff + A[i] * c;

    c := c * x;

  end;

  func := ff;

end; 

procedure TryCheck;

var

  i:Integer;

  x,y,f,

  delta:real;

begin

  delta := 0;

  for i:=1 to 40 do begin

    x:=xy[i,1]; y:=xy[i,2]; f:=func(x);  

    delta := delta + sqr(y - f);

  end;

  delta := sqrt(1/40*delta);

  Writeln('delta = ',delta:0:6);

end; 

begin

  n := 2; { степень полинома }

  writeln('Степень  полинома n=',n);

  m := 2*n; { количество  коэффициентов }

  FillMatr;

  TryFind;

  TryCheck;

  writeln;

    n := 4; { степень полинома }

  writeln('Степень  полинома n=',n);

  m := 2*n; { количество коэффициентов }

  FillMatr;

  TryFind;

  TryCheck;

  writeln;

    n := 5; { степень полинома }

  writeln('Степень  полинома n=',n);

  m := 2*n; { количество  коэффициентов }

  FillMatr;

  TryFind;

  TryCheck;

end.

     

     Перечень использованных в программе идентификаторов и их описание:

xy: array [1..40,1..2] of real - массив исходных значений;

matr: array [1..20,1..20] of real - основная матрица для нахождения решения;

sum: array [1..21] of real - массив коэффициентов для матрицы поиска;

B: array [1..20] of real - массив вектора свободных членов;

A: array [1..20] of real - массив для искомых коэффициентов полинома;

n: intege r - переменная для размерности матрицы поиска;

m: integer - переменная для размерность массива коэффициентов.

     Переменные, встречающиеся в большинстве  процедур:

i, j, z: integer - переменные  итераций для перебора элементов массива.

      Переменные  в процедурах:

      Процедура FillMatr

ss: real - переменная для подсчета коэффициентов матрицы;

      Процедура TryFind

Q: real - переменная для вычисления множителя;

     Функция func

c: real - переменная для вычисления степени переменной x(от 0-ой до n-ой);

ff: real - переменная для вычисления значения функции-полинома;

     Процедура TryCheck

x: real - заданное значение независимой переменной;

y: real - заданное значения неизвестной функции в точке x;

f: real - переменная для значения вычисленной функции-полинома в точке x;

delta: real - среднеквадратичное отклонение.

     Программа, разработанная  в среде PascalABC дает следующие результаты:

 

Описание  решения задачи в  среде MathCAD

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

      

     Функция regress является вспомогательной, она подготавливает данные, необходимые для работы функции interp. Вектор vs содержит, в том числе, и коэффициенты полинома

      

      Функция interp возвращает значение полинома в точке z. Определив новые функции f2, f4 и f5, получим возможность находить значение полинома в любой заданной точке.

     

     

     

     Коэффициенты задаем следующим образом:

     

     

     

     Коэффициенты  имеют вид:

     

     

     

     Среднеквадратичные отклонения для каждого из трех случаев:

    

    

    

    

    

 

Результаты вычислений и их анализ

 

Результаты вычислений и их анализ

     Представим результаты  в виде сравнительных таблиц для применявшихся сред вычисления.

     Для степени полинома n=2 получены следующие коэффициенты и среднеквадратичное отклонение (с точностью 0,001):

     Итак, коэффициенты полинома второй степени, полученные в двух средах вычислений фактически не отличаются (отличается лишь А₀ на 0,002). Получен следующий полином второй степени, аппроксимирующий исходную табулируемую функцию: