Лекция_3(Интерполяция)

3.3.3. Аппроксимация данных многочленом заданной степени

Пусть, например, функция задана в виде набора значений . Задача состоит в аппроксимации неизвестной функциональной зависимости между и многочленом заданной степени :

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

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

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

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

В данном случае заданную табличную функцию требуется аппроксимировать многочленом -й степени . В этом случае нормальная система уравнений относительно коэффициентов и имеет вид

Эту систему можно представить в матричном виде:

Здесь как обозначены суммы соответствующих степеней -координат заданных точек:

а , , используются для обозначения сумм произведений -координат заданных точек на соответствующие степени их -координат:

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

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

Кроме абсолютной и относительной погрешности необходимо также определить среднее квадратическое отклонение (СКО) по формуле

CKО характеризует разброс параметров, или ширину интервала рассеяния, в котором находится большинство точек исходной зависимости.

Пример 1.

{$N+,E+,F+}

PROGRAM MetNKv;

{Аппроксимация функции по методу наимньших квадрвтов}

{Полиномом заданной степени}

Uses CRT;

Const

N=21; {Количество заданных точек кривой намагничивания}

{NN=4;} {Степень полинома}

NameInt = 'mnk.inp'; {Имя файла с данными}

NameOut = 'mnk4-pu.out';} {Имя выходного файла

Type

RealType = Extended; {базовый вещественный тип}

IntType = Integer; {базовый целый тип}

Vec = Array [1..N] of RealType;

TypeS = Array [1..28] of RealType;

TypeA = Array [1..29,1..29] of RealType;

TypeB = Array [1..29] of RealType;

TypeBool = Boolean;

{Глобальные переменные}

Var

S : TypeS;

h,SS,OS,AS : RealType;

Y,X : Vec;

B,P : TypeB;

A : TypeA;

Flag : TypeBool;

F1, F2 : Text;

i, j, k, NN : IntType;

Procedure KoeffA(NN: IntType; Var S: TypeS);

{Вычисление коэффициентов матрицы A}

Var

i, j : IntType;

BEGIN

For j:=1 to 2*NN do

begin

S[j]:=0;

For i:=1 to N do

S[j]:= S[j] + exp(j*ln(X[i]));

end;

END; {KoeffA}

Procedure KoeffB(NN,N: IntType; Y,X: Vec;Var B: TypeB);

{Вычисление правой части системы уравнений B}

Var

i, j : IntType;

BEGIN

For j:=1 to NN+1 do

begin

B[j]:=0;

For i:=1 to N do

B[j]:= B[j] + Y[i]*exp((j-1)*ln(X[i]));

end;

END;{KoeffB}

Procedure GaussKlass(N: IntType; A: TypeA; B: TypeB; Var X: TypeB);

{Решение системы уравнений классическим методом Гаусса}

{Прямой ход исключения переменных}

Var

i, j, k : IntType;

BEGIN

For i:= 1 to N-1 do

begin

For j:= i+1 to N do

begin

A[j,i]:= -A[j,i]/A[i,i];

For k:= i+1 to N do

A[j,k]:= A[j,k] + A[j,i]*A[i,k];

B[j]:= B[j] + A[j,i]*B[i];

end;

{Обратный ход: последовательное нахождение корней}

X[N]:= B[N]/A[N,N];

For i:= N downto 1 do

begin

h:= B[i];

For j:= i+1 to N do h:= h - A[i,j]*X[j];

X[i]:= h/A[i,i];

end;

END; {GaussKlass}

{Головная программа}

BEGIN

ClrScr;

WriteLn;

Write('Введите степень полинома, NN = ');

ReadLn(NN);

Assign(F1,NameInt);

Reset(F1);

For i:=1 to N do Read(F1,X[i]);

For i:=1 to N do Read(F1,Y[i]);

Assign(F2,NameOut);

Rewrite(F2);

WriteLn(F2);

{Вычисление коэффициентов матрицы A}

KoeffA(NN,S);

WriteLn(F2,'МАТРИЦА СИСТЕМЫ УРАВНЕНИЙ A:');

{Заполнение коэффициентов матрицы}

For i:=1 to NN+1 do

begin

For j:=1 to NN+1 do A[i,j]:=S[i+(j-2)];

A[1,1]:= N;

end;

For i:=1 to NN+1 do

begin

For j:=1 to NN+1 do Write(F2,A[i,j]);

WriteLn(F2);

end;

WriteLn(F2);

WriteLn(F2,'ВЕКТОР-СТОЛБЕЦ ПРАВОЙ ЧАСТИ B:');

{Вычисление правой части системы уравнений B}

KoeffB(NN,N,Y,X,B);

For j:=1 to NN+1 do WriteLn(F2, B[j]);

WriteLn(F2);

GaussKlass(NN+1,A,B,P);

WriteLn(F2);

WriteLn(F2,'КОЭФФИЦИЕНТЫ ПОЛИНОМА ',NN:2,' СТЕПЕНИ:');

j:=0;

For i:=1 to NN+1 do

begin

WriteLn(F2,'P',j,'=',P[i],' ');

j:= j+1;

end;

WriteLn(F2);

WriteLn(F2,'ТАБЛИЦА ЗНАЧЕНИЙ');

WriteLn(F2,' N ',' ЗНАЧ.АРГ.',' ИСТ.ЗН.ФУНК.', ' АПП.ЗН.ФУНК.',' АБС.ПОГР.', ' ОТН. ПОГР.');

WriteLn(F2,' ',' X ',' Y ', ' Y* ',' (Y-Y*) ', ' (Z) ');

SS:=0;

For i:=1 to N do

begin

{Значения аппроксимированной функции}

h:=0;

For j:=1 to NN+1 do

h:= h + P[j]*exp((j-1)*ln(X[i]));

{Абсолютная погрешность}

AS:= Abs(Y[i] - h);

{Относительная погрешность}

OS:= Abs(AS/Y[i])*100;

SS:= SS + AS*AS;

WriteLn(F2,i:3,' ',X[i]:8:2,' ',Y[i]:9:3,' ',h:11:4,' ',AS:9:5,' ',OS:9:5);

end;

{Среднеквадратическое отклонение}

SS:=Sqrt(SS/(N-1));

WriteLn(F2);

WriteLn(F2,'СРЕДНЕКВАДРАТИЧЕСКОЕ ОТКЛОНЕНИЕ:',SS);

Close(F1);

Close(F2);

ClrScr;{Очищает экран, модуль Crt}

END.

Пример 2.

const nmax=25;

Содержание