Алгоритмы приближения функций

Если функция y = f (x) задана, то тем самым любому допустимому значению x сопоставляется некоторое значение y.

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

Различают два основных способа выбора приближения функций : интерполяцию и аппроксимацию.

Использование интерполяции для решения задачи приближения функции

Постановка задачи в этом случае может быть сформулирована следующим образом.

Пусть функция f(X) задана таблицей ее значений Y0,Y1,...,Yn в точках X0,X1,...,Xn, называемых узлами интерполяции.

Задача интерполяции состоит в выборе такой функции

j (X), которая в узлах Xi принимала бы те же значения, что и f (X), то есть j (Xi) = Yi (i = 0,1,2,…,n ).

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

Степень интерполяционного многочлена будет равна n (так как число узлов n + 1).

Если значения функции f (X) получены экспериментально и содержат погрешности, то нет необходимости точного совпадения значений j (X) и f (X) в узловых точках. Более того, для практики можно пользоваться многочленом меньшей степени аналитической зависимостью другого вида, но график, которой проходил бы достаточно близко от точек (Xi, Yi) на всем протяжении интервала изменений.

Интерполяционная формула Ньютона

Введем понятие разделенных разностей.

Отношения

A1i = A (Xi, Xi+1) = [f(Xi+1) – f(Xi)] / (Xi+1 - Xi), i = 0,1,2,…, n –1

называются разделенными разностями первого порядка.

Отношения

A2i = A (Xi, Xi+1, Xi+2) = [A (Xi+1, Xi+2) - A (Xi, Xi+1)] / (Xi+2 - Xi) =

= (A1, i+1 – A1i)/ (Xi+2 - Xi) называются разделенными разностями второго порядка.

Разделенными разностями m-го порядка имеют вид:

Ami = (Am-1, i+1 – Am-1, i)/ (Xi+m - Xi).

Рассмотрим многочлен второй степени:

P2(X) = f(X0) + A10 (X0,X1)(X – X0) + A20 (X0,X1,X2)(X – X0)(X – X1).

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

P2 (X0) = f(X0);

P2 (X1) = f(X0) + A10(X0,X1)(X1 – X0) = f(X0) +[(f(X0)-f(X1))/(X1-X0)]*(X1-X0)=f(X0);

P2(X2) = f(X0) + A10(X0,X1)(X2 – X0) + A20(X0,X1, X 2)(X2 – X0)(X2 -- X1) = f(X2).

Таким образом, многочлен, построенный с помощью первых разделенных разностей первого и второго порядка, в точках X0,X1,X2 принимает те же значения, что и функция f(X).

Имеет место следующее утверждение:

Многочлен степени n Pn (X) = f (X0) + A10 (X0, X1)(X – X0) + A20 (X0, X1, X 2)(X – X0)(X - X1) + An0 (X0, X1, Xn)(X – X0)(X - X1)...(X – Xn-1) в точках X0,X1,...,Xn принимает те же значения, что f(X) .

Это выражение называется интерполяционным многочленом Ньютона

Для равноотстоящих узлов разделенные разности вычисляются по формуле:

Ami = (Am-1,i+1 – Am-1,i) / )(m * h) , где m = 1,2,…, n-1.

Тогда

Pn(X) = (...((An0(Z- Xn-1) + An-1,0) * (Z – Xn-2) + An-2,0)...)(X –X0) + A0,

где A0 = Y1.

Исходными данными для алгоритма вычисления интерполяционного многочлена Ньютона в случае равноотстоящих узлов являются следующие величины:

X – начальный узел интерполяции;

h – шаг интерполяции;

n – число узлов;

Z - значение аргумента, для которого вычисляется значение многочлена;

k – вспомогательная переменная (k = 0 при вычислении первого значения интерполяционного многочлена, k = 0 при последующих вычислениях )

Y(n) – массив значений функции (при k = 0) или массив разделенных разностей (при k = 1);

P – значение интерполяционного многочлена.

Замечание:

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

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

Интерполяционный многочлен Лагранжа.

Составление этого многочлена основано на одновременном введении в вычисления всех узлов интерполяции..

Коэффициентами Лагранжа называются функции:

(n)

Li(X) = ((X – X0) * (X – X1) * ...*(X – Xi-1) * (X - Xi+1) *...*(X–Xn)) / / ((X i – X0) * (Xi – X1) * ...*(Xi – Xi-1) * (Xi - Xi+1) *(Xi – Xn));

(n) (n)

Причем Li (Xi) = 1, Li (Xj) = 0 при i ¹ j.

Интерполяционным многочленом Лагранжа называется многочлен вида:

n (n)

Ln (X) = S Yi * Li (X).

i = 0

Недостатком интерполяции с использованием многочлена Лагранжа является тот факт, что при изменении числа узлов интерполяции n коэффициенты Лагранжа придётся вычислять заново.

Погрешность интерполяции

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

f (X) = Ln (X) + Rn (X), (1)

где Rn (X) – остаточный член, то есть погрешность интерполяции.

Получим некоторое выражение остаточного члена Rn (X) в предположении, что fÎ Cn+1[a,b] (принадлежит классу функций, определенных на [a,b] в (n + 1)-о й точке), где [a,b] - отрезок, содержащий узлы интерполяции Xi и точку X.

Будем искать остаточный член Rn (X) в следующем виде:

Rn (X) = Wn (X) * rn (X), (2)

где Wn (X) = (X – X0) * (X – X1) * ...*(X - Xn), (3)

rn (X) – некоторая функция, значения которой в узлах интерполяции Xi можно задавать какие угодно, так как

Rn (X) = 0 и Wn (Xi) = 0, i = 0,1,2,…, n.

Зафиксируем произвольное X Î [a,b] X ¹ Xi и рассмотрим следующую функцию от t:

j(t) = Ln (t) + Wn (t) * rn (X) – f(t). (4).

Эта функция на основании того, что Ln (Xi) = 1, обращается в ноль при t = Xi и t = X, то есть во всяком случае в (n + 2) –х точках отрезка [a,b] , на котором изменяется t.

По теореме Ролля, которая гласит:

если f (X) непрерывна на отрезке [a, b] и дифференцируема внутри этого интервала и f (a) = f(b), то найдется точка

c Î [a, b] такая, что f ’(c) = 0,

j’t обращается в ноль, по крайней мере, в (n+1) - Й точке интервала (a,b), j’’t равна нулю минимум в n точках этого интервала и так далее.

Таким образом, найдется хотя бы одна точка x Î (a,b), в которой

(n+1)

j t (x) = 0.

Отсюда и из (4), учитывая, что

(n+1) (n+1)

Ln (X) (x) = 0, Wn (x) = (n+1)!, получим

(n+1)

(n+1)! * rn (X) – f (x) =0.

Следовательно,

(n+1)

rn (X) = f (x) / (n +1)! и в соответствии с (1) и (2)

(n+1)

Rn (X) = Wn (X) * f (x) / (n + 1)!, (5)

(n+1)

f (X) = Ln(X) + (Wn (X) * f (x) / (n+1)!) ,где (6)

x = x(X) Î (a, b) – некоторая неизвестная точка.

Из равенства (6) вытекает оценка погрешности интерполяции в текущей точке XÎ [a, b]:

| f (X) – Ln (X)| £ ( Mn+1 / (n+1)! )* |Wn(X)| (7)

и оценка максимальной погрешности интерполяции на всем отрезке [a, b]:

max | f (X) – Ln (X)| £ ( Mn+1 / (n+1)! ) * max |Wn(X)|, где

[a, b] (n+1) [a, b]

Mn+1 = max f (X).

[a, b]

Рассмотрим пример.

Необходимо оценить погрешность приближения f(X) = ÖX в точке X =116 на отрезке [a, b], где a = 100, b = 144 с помощью интерполяционного многочлена Лагранжа L2(X) , построенного с узлами

X0 = 100, X1 = 121, X2 = 144.

Решение:

f ’(X) = ½ X^(- ½ ), f ’’(X) = - ¼ X^(-3/2), f ’’’ = 3/8 X^(-5/2).

M3 = max | f ’’’ (X)| = 3/8 * 100^(-5/2) = 3/8 * 10^(-5).

[a, b]

На основании (7):

|Ö116 - L2 (116)| £ (3/8) * 10^(-5)*(1/(3!)) * |(116 - 100) * (116 - 121) * *(116 - 144)| = 1/16 * 10^(-5) * 16 * 5 * 28 = 1.4 * 10^(-3).

В силу оценки (8):

Max |ÖX - L2 (X)| £ (10^(-5) / 16) * max |X -100| * |X -121| * |X -144| @

[a, b]

@ 2.5 * 10^(-3).

Аппроксимация при решении задачи приближения функций

Задача интерполяции состояла в построении многочлена Pn (X) степени n, значения которого в (n+1) точках интерполяции соответствовали рассматриваемой функции f (X), то есть

Pn (Xi) = f (Xi) (8)

С увеличением числа узлов интерполяции n трудности построения такого многочлена возрастают.

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

Отсюда следует задача об определении многочлена более низкой степени m (m << n)

m

Pm (X) = amX + ... + a1X + a0, «расстояние» между которым и функцией f (X) в некотором смысле минимально.

Такой многочлен называется аппроксимирующим. Он сглаживает локальные особенности экспериментальной таблицы и отражает общее поведение функции f (X).

Метод наименьших квадратов

Пусть задана таблица значений функции Y = f (X) в виде пар значений {Xi, yi}, i = 1,2,…, n.

Аппроксимирующуя функцию будем разыскивать в виде:

Y = j (X, a0, a1, a2,..., am), (1)

где a0, a1, a2,..., am – подлежащие определению параметры.

Вычислим значения функции (1) при X = Xi:

Yi = j (Xi, a0, a1, a2,..., am), i = 1,2,…, n.

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

n 2 n 2

S(a0,a1,a2,...,am)=S(Yi-yi)=S(j(Xi,a0,a1,a2,...,am)-yi). (2)

i=1 i=1

Параметры a0, a1, a2,..., am находятся из условия минимума ошибки аппроксимации (2):

n

¶S/¶ ak = 2 S [j ( Xi, a0, a1,..., am)– yi] (¶/¶ ak) j (Xi, a0, a1,...,am) = 0 (3).

i =1

В результате получаем систему из (m+1) уравнений для определения (m+1) неизвестных параметров.

R0a0 + R1a1 + ... +Rmam = B0

R1a0 + R2a1 + ... +Rm+1am = B1

....................................................

Rma0 + Rm+1a1 + ... +R2mam = Bm, где

n k n j

Rk = S Xi / n , k = 0,1,…, 2m, Bj = S yi Xi , j = 0,1,…, n

i=1 i=1

Алгоритм аппроксимации функции методом наименьших квадратов.

Исходные данные:

- n – количество точек, где известны значения функции;

- m – степень аппроксимирующего многочлена;

- X, Y- массивы табличных значений аргументов и функции;

- A – массив коэффициентов системы линейных уравнений;

- R, B – массивы коэффициентов и свободных членов.

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

clip_image002clip_image003clip_image004 n

clip_image005D = S (f (Xi) – Pn (Xi))^2 / (n + 1)

i=1

Сглаживание наблюдений

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

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

Некоторые формулы для сглаживания:

`fi = 1/21 (- 2fi-3 + 3fi-2 + 6fi-1 + 7fi + 6fi+1 + 3fi+2 –2fi+3) (1)

k

`fi = (1 /(n +1)) * (fi + S (f i-j + f i+j ), где n = 21 (2)

i=1

`fi=1/35(-3fi-2+12fi-1 + 17fi + 12fi+1 – 3fi+2) (3)

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

1) Если наблюдаемые значения содержат случайные ошибки со свойствами:

M [hi] = 0, i = 0,1,…,n, где M [hi] – математическое ожидание ошибки в i- ом наблюдении

0, i ¹ j

и дисперсия не зависит от I,то есть M [h, hj] =

s^2, i = j,

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

В частности, при использовании (1)

`s^2 = 1/3 s^2 (`s = (1 / 3) * s), где

`s - среднеквадратичные ошибки сглаженных значений,

s - среднеквадратичные ошибки наблюдений.

2) Ошибки в сглаженных значениях `fi становятся сильно коррелированными.

3) Если функция не является алгебраическим многочленом степени n , то происходит “размазывание” ее графика.

Повторное сглаживание делать не рекомендуется.

Предлагаю ознакомиться с аналогичными статьями: