Аппроксимация МНК (стр. 2 из 4)

Обычно указывают класс функций (например, множество линейных, степенных, показательных и т.п.) из которого выбирается функция

, и далее определяются наилучшие значения параметров.

Если в эмпирическую формулу (1) подставить исходные

, то получим теоретические значения

, где

.

Разности

называются отклонениями и представляют собой расстояния по вертикали от точек

до графика эмпирической функции.

Согласно методу наименьших квадратов наилучшими коэффициентами

считаются те, для которых сумма квадратов отклонений найденной эмпирической функции от заданных значений функции

(2)

будет минимальной.

Поясним геометрический смысл метода наименьших квадралтов.

Каждая пара чисел

из исходной таблицы определяет точку

на плоскости

. Используя формулу (1) при различных значениях коэффициентов

можно построить ряд кривых, которые являются графиками функции (1). Задача состоит в определении коэффициентов

таким образом, чтобы сумма квадратов расстояний по вертикали от точек

до графика функции (1) была наименьшей.

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

Если неизвестен характер зависимости между данными величинами x и y, то вид эмпирической зависимости является произвольным. Предпочтение отдается простым формулам, обладающим хорошей точностью. Удачный выбор эмпирической формулы в значительной мере зависит от знаний исследователя в предметной области, используя которые он может указать класс функций из теоретических соображений. Большое значение имеет изображение полученных данных в декартовых или в специальных системах координат (полулогарифмической, логарифмической и т.д.). По положению точек можно примерно угадать общий вид зависимости путем установления сходства между построенным графиком и образцами известных кривых.

Определение наилучших коэффициентов

, входящих в эмпирическую формулу производят хорошо известными аналитическими методами.

Для того, чтобы найти набор коэффициентов

, которые доставляют минимум функции S , определяемой формулой (2), используем необходимое условие экстремума функции нескольких переменных - равенство нулю частных производных. В результате получим нормальную систему для определения коэффициентов

:

(3)

Таким образом, нахождение коэффициентов

сводится к решению системы (3).

Эта система упрощается, если эмпирическая формула (1) линейна относительно параметров

, тогда система (3) - будет линейной.

Конкретный вид системы (3) зависит от того, из какого класса эмпирических формул мы ищем зависимость (1). В случае линейной зависимости

система (3) примет вид:

(4)

Эта линейная система может быть решена любым известным методом (методом Гаусса, простых итераций, формулами Крамера).

В случае квадратичной зависимости

система (.3) примет вид:

(5)

Метод Гаусса решения систем линейных алгебраических уравнений

Метод Гаусса состоит в последовательном исключении неизвестных до тех пор, пока не останется одно уравнение с одним неизвестным.

При этом матрица СЛАУ приводится треугольному виду, где ниже главной диагонали располагаются только нули.

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

Имеем Ax=b (6), где A=[a_ij] - матрица размерности n*n, det A

0, b=(b₁, b₂, …, b_n)^T.

В предположении, что a₁₁

0, первое уравнение системы (6):

делим на коэффициент a₁₁, в результате получаем уравнение

.

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

Первое неизвестное оказалось исключенным из всех уравнений, кроме первого. Далее предполагаем, что 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];