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;
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,'ТАБЛИЦА ЗНАЧЕНИЙ');
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;
- Содержание
- 3. Интерполялия, экстрополяция, аппроксимация, сглаживание 5
- 3. Интерполялия, экстрополяция, аппроксимация, сглаживание
- 3.1. Введение
- 3.2. Интерполяция
- 3.2.1. Полиномиальная интерполяция
- Аппроксимационная теорема Вейерштрасса.
- 3.2.2. Вычисление значений многочлена. Схема Горнера
- 3.2.3. Линейная интерполяция
- 3.2.4. Квадратичная интерполяция
- 3.2.5. Построение других базисных функций
- 3.2.6. Многочлены Тейлора
- 3.2.7. Лагранжева интерполяция
- I, j, n : Integer;
- 3.2.8. Ошибки полиномиальной интерполяции
- 3.2.9. Кусочно-линейная интерполяция
- Var X,y : Array[0..N] of Real;
- I,j : Integer;
- Var f:Real;
- 3.2.10. Кусочно-кубическая интерполяция
- 3.2.11. Эрмитов кубический интерполянт
- 3.2.12. Кубические сплайны
- Var r, s, l : Vect;
- Var l, I, j : Integer;
- 1 : Begin
- 0 : Begin
- Var XX:RealType;
- 3.2.13. Кривые Безье. Сплайны
- 3.2.14. Итерационный способ вычисления интерполяционного полинома (способ Эйткена)
- 3.2.15. Интерполяционный многочлен Ньютона
- 3.2.16. Интерполяционный многочлен Гаусса
- 3.2.17. Интерполяционный многочлен Стирлинга
- 3.2.18. Интерполяционный многочлен Эверетта
- 3.3. Аппроксимация данных методом наименьших квадратов
- 3.3.1. Аппроксимация данных методом наименьших квадратов
- 3.3.2. Аппроксимация данных с другими нормами
- 3.3.3. Аппроксимация данных многочленом заданной степени
- Var X,y:array[1..Nmax] of real;
- I,n:integer;
- Литература
- Простейшие способы интерполяции
- Интерполяционные полиномы
- Сплайн-интерполяция
- Тригонометрическая интерполяция
- Неклассические методы интерполяции
- Реконструкция функций
- Всюду гладкая интерполяция