Интерполяция, сглаживание и аппроксимация

Программа PL10 вычисляет значение функции y(x) в заданной точке x с помощью интерполяционного процесса Эйткена-Лагранжа. Функция y(x) задана таблично массивами Xi и Yi, i=1, 2, . n. Массив Xi должен быть упорядочен в возрастающем порядке, т.е. для всех его элементов xi<xi+1. Максимальное число узлов m, используемых для интерполяции, задается при вызове программы. Процесс интерполяции построен по итерационной схеме. В качестве первого приближения используется значение функции в ближайшем узле интерполяции. Для получения более точных приближений используются последовательно вычисляемые значения многочлена Лагранжа L (k) (x), степень которого возрастает на каждом шаге итерации:

где ω(k)(x) = (x-x1)(x-x2) . (x-xk+1), k=1, 2, . n-1. Итерационный процесс продолжается до получения заданной точности интерполяции. Вычисления проводятся по треугольной интерполяционной схеме Эйткена: x1 y1 x2 y2 L12 x3 y3 L13 L123 x4 y4 L14 L124 L1234 . . . . . . . . . . . . . . . . . . xm ym L1m L12m L123m . L123. m Каждый многочлен Лагранжа более высокой степени получается из двух предыдущих многочленов по рекуррентным формулам:

Элементы L12. k нисходящей диагонали этой схемы представляют собой интерполяционные многочлены Лагранжа L (k) (x) последовательно возрастающих степеней (k=2, . m). Для определения значений многочлена L (k) (x) при росте k в вычислительный процесс вовлекаются все более удаленные от x узлы интерполяции. Поэтому значения xk и yk (k=1, 2, . m) предварительно упорядочиваются таким образом, чтобы модули разностей |xk - x| возрастали при росте индекса k: |xk - x| ≤ |xk+1 - x| .

CALL PL10(X, Y, XX, M, Eps, YY, Error)

X, Y, XX, M, Eps - входные параметры; YY, Error - выходные параметры;

Real X(1:N), Y(1:N) - массивы, содержащие значения аргумента и функции. Массив X должен быть упорядочен в возрастающем порядке; Real XX - значение аргумента, в котором требуется вычислить значение функции; Integer M - число узлов, которое разрешается использовать при интерполяции, 1 PL12 вычисляет значение функции y(x) в заданной точке x с помощью интерполяционного процесса Эйткена-Эрмита. Функция y(x) задана таблично массивами Xi и Yi, i=1, 2, . n. Массив Xi должен быть упорядочен в возрастающем порядке, т.е. для всех его элементов xi<xi+1. Максимальное число узлов m, используемых для интерполяции, задается при вызове программы. Процесс интерполяции построен по итерационной схеме. В качестве первого приближения используется интерполяционный многочлен Эрмита первой степени. На каждой новой итерации степень многочлена увеличивается на 1. Итерационный процесс продолжается до получения заданной точности. Для построения многочленов Эрмита нечетной степени 2k-1 используются значения функции и ее производной во всех k узлах, для построения многочленов четной степени 2k используется значение функции в k узлах и значение производной в k-1 узле. Вычисления производятся по треугольной итерационной схеме Эйткена: x1 y1 V1 V12 V13 V14 V15 . V1,2m-1 y1' V2 V23 V24 V25 . x2 y2 V3 V34 V35 . y2' V4 V45 . x3 y3 V5 . y3' . . . . . . . . . . . . . . . . . . . xm ym V2m-1 ym' Элементы V1, V12, . V2m-1 первой строки этой схемы представляют собой значения многочленов Эрмита последовательно возрастающих степеней: V1 = H (1) (x), V12 = H (2) (x), . V1,2n-1 = H (2m-1) (x) . Элементы первого столбца V1, V2, . V2m-1 вычисляются по формулам: Vi = yj + yj'(x - xj) Vi+1 = yj + (yj+1 - yj)[(x - xj)/(xj+1 - xj)] , где i=1, 3, 5, . 2m-1, и j=(i+1)/2. Остальные элементы этой схемы вычисляются (вверх по диагонали) по следующей рекуррентной формуле:

где j=(i+1)/2, l=(k+3)/2. Для определения значений многчлена H (k) (x) при росте k в вычислительный процесс вовлекаются все более удаленные от x узлы интерполяции. Поэтому значения xk, yk и yk' (k=1, 2, . m) предварительно упорядочиваются таким образом, чтобы модули разностей |xk - x| возрастали при росте индекса k: |xk - x| ≤ |xk+1 - x| .

CALL PL12(X, Y, DY, XX, M, Eps, YY, Error)

X, Y, DY, XX, M, Eps - входные параметры; YY, Error - выходные параметры;

Real X(1:N), Y(1:N), DY(1:N) - массивы, содержащие значения аргумента, функции и производной. Массив X должен быть упорядочен в возрастающем порядке; Real XX - значение аргумента, в котором требуется вычислить значение функции; Integer M - число узлов, которое разрешается использовать при интерполяции, 1 SPL вычисляет коэффициенты кубического интерполяционного сплайна для функциональной зависимости yi = y(xi), i=1, 2, . n. Кубический сплайн представляет собой многочлен третьей степени на каждом из отрезков [xi, xi+1], i=1, 2, . n-1: │s1(x) = y1 + b1(x - x1) + c1(x - x1) 2 + d1(x - x1) 3 │ . . . . . . . . . . . . . . . . . . . . (1) │sn-1(x) = yn-1 + bn-1(x - xn-1) + cn-1(x - xn-1) 2 + dn-1(x - xn-1) 3

Для нахождения коэффициентов сплайна bi, ci, di используют условия сшивания для самих многочленов и для их первой и второй производной в точке xi:

и дополняем уравнением cn-1 + 3dn-1hn-1 = cn. Вводим обозначения

Имеем с помощью (2):

Подставляем в (3) с учетом (2):

Мы получили систему из n-2 уравнений и n неизвестных. Дополним полученную систему (5) уравнениями s1"(x1) = 0 и sn-1"(xn) = 0, или, что то же самое, σ1=σn=0. Тогда получаем полную систему для определения неизвестных bi, ci, di коэффициентов сплайна (1). Матрица полученной системы является симметричной и с диагональным преобладанием. Построенный сплайн минимизирует функционал b J(s) = ∫(s"(x)) 2 dx , a где a=x1, b=xn. Применение краевых условий s1"(x1)=0 и sn-1"(xn)=0 не является оптимальным, поскольку в этом случае значение сплайна вблизи граничных точек может сильно отличаться от значения интерполируемой функции. Поэтому в программе применяются следующие типы краевых условий: - краевые условия 1-го типа. s'(a) = y'(a), s'(b) = y'(b) - на концах промежутка [a,b] задаются значения 1-ой производной интерполируемой функции; - краевые условия 2-го типа. s"(a) = y"(a), s"(b) = y"(b) - на концах промежутка [a,b] задаются значения 2-ой производной интерполируемой функции; - краевые условия 3-го типа. s'(a) = s'(b), s"(a) = s"(b) - периодические краевые условия, которые применяются в случае, если интерполируемая функция является периодической с периодом T=b-a; - краевые условия 4-го типа. Применяются, когда неизвестны ни первые, ни вторые производные в граничных точках и функция не является периодической. В этом случае недостающие два уравнения получаются из аппроксимации третьих производных на концах промежутка по четырем первым и четырем последним из заданных точек. Недостающие уравнения выглядят следующим образом: для краевых условий 1 типа 2h1σ1 + h1σ2 = ∆1 - y'(x1) , hn-1σn-1 + 2hn-1σn = y'(xn) - ∆n-1 ; для краевых условий 2 типа 6σ1 = y"(x1), 6σn = y"(xn) ; для краевых условий 3 типа 2(h1 + hn-1)σ1 + h1σ2 + hn-1σn-1 = ∆1 - ∆n-1 , σ1 = σn ; для краевых условий 4 типа -h1σ1 + h1σ2 = h1 2 ∆1 (3) , hn-1σn-1 - hn-1σn = hn-1 2 ∆n-3 (3) , где ∆1 (3) = [(∆3 - ∆2)/(x4 - x2) - (∆2 - ∆1)/(x3 - x1)]/(x4 - x1), ∆n-3 (3) = [(∆n-1 - ∆n-2)/(xn - xn-2) - (∆n-2 - ∆n-3)/(xn-2 - xn-3)]/(xn - xn-3) - аппроксимации третьих производных.

CALL SPL(X, Y, ib, D1, DN, B, C, D)

X, Y, ib, D1, DN - входные параметры; B, C, D - выходные параметры.

Real X(1:N), Y(1:N) - массивы, содержащие значения аргумента и функции. Массив X должен быть упорядочен в возрастающем порядке; Integer Ib - параметр, который задает тип граничных условий, 1 SPL используется программа LA30 .

Программа SPF вычисляет значение интерполяционного ( SPL ) или сглаживающего ( SPS ) сплайна в заданной точке X. Для поиска прмежутка [xi..xi+1], в который попадает точка X, применяется двоичный поиск. Программа запоминает во внутренней переменной L номер промежутка, в котором вычислялся сплайн при предыдущем вызове. Если при последующем вызове точка X попадает в тот же промежуток, то новый промежуток можно не искать и сплайн вычисляется быстрее.

CALL SPF(X, Y, B, C, D, XX, Id)

X, Y, B, C, D, XX, Id - входные параметры; Real Spf - возвращаемое значение.

Real X(1:N), Y(1:N) - массивы, содержащие значения аргумента и функции; Real B(1:N), C(1:N), D(1:N) - массивы, содержащие значения коэффициентов сплайна; Real XX - значение аргумента, в котором требуется вычислить значение функции; Integer Id - параметр, задающий режим вычисления: id=0 - программа вычисляет значение функции; id=1 - программа вычисляет значение первой производной; id=2 - программа вычисляет значение второй производной.

Программа SPS вычисляет коэффициенты кубического сглаживающего сплайна для функциональной зависимости yi = y(xi), i=1, 2, . n. Предполагается, что точные значения yi неизвестны, а заданы с некоторой погрешностью ρi, 0 ≤ ρi ≤ 1. В случае, если все ρi=0, то коэффициенты ai совпадают с yi и сглаживающий сплайн превращается в интерполяционный. Это означает, что чем точнее заданы величины yi, тем меньше должны быть соответствующие значения ρi. Если же требуется, чтобы сплайн прошел в точности через точку (xi,yi), то соответствующую ей погрешность ρi следует положить равной нулю. Кубический сплайн представляет собой многочлен третьей степени на каждом из отрезков [xi, xi+1], i=1, 2, . n-1: si(x) = ai + bi(x - xi) + ci(x - xi) 2 + di(x - xi) 3 . Для нахождения коэффициентов сплайна ai, bi, ci, di используются условия сшивания для самих многочленов и для их первой и второй производной в точках xi:

и выполнив элементарные преобразования (см. прогр. SPL ), получаем следующую систему уравнений:

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

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

Коэффициенты сплайна находятся из уравнений:

CALL SPS(X, Y, P, ib, D1, DN, A, B, C, D)

X, Y, P, ib, D1, DN - входные параметры; A, B, C, D - выходные параметры.

Hermite

Задано некоторое множество точек (вершин) P1 . Pn на плоскости или в пространстве, n>1. Программа Hermite строит составную аппроксимирующую кривую Эрмита γ (1) ᑌ . ᑌγ (n) , которая является объединением элементарных кубических кривых γ (i) , i = 1, . n-1. Элементарная кубическая кривая Эрмита определяется двумя вершинами Pi и Pi+1 и ненулевыми касательными векторами в этих точках Qi и Qi+1 (Рис 1).

Координаты точек кривой между вершинами Pi и Pi+1 можно найти из векторного уравнения

R(t) = (1 - 3t 2 + 2t 3 ) Pi + t 2 (3 - 2t) Pi+1 + t(1 - 2t + t 2 ) Qi - t 2 (1 - t) Qi+1, 0 ≤ t ≤ 1

Это уравнение в координатной (матричной) форме может быть записано как

R(t) = GMT, 0 ≤ t ≤ 1

где (для точек P1 и P2)

Матрица M называется базисной матрицей кубической кривой Эрмита, а матрица G - её геометрической матрицей. Здесь x, y и z - пространственные координаты вершин P1 и P2, u, v, w - координаты касательных векторов Q1 и Q2.

Составная кривая Эрмита (Рис. 2) является объединением элементарных кривых Эрмита и определяется массивом вершин P1 . Pn, n>1. Каждая i-ая кривая описывается уравнением

R (i) (t) = (Pi Pi+1 Qi Qi+1)MT, 0 ≤ t ≤ 1

где M - базисная матрица, а промежуточные векторы Q2 . Qn-1 можно найти из условий равенства вторых производных в базисных точках P2 . Pn-1, которые сводятся к матричному уравнению

Свойства составной кривой Эрмита описаны в [Б5].

Программная реализация

Для начала вычислений нужно задать вершины, через которые должна пройти кривая, и векторы касательных в конечных вершинах. Программа HermiteTangents вычисляет значения касательных в промежуточных вершинах. Программа Hermite вычисляет координаты заданной точки на построенной кривой. Точка на кривой задаётся параметром t, который должен лежать в пределах от 0 до n-1. Целая часть параметра определяет номер i элементарной кривой Эрмита, а дробная часть задаёт положение точки на этой кривой.

Вызов программы HermiteTangents

call HermiteTangents(XP, xq1, xqN, XQ)

XP(1:n), xq1, xqN - входные параметры; XQ(1:n) - выходной параметр.

Real XP(1:n) - массив, содержащий x-, y- или z-координаты точек Pi, i=1, . n; Real xq1, xqN - x-, y- или z-координаты касательных векторов в точках P1 и Pn: Real XQ(1:n) - массив, содержащий x-, y- или z-координаты касательных векторов Qi в точках Pi, i=1, . n; integer n - задаёт количество точек P и размерность массивов XP и XQ;

Вызов программы Hermite

Hermite(XP, XQ, t) Программа Hermite должна быть вызвана после вызова программы HermiteTangents

XP(1:n), XQ(1:n), t - входные параметры; Real Hermite - возвращаемое значение;

Real XP(1:n) - массив, содержащий x-, y- или z-координаты точек Pi, i=1, . n; Real XQ(1:n) - массив, содержащий x-, y- или z-координаты касательных векторов Qi в точках Pi, i=1, . n (вычислены с помощью программы HermiteTangents); integer n - задаёт количество точек P и размерность массивов XP и XQ; Real t - значение параметра, в котором нужно найти x-, y- или z-координату составной кривой Эрмита, 0 ≤ t ≤ n-1

График построен программой Advanced Grapher

Программа ORT осуществляет разложение таблично-заданной зависимости yi=y(xi), i=1, . n в ряд по системе многочленов, ортогональных на множестве точек : m y(x) = ∑ αkuk(x) , k=0 где αk - коэффициенты этого разложения, а uk(x) - ортогональные многочлены последовательно возрастающих степеней: k uk(x) = ∑ βj (k) x j k=0, . m . j=0 Предполагается, что значения yi заданы с некоторой неопределённостью (с весовыми коэффициентами) pi, 0 < pi ≤ 1. Значение pi=1 означает, что соответствующее значение yi задано точно, и наоборот, чем ближе pi к нулю, тем менее учитывается данная точка при аппроксимации. Весовые коэффициенты легко вычислить, если известна погрешность ∆i величин yi, т.е. yi ± ∆i. Тогда весовые коэффициенты будут равны

где ∆max - максимальное значение погрешности. Система ортогональных многочленов строится по рекуррентному соотношению

где A=∑pixiuk 2 (xi), B=∑pixiuk(xi)uk-1(xi). A * находится из условия нормировки n ∑ piuk+1 2 (xi) = 1 . i=1 Коэффициент разложения αk определяется из формулы n αk = ∑ piyiuk(xi) . i=1 Программа последовательно увеличивает количество ортогональных многочленов, участвующих в аппроксимации заданной дискретной зависимости y(xi), пока степень старшего многочлена не достигнет заданного значения m (m≥0) или пока из-за накопления вычислительной погрешности среднеквадратичное отклонение

не перестанет уменьшаться. Достигнутая программой степень старшего многочлена возвращается в переменной Order. Для упрощения дальнейших расчетов коэффициенты всех многочленов при одинаковых степенях x суммируются и результатом выполнения программы является массив коэффициентов c0, c1, . COrder. Значение функции в точке x (x1 ≤ x ≤ xn) можно вычислить по формуле Order y(x) = ∑ ckx k . k=0 Краткое описание алгоритма построения многочленов, ортогональных на множестве точек, можно найти в [А4]. Данный метод более устойчив к накоплению ошибок вычислений, чем метод наименьших квадратов.

CALL ORT(X, Y, P, C, M, Order, Sigma)

X, Y, P, M - входные параметры; C, Order, Sigma - выходные параметры.

Real X(1:N), Y(1:N) - массивы, содержащие значения аргумента и функции. Массив X должен быть упорядочен в возрастающем порядке; Real P(1:N) - массив, содержащий значения весовых коэффициентов, 1>=P(i)>0. Значение P(i) должно быть равно 1, если значение ординаты Y(i) известно точно и стремиться к 0, если значение Y(i) можно игнорировать; Real C(0:M) - массив, содержащий коэффициенты многочлена, аппроксимирующего исходную зависимость; Integer M - наибольшая степень многочлена, которую разрешается достигнуть программе; Integer Order - наибольшая степень многочлена, которую удалось достигнуть программе, Order ≤ M; Real Sigma - среднеквадратичное отклонение, достигнутое при аппроксимации.

📎📎📎📎📎📎📎📎📎📎