CC BY
54удк: 519.652+550.3 msc2010: 65d07+86-08
О ПРИНЦИПАХ СПЛАЙН-ЭКСТРАПОЛЯЦИИ
© А. С. Костинский
Крымский федеральный университет им. В. И. Вернадского Институт сейсмологии и геодинамики лаборатория региональной сЕймичности и процессов в очагах ул. Гагарина, 20, Симферополь, 295026, Российская Федерация e-mail: kostinsky@yahoo.com
About principles of spline-extrapolation.
Kostinsky A. S.
Abstract. Possible applications of spline mathematics is discussed for situations typical for geophysical observations when only numerical values of time series of data are known, to build a physical dynamic model is either impossible or too complicated, unreasonable mainly because of complexity of geological "scene" on which the events occur. Dipmeter survey, systematic measurements of varying level and temperature of ground water, radon concentration in wells are associated traditionally with search for possible precursors of earthquakes, forecasting of the nearest following value at the sequences of such kind does not presuppose from the outset knowledge or even existence some dynamic connection between phenomena. As a consequence, an interpolation on the set of experimental points proves not approximation, in sense of convergence to something a priori existing, as well as an extrapolation it should be treated as a designing, modeling of dependence for given sum-total of the points at plane. "Principle of maximum simplicity", or Occam's Razor, may mean minimal distinction of model from polygonal line joining points, then cubic splines appear as extremals for functional of the second derivative norm. Simple idea of transferring of properties of the set of the points beyond the bounds of given interval is being found as an only way: net of knots on specified segment is supplemented by a potentially predictable point, "prognostic" spline on the augmented net is built, we must ensure minimum of integral of quadratic deviation depending on ordinate of the add-on point as a parameter. The first key moment is an expansion of "prognostic" spline in terms of system of fundamental splines, a way is opened up to simple simultaneous linear equations expressing minimum condition, analytical solution can be written in explicit form. In the second place, very essential problem is end conditions for spline interpolating experimental points, as a rule, nothing is known about the values of the first and, all the more, second derivatives at the moments of regular measurement, for example, temperature of air. One can make the choice of end values of the first derivatives not depending on interpreter if it is realized as a result of complementary optimization: the main spline must "differ to the least degree" from a cubic polynomial. For uniform net and end conditions of continuity of the third derivative of "prognostic" spline structural units of the extrapolation algorithm are represented in form of sequence of expansions in terms of coordinates
of the specified points, expansion coefficients are available analytically. In wrap-up result ordinate of the forecasted point does not depend on uniform net spacing, that is very essential applying to situations of forecasting concerning regular measurements when principal is not a value of interval between measurements but its sameness.
Keywords: interpolation, principles of extrapolation; cubic splines.
Введение: о смысле результатов интерполяции
Наша цель есть математическая схема сплайн-экстраполяции применительно к некоторой временной последовательности, предположим, что об исходных данных неизвестно ничего, кроме числовых значений. Такая ситуация обыденна в геофизике, пример — регулярные измерения химического состава, уровня и температуры грунтовых вод, концентрации радона в скважинах, наклономерные наблюдения, развернутые в районе, не исследованном с точки зрения геологии, с целью поиска предвестников ожидаемого землетрясения.
Здесь необходимы некоторые пояснения слов "неизвестно ничего". Мысль о единственной неизвестной функции, аппроксимирующей результаты измерения, естественно, рождается "тенью", сопутствующей многим и многим тысячам экспериментов доквантовой эры естественных наук. Физик, исследующий зависимость тока от приложенного напряжения для заданного проводника, получает экспериментальные точки, картина их распределения устойчива, демонстрирует те же самые особенности, если условия эксперимента меняются, например, проводник нагревается. Кажется очевидным обнаружить незначительное изменение измеренного тока при малом изменении напряжения, в сознании сама собой возникает полупрозрачная воображаемая линия графика, связь точек. В данном случае физическая картина того, что происходит, доступна для нас, анализ кинетики электронов в проводнике приводит к закону Ома, единственному варианту "ткани" связывающей точки. Точки ложатся вблизи прямой, и мы даже заранее знаем степень гладкости аппроксимирующей линии, но этот хрестоматийный пример даже для классической физики скорее исключение, с тех пор, как стало повседневностью понятие динамического хаоса.
В геофизических наблюдениях, как правило, нет опоры, подобной непререкаемости теорий физики времен Максвелла и Герца. Геофизическая картина Земли — это только "островки" связей между явлениями, "островки" как в пространстве, так и во времени, они могут появляться внезапно и исчезать навсегда и не сливаются в исследованный и освоенный "материк". Как иллюстрация, измерения крипа в точке разлома рядом с очаговой областью не показывали никаких закономерных
изменений перед землетрясением с магнитудой 5,9 в центральной Калифорнии 6 августа 1979 г., но в коровых движениях перед землетрясением Тонанкай 7 декабря 1944 г. составляющая-предвестник несомненна [1]. Только физики или химии процессов здесь недостаточно для полноты динамического описания, описания как некоторой эффективной динамической системы, приемлемо сузить набор динамических переменных не удается главным образом из-за сложности геологической "сцены", на которой происходят события. Динамическая связь означает единственное следствие для единственной причины, однозначно определенное состояние эффективной системы для заданного набора значений управляющих параметров. Как следствие, должна наблюдаться универсальная воспроизводимость результатов в некотором диапазоне меняющихся условий, но этого не происходит из-за "полифонии" интерференции причин и следствий, среди которых невозможно выделить (и, возможно, не существует) доминанты. К тому же мыслимая эффективная система наверняка не только очень сложна, но и нелинейна, существенную роль могут играть состояния перехода к хаосу. В то же время и классическое статистическое описание, и статистическая вероятность прогнозных, наиболее практически важных, оценок кажутся некорректными из-за невозможности выделить соответствующий статистический коллектив, совокупность однотипных испытаний в одинаковых условиях. Итак, мы должны с самого начала признать существование множества, ансамбля функций ф, с равной эффективностью интерполирующих последовательность исходных данных, на этом этапе в ансамбле нет функций более предпочтительных с точки зрения соответствия эксперименту, ни одна не может считаться "истинной".
Теперь, строго говоря, поиск модели, "наилучшего" варианта, есть не совсем аппроксимация. Аппроксимация в первоначальном смысле слова подразумевает некоторую сходящуюся последовательность, отражает наше стремление к определенной цели, в данном случае стремление восстановить некое "ядро" информации, скрытое в последовательности точек. Теперь перед нами скорее проблема выбора. Интерполяция-восстановление заменяется интерполяцией-конструированием, и решающую роль приобретают общие "первичные" принципы, позволяющие отсечь варианты неподходящие. В первую очередь подразумевается в определенном смысле понимаемая оптимальность, отражение всеобщности "бритвы Оккама". Нам дано только ощущение того, что в объяснении вещи не должно быть сделано больше предположений, чем необходимо, но в применении к биологическим системам это принимается как общее свойство: если существует цель наилучшего управления биологическим объектом, то он, объект, должен быть устроен максимально просто, при условии выполнения заданной функции [2].
Фактически именно "принцип максимальной простоты" [2] в чисто математической задаче конструирования "с нуля" зависимости у = f (х) по значениям ординат в узлах сетки приводит к кубическим сплайнам [3]. Действительно, пусть на отрезке [а, Ь] в узлах сетки А : а = х0 < х1 < ... < хN = Ь заданы значения Уг, г = 0,.., N, и речь идет о том, чтобы провести гладкую кривую через заданные точки (хг , Уг), то есть имеется в виду простая интерполяция без учета ошибок, интерполяция, повторимся, понимаемая как конструирование. Предположим для f (х) простейший вариант "физически разумного" графика, непрерывную дифференцируемость, то есть класс гладкости С1 на [а,Ь]. Стараясь обеспечить минимальное отличие модели от ломаной, мы находим множество экстремалей Б; х) функционала нормы второй производной
ь
з (а ) = / 1Г(ж)12 ¿х, (1)
а
и теперь, располагая произвольными коэффициентами кубических полиномов, необходимо удовлетворить условиям интерполяции
Б{М; х)|л=л. = Уг, г = 0,..., N
и достичь минимума интеграла (1). Вторая производная экстремали, реализующей минимум, должна быть непрерывна [3],
; х) е С [а, Ь],
так возникает мысль о сшивании решений и записи [3]
Б- (а; х) = Уг(1 - ¿)2(1 + 2*) + Уг+1^2(3 - 2*) + - *)2 - Шг+1^г^2(1 - *),
х е [хг,хг+1], * = (х - хг)/Нг, Нг = хг+1 - хг, г = 0,.., N - 1, (2)
= тг, г = 0,.., N, (3)
И
Б(м)(а;х)|х=х• = Уг, (а;х)
х=х Их
где последовательность тг, г = 0, ..,N определяется из системы уравнений с трех-диагональной матрицей
ЛгШг-1 + 2тг + ^¿^¿+1 = 3 ( ^г-У + ЛгУ , У'-М , { = 0, .., N - 1, (4)
\ Нг Нг-1 У
^г = Нг_ 1 (Нг-1 + Нг)-1, Лг = 1 - ^г, замкнутой, если заданы два дополнительных краевых условия. Как один из вариантов, для известных концевых значений производной
/(а) = га, а'(Ь) = 4,
х=х
минимум 3(/) приходится искать на экстремалях таких, что
= 4, (5)
-fs(N)(/; x) dx
za , is (N)(/; x)
x=b
и система (4) имеет единственное решение.
Логика, приведенная выше, оказывается простейшей из мыслимых в ситуации, когда выбор должен осуществляться в ансамбле функций, удовлетворяющих условиям интерполяции. Представим себе определенную последовательность, ряд, каждый элемент которого есть класс сложности возможных моделей ф . В иерархии этих классов д-й уровень, д > 2, есть предположения ф € С(9-1), ф(9) е Ь2[а,6], а вариационный принцип
J (q)(0)
ь 2 -q 2
-xqф(х)
dx = min, q > 2, (6)
означает минимальное интегральное отличие функции-модели ф от полинома (д— 1)-й степени. Для произвольного д уравнения Эйлера экстремалей имеют порядок 2д, как результат сшивания экстремалей степени 2д — 1 возникает сплайн 5(д'м\ Б(0'м) = 5(м\ реализующий минимум только в случае дефекта 1, Б(д'м) € С29-2. Для доказательства, аналогично книге [3], следует рассмотреть функционал (6), например, на множестве функций /(х) € Ж2[а,Ь] е [а,Ь], то есть функций, имеющих на [а, Ь] абсолютно непрерывную производную порядка д — 1, д-ю производную из Ь2[а,6] и удовлетворяющих краевым условиям
/(р)(а) = га, /(р)(Ь) = р =1,..,д — 1, (7)
а с другой стороны, сплайн степени 2д — 1 дефекта 1,
Б = Б2,-1,1, Б е С29-2[а, Ь],
построенный по тем же условиям (7), и д — 1 раз интегрировать по частям.
Запись Б(<?'М) формально может совпадать с записью общего эрмитового сплайна нечетной степени 2т +1 = 2д — 1, т = 1, 2,.., дефекта т + 1 = д и гладкости т = д — 1 [3]. Поскольку для Б(д'м) наивысший порядок непрерывной производной, 2(д — 1), есть число четное, в нашей ситуации гладкость удвоена, равна 2т, и коэффициенты такого эрмитового сплайна, то есть значения производных до порядка д — 1 в узлах сетки, не считаются заданными, а должны определяться из условий непрерывности следующих производных и дополнительных д — 1 краевых условий. Как вариант, для условий 1-го рода известны или считаются известными значения производных на концах промежутка (7).
x=a
1. Общая логика сплАйн-экстРАполяции: воспроизведение свойств последовательности "экспериментальных" точек
Пусть в есть множество, ансамбль функций ф с "физически разумным" графиком, ф е С1, интерполирующих, без учета "коридора" ошибок, последовательность данных на отрезке [а, Ь], в = {ф}, ф(хг) = Уг, г = 0, ..,N. Сплайн Б-(х) есть, несомненно, одна из точек множества в, в то же время это выделенная, "реперная" точка, обеспечивающая минимум функционала 3(ф), ф е в. Сплайн нелокален, его свойства в каждой точке зависят от последовательности в целом (хг,Уг), г = 0, ..,N, Б-(х), следовательно, есть функция-характеристика распределения узловых точек. Заметим, что при такой постановке вопроса сплайн не возникает как попытка "воспроизвести" или приблизить какую-то из функций ф, с таким же успехом любая из функций ф может рассматриваться как некоторая аппроксимация
Можно пойти еще дальше, именно — с самого начала видеть в функции Б-(х) только характеристику распределения исходной последовательности точек, словно исчезает или делается невидимым "второй план" в виде вариантов интерполяции ф. Тогда для некоторой точки х = хс, У = Ус справа от [а, Ь] кубический сплайн Б+1- (х) на расширенной сетке 8 : а = хо <х1 < ... < < +1 = хс,с дополнительным условием интерполяции в узле ж = жс,
Б^ +1)(ж)|х=хс = Ус,
есть такая же характеристика, но уже для дополненной последовательности (жг,Уг), г = 0,.., N, N + 1, Ум+1 = Ус. Задача экстраполяции — в том, чтобы в максимальной степени продолжить, перенести свойства распределения "экспериментальных" точек за пределы интервала [а,Ь], успешная в максимальной степени экстраполяция, следовательно, должна соответствовать минимальному различию характеристик Б+1-и Б(м
"Расстояние" между Би Б(м) на измеряется, строго говоря, тремя интегралами квадратичных отклонений
ь
Jifc)
dk S+1)(/; x) - S>(/; x)'2
dxk ' dxk
dx, k = 0,1, 2,
зависящими от yc как от параметра, для определенности договоримся рассматривать пока только Je(0) = 1е. Общий замысел решения задачи экстраполяции в смысле поиска минимума интеграла Ie заключается в разложении сплайна S(N+1) по системе фундаментальных сплайнов на сетке 8 (альтернатива, конечно, B-сплайны с конечным носителем [3]). Затем: пакет Maple имеет встроенную подпрограмму символического решения системы рекуррентных уравнений, ее использование применительно
к расчету сплайн-коэффициентов возможно, к сожалению, только для случая равномерной сетки. Как следствие, мы можем видеть "в буквах" решения системы с трехдиагональной матрицей для коэффициентов фундаментальных сплайнов, но с самого начала договариваясь рассматривать только ряды равноотстоящих во времени данных. Это ограничение не является принципиальным. Если в качестве приложения теории приходится иметь дело с событиями и измерениями в произвольные неравноотстоящие моменты времени, потребуется построить "рамочный" сплайн и оцифровать его с постоянным наиболее приемлемым шагом, стараясь не пропустить особенности поведения "эскизной" кривой. Итак, в качестве первого шага предстоит записать в явном виде систему линейных уравнений относительно параметров и ее аналитическое решение, краевые условия для Бдиктуют характер разложения по базису фундаментальных сплайнов и прогнозируемые параметры в точке х = хс .
2. ПОПЫТКА ПРОГНОЗА ПО ОРДИНАТЕ И ПРОИЗВОДНОЙ: НЕДООПРЕДЕЛЕННЫЕ СИСТЕМЫ
Наиболее простая и естественная схема прогноза возникает, если система уравнений для коэффициентов "прогностического" сплайна Б(м+1) замыкается краевыми условиями 1-го рода
= Ус
где у^, УС есть некоторые исходные параметры. Сплайн Б+1) представляется в виде линейной комбинации
N
S(N+1) = £ у^(х) + ycFN+i(x) + ya Фа (x) + yC Фс(х)
(8)
г=0
где ^¿(х), г = 0,.., Ж, N +1, - фундаментальные сплайны на сетке ¿, удовлетворяющие условию нормализации
N +1
^#(х) = 1, х е [а,Хе], (9)
г=0
условиям интерполяции вида
) = , г = 0,.., Ж, N +1, ] = 0,.., Ж, N +1 (10)
и нулевым краевым условиям
^(х)
= 0, —Fi(x)
dx
ж=а
0
x = tt
X=Xc
Дополнительные базисные функции ФО(ж), Фс(ж) подчинены условиям [4]
Фо(ж.) = Фс(ж.) = 0, i = 0,.., N. N +1,
ахФ(х)
1, Фа(х)
аж
ЗхФс(х)
0,
сплайновые конструкции ^¿(ж), ФО(ж), Фс(ж) определяются только расширенной сеткой
Введем обозначения варьируемых переменных
6 = Ус С2 = уО , Сэ = уС ,
и присвоим номера вспомогательным сплайнам, как указано в таблице.
Таблица 1
N ЕУгВД — Б ^ )(Х) ¿=0 +1(ж) Фо(х) Фс (ж)
номер 0 1 2 3
Подставляя разложение (8) в выражение для /е, можем записать /е в виде формы переменных , С2, СЭ :
э э э
!е = «0,0 + 2 ^ «0,к(к + ^ «т,пСт(п, (11)
к=1 т=1п=1
где коэффициент «^, 0 < ^ < 3 есть интеграл по отрезку [а,Ь] от произведения сплайнов с номерами ^. Условие минимума /е по ,С2,С3 дает систему линейных уравнений, решение записывается по правилу Крамера
а = сктт) = —^, к = 1, 2, 3,
А
Ад
(12)
Ад
«11 «12 «13 «21 «22 «23 «31 «32 «33
определители А2, получаются, как обычно, заменой столбца. Подставляя (12) в выражение (11) для /е, получаем минимальное значение интегрального отклонения в виде отношения определителей:
3
/етт)=«оо+^ «о,к &
й=1
А.
Ад;
Х = О
Х=Хс
Х = Хс
Х=О
(13)
«00 «01 «02 «03 ао1 ап «12 «13 «02 «21 «22 «23 «03 «31 «32 «33
Заметим, что для двух произвольных сплайнов вида (2) с коэффи-
циентами т(А), г = 0,..,Ж, интерполирующих нулевые значения на сетке А,
интеграл от произведения Бпо отрезку [а,Ь] выражается суммой
1 N—1
105 ^
i=0
(A,m(B, i "Ч
(A, (B,
3 4l
^ / (A, (B) (A, (B, \ 1
(A, (B,
эта запись служит основой аналитического вычисления в программе Maple коэффициентов а и затем определителей Дд и . Выясняется, однако, что не только для обычных кубических сплайнов оба этих определителя равны нулю, но попытка "раскрыть" неопределенность путем введения параметра в числителе и знаменателе (13) не меняет ситуации. Кубический нелокальный сплайн класса C1 с единственным ненулевым в узле xN параметром а, рациональный сплайн с произвольными параметрами p, q, отличными от нуля только на отрезке [xn , xc] [3], сохраняют в качестве S(N+^ Дд = 0. Возможно, сама конструкция кубического, даже обобщенного, сплайна оказывается слишком простой, "не вмещающей" информацию о производной за пределами отрезка [а,6].
3. СГЛАЖИВАНИЕ ТРЕТЬЕЙ ПРОИЗВОДНОЙ: СХЕМА ПРОГНОЗА ПО
ОРДИНАТЕ
Корректный, хотя и несколько искусственный, вариант, прогноз только по ординате, возникает, если потребовать непрерывности третьей производной S(Nв узлах x1,xn [3]. Мы имеем одну варьируемую переменную Z = Ус, сплайн S(N+^ в этом случае представляется в виде линейной комбинации
N
S(N= ^ yiFi(x) + ycFN+1(x), (14)
i=0
где Fi(x), i = 0,.., N, N +1, — фундаментальные сплайны на сетке удовлетворяющие условию нормализации (9), условиям интерполяции (10) и, опять-таки, условиям непрерывности третьей производной в узлах x1,xn .
Присвоим номера вспомогательным сплайнам, как указано в таблице.
Таблица 2
N ЕУгВД — Б ^ )(Ж) ¿=0 +1(ж)
номер 0 1
Подставляя разложение (14) в выражение для /е, записываем минимизируемую функцию переменной £,
Те = «00 + 2«01( + «п(2,
минимум по £
достигается в точке
и всегда существует, «11 = 0.
2
. / т ч «00«11 — «01
тт(7е) =-
С «11
£ (тт) = — «01 «11
Обсудим краевые условия для сплайна с номером 0. В первую очередь следует договориться о выборе краевых условий для сплайна Б^), это вопрос существенный, если преследуется цель аналитически оценить зависимость прогнозируемой ординаты от характерного масштаба задачи, доминирующего шага "основной" сетки А. Пусть предпочтение отдано естественным, наиболее простым по смыслу условиям (5), где числовые значения ,2^ получаются, например, конечно-разностной аппроксимацией или как результат дополнительного "сглаживания": Б^) должен в некотором смысле "наименее отличаться" от полинома второй или третьей степени [4]. Задача построения сплайна Б^) разрешима для любых вещественных гО,2^, так или иначе, будем пока считать ,2^ известными параметрами. Затем, без ограничения общности, можно выбрать сетку А равномерной с шагом к и записать жс — ^ = ккс, безразмерный множитель кс исчисляется в долях к.
(р. )
Можно утверждать, что коэффициенты т( , з = 0,.., N N +1 каждого из фундаментальных сплайнов •¿(ж), i = 0,.., N удовлетворяют системе N уравнений
1
13
1 + 2т(Рг) + -т<+ = — ■ (^¿,^+1 — ¿„-1) , з = 1, 2,.., N — 1, (15)
2т-1 +2т} - + ^т}+1 = ' (^¿,^+1 — ^¿,^-1), з =
для з = N
кС (Р) , о (Р) , 1 (Р)
кСГГ + ) + кт
3 ( íг,N
к(кс +1) V к,
+ М^ — ^-1) , (16)
а условия непрерывности 3-й производной в узлах х1, ХN переписываются в виде [3]
2
- = - ¿¿,0 - ¿¿,2), (17)
+ (1 - -2)ш№ + ^N+1 = - - • (+ -2 (¿^ - ¿^-1^ . (18)
Первое из уравнений (15)
2Щ™ + + 2ШТ0 = 2-(¿¿,2 - ¿¿,0)
совместно с (17) дает
mO^ + 2m(lFi) = h ( 25i,1 - 2й,о + 2^ ) , (19)
а подставляя т^у+к из (18) в (16), получаем
,/1 , и _ 1 ¿ + -с(3 + 2-с)
+1) ^ + -(-с +1)
hcmNvi,1 + (1 + h^mN^ = , ^&,N + ^ , ,C ■ (&,N - &,N—1). (20)
N
Линейная комбинация ^ у^Дх), то есть для х е [х,,х,+1], 0 < ] < N-1 функция
¿=0
N N
(1 - t)2(1 + 2t) + J] yi5i,,-+1t2(3 - 2t) +
i=0 i=0 NN
+h £ yimjFi,t(1 - t)2 - h £ yimj+1t2(1 - t) =
imf^d - t)2 - hVyimf^2'
i=0 i=0
NN
+у,(1 - *)2(1 + 2*) + у,+1 *2(3 - 2*) + +-- *)2 - - *)
¿=0 ¿=0 есть сплайн на сетке А со значениями интерполируемой функции у, и коэффициен-
N №) N
тами У] у^т, , а разность ^ у^(х) - Б^) есть сплайн на той же сетке с нулевыми
¿=0 ¿=0
значениями интерполируемой функции и коэффициентами т(0) = N УгШ,^ - ш,5 ),
¿=0
которые удовлетворяют системе уравнений
1 ш,-1 + 2ш,0) + 2ш,0)1 = 0, ^ = 1, 2,.., N - 1 и уравнениям, следующим из (19), (20),
ш00) + ЩТ)) + 2 (ш<0) + шГ,)) = - (2У1 - 2У0 + 1У2) ,
кс (тХ-1 + т^.?^ + (1 + кс) ^ + т^^ =
1 , кс(3 + 2кс) . .
УN + 7/7 , -П ■ (УN — УN-1),
ккс(кс + 1) к(кс + 1)
которые можно переписать в виде
т00) + 2т(10) + ^О = 0, кстХ-1 + (1 + кс) (т$ + ^ = 0, (21)
где обозначено
*о - 20 + 2тГ)) — к (2у, — 5У0 +1 у2) ,
_ кс (^)) . , . 1 кс(3 + 2кс) , Ч
^ = ТТк:+ + ккс(1 + кс)2 УN — к,(1 + кс)2 ■ (УN — УN-1).
4. на пути к Алгоритму, независимому от шага сетки:
ОПТИМАЛЬНЫЕ КРАЕВЫЕ УСЛОВиЯ "БАЗОВОГО" СПЛАЙНА
Обозначим через , к = 0, ..^, разрывы третьей производной сплайна Б(N) в узлах равномерной сетки А, А0, — DN принимаются равными значениям третьей производной соответственно в точках ж0, ,
6 У1 - У0
А0 = т^т т1 + т0 — 2
к2 1 0 к 6 , _ оУ^+1 — УЛ 6 , _ оУ^ — У^-1
А = + т— 2—к—) — к2 + т-1— 2 к
3 = 1,.., N — 1,
п 6 ^ oУN — УN-1
^ = — ^ т« + т^- 1 — 2
к2 V "-1 к
(для упрощения записи, пока речь идет только о сплайне Б^), не будем указывать на это в обозначениях коэффициентов т). Потребуем, чтобы сплайн Б ^) доставлял
минимум по переменным сумме
N к=0
Ё (22)
иначе говоря, Б ^) в смысле гладкости выбирается как наименее отличающийся от полинома третьей степени. Представим Б^) в виде [4]:
N
Б(^ = Ё У^(ж) + •(;) + 4 •ь(ж), (23)
¿=0
где •¿(ж), i = 0,..,N, — фундаментальные сплайны на сетке А, удовлетворяющие условиям интерполяции
•¿(ж,-) = ^, i = 0,.., N, з = 0,.., N
и краевым условиям
• (а) = ^(6) = 0, появляются базисные функции •л(ж), •Ь(ж), подчиненные условиям
•о(ж) = №) = 0, i = 0,.., N,
•0(а) = •Ь(б) = 1, •(6) = •Ь(а) = 0.
Сплайновые конструкции •¿(ж), •л(ж), •Ь(ж) определяются только сеткой А, фундаментальные сплайны теперь нормализованы так:
N
5>(ж) = 1, ж € [а, 6]. ¿=0
Из (23) следует выражение через разрывы базисных сплайнов
N
Ас (Б ^)) = Ё У¿Dfc (•¿) + 20 Вк (^о) + 4 А (•), к = 0,.., N (24)
¿=0
подставляя в (22), получаем квадратичную форму по . Условие минимума формы дает систему линейных уравнений, из которой, в свою очередь, получаются "оптимальные" краевые значения первых производных как линейные функции ординат y¿. Решающую роль играет то, что разрывы в правой части (24) могут быть просто вычислены в аналитическом виде для равномерной сетки с шагом к, и отсюда можно начать отслеживание "движения" параметра к по алгоритму экстраполяции. Именно, если записать условную общую формулу для фундаментальных сплайнов •¿(ж), i = 0,.., N,
А (•¿) = — х числитель х
h3 (-2 -V3)N - (-2 +
i = 0,1,.., N - 1,N, k = 0,.., N,
и аналогично для дополнительных базисных функций •л(ж), •Ь(ж),
(•!, •) = х числитель х -1-, к = 0,.., N
к2 (—2 — ^ — (—2 + ^
следующая таблица представляет первое резюме.
Таблица 3
индекс г индекс к множитель а числитель выражения для разрыва
0 0 6 (—4 + 3/3) (—2 — /3^ + (4 + 3/3) (—2 + /3^
0 1 6 (19 — 12/3) (—2 — /3^ — (19 + 12/3) (—2 + /3^
0 2,..^ — 1 36/3 (—2 —/3^-к + (—2 + /3^-к
0 N 36/3 1
1 0 6 (19 — 12/3) (—2 — /3^ — (19 + 12/3) (—2 + /3^
1 1 96 (—5 + 3/3) (—2 — /3^ + (5 + 3/3) (—2 + /3^
1 2 6 (289 — 168/3) (—2 — /3^ — (289 + 168/3) (—2 + /3^
1 3,.., N — 1 —144/3 (—2 —/3^-к + (—2 + /3^-к
1 N —144/3 1
Формулы для ^(ж), 2 < г < N — 2, запишем последовательно:
36/3 (—2 — у/3)^ + (—2 + /3) ^ (—2 —/3^ — (—2 + /3^
-г+М
А)(Я)
, ч 6 I Г (—2 —/3) -г+" + (—2 + /3) -г+" = — 24/3 А-У-Л-(-^-+ $
Л 1 (—2 —/3^ — (—2 + л/3)
. -г+М
. N
<*,=^ |б/51(-2 /;+(—2+/> +
Н ' ' (—2 —/3)М — (—2 + /3)
N-г+к N
(—2 —/3^-г+к + (—2 + /3)
^г-к
(—2 —/3^ — (—2 + /3)
N
к0=0
£ (—2 —/э)к-к0 —(—2 + /а)
к к0
^к0,г
+ $к+1,г — 8$к,г + 1,
k = 2,.., Ж — 1
dn(f) = 36^3 и -^3)г + (-2 + v3)
h3 (-2 -V3)N - (-2 +
и удобно представить заключительное резюме снова как две таблицы:
Таблица 4
индекс i индекс k множитель ст числитель выражения для разрыва Д (Я)
N - 1 0 -144^3 1
N - 1 1,.., N - 3 -144^3 (-2 -^3)' + (-2 + ^3)'
N - 1 N - 2 6 (289 - 168^3) (-2 - ^3^ - (289 + 168^3) (-2 + ^3^
N - 1 N - 1 96 (-5 + 3^3) (-2 - ^3^ + (5 + 3^3) (-2 +
N - 1 N 6 (19 - 12^3) (-2 - ^3^ - (19 + 12^3) (-2 + ^3^
N 0 36^3 1
N 1,.., N - 2 36^3 (-2 -^3)к + (-2 + ^3)к
N N - 1 6 (19 - 12^3) (-2 - ^3^ - (19 + 12^3) (-2 + ^3^
N N 6 (-4 + 3^3) (-2 - ^3^ + (4 + 3^3) (-2 + ^3^
Мы вычисляем коэффициенты линейного разложения
N N
4 = Ё ' = Ё ^' (25)
г=0 г=0
1 N
Рг = 1 Ё ^ №) ) - "22^,
Q
k=Q
1 N
Ф = 1 Ё ^- №)] , г = 0, .., Ж,
к=0
упрощает ситуацию свойство = — ^-г- Элементы матрицы системы, обозначенные как , зависят только от разрывов базисных функций (ж), ^ь(ж),
NN N
"11 = Ё №(^„)|2 , "22 = Ё №(^,)}2 , "12 = Ё ^ к=0 к=0 к=0
все это достаточно громоздкие выражения. Как промежуточный итог выделим зависимость от шага сетки,
^11,^22,^12 ~ П = (^) = ^11^22 — ~ , Рг, <?г ~ 1/^.
Таблица 5
индекс к множитель а числитель выражения для разрыва Д (Fa,Fb)
Fa 0 6 (—1 + /3) (—2 — /3^ + (1 + /3) (—2 + /3^
Fa 1,.., N — 1 12/3 (—2 — /3^-к + (—2 + /3^-к
Fa N 12/3 1
Fb 0 —12/3 1
Fb 1,.., N — 1 —12/3 (—2 —/3)к + (—2 + /3)к
Fb N -6 (—1 + /3) (—2 — /3^ + (1 + /3) (—2 + /3^
Из (23) и (25) следует линейное разложение коэффициента т1 .
N
)) = ^ + 4Ш1(^0) + 4Ш1(^Ь)
г=0
N
У^Уг {т^) + РгШ^К) + ^Ш^ь)} ,
г=0
и аналогично для коэффициента т)у_ 1 •
N
т5^-_1)) = ^ Уг {mN-1^) + Ргт^-1 (^о) + ^г"^-1(^,)}, г=0
это, в свою очередь, дает возможность линейного разложения по ординатам параметров ^0в соотношениях (21), замыкающих систему уравнений с трехдиагональной матрицей для первых производных сплайна с номером 0. Коэффициент а01 в разложении по £ интегрального отклонения /е есть линейная комбинация Z'a^,
«01 =
2 1 1 = —т,--Х
1260^с (5йс — 3йс/3 + /3 — 2 + 5йси2 + 3и2йс/3 — и2/3 — 2и2)2
I 9 + ^л/З Л90Ж - ножV3 - 20^hcu2 - 520ЖС - 10Жи2+ [ 1 + hc V
+ 10^u2V3 + 300^hcV3 +11 ^5 - 3^3) (1 + hc)(1 - u2)) uZ^+ + ^132 - 77^3 - 120Жи2 + 132u4 - 264u2 + 77u4V^ Zb } ,
u = (-2 + ^N ,
коэффициент aii не зависит от ординат,
= 1 132 - 77V3 - 120Жи2 + 132u4 - 264u2 + 77u4 V3
1260h2(1 + hc)2 (5hc - 3hcV3 + V3 - 2 + 5hcu2 + 3u2hcV3 - u2 V3 - 2u2)2 '
в итоге искомая координата Z(mm) точки минимума не зависит от шага сетки h.
Заметим по поводу полученного результата, что, например, регулярные метеорологические наблюдения, в сущности, не предполагают точного знания времени измерения. Не подвергается сомнению только, что промежуток времени между последовательными измерениями одинаков. В этой ситуации всегда ли допустимо предположение, что правдоподобный прогноз ближайшего следующего измерения будет зависеть от величины промежутка не явно, а только через значения ординат? Сколько-нибудь убедительного ответа на этот вопрос применительно к геофизике некоторого региона, скорее всего, не существует.
Заключительные замечания: тестирование алгоритма и
проблема CETI
Описанный выше алгоритм исчерпывает чисто математическую сторону одного из простейших вариантов экстраполяции, но оставляет открытым вопрос, что за информацию, скрытую в распределении экспериментальных точек, мы стараемся "перенести" за пределы исходного отрезка [а,6]. Сколько-нибудь полный ответ, если подразумевать конкретные обстоятельства геофизических наблюдений, — за пределами данной заметки. Тем не менее, что касается практических приложений, если исходная последовательность данных с самого начала предполагает некоторую определенную, "шифрующую" функцию, аналитическое решение оказывается имеющим некоторое отношение к проблеме CETI (Connection with Extraterrestrial Intelligence). Именно, пусть некто, назовем его Разум 1, имеет в своем распоряжении график некоторой, известной ему, функции и, измеряя на произвольно выбранном отрезке значения абсцисс и соответствующих ординат, составляет последовательный ряд пар чисел. Полученный ряд через расстояния, преодолимые только световым лучом,
передается Разуму 2, не знающему ничего о процессе конструирования данных и пребывающему в полнейшем неведении даже относительно того, существует ли хотя бы одна "физически разумная" функция, на график которой могли бы уложиться полученные им точки. Разум 1, в свою очередь, не знает ничего о месте получения сигнала и далеко не уверен, есть ли там разумные существа. Что, с его точки зрения, могло бы быть свидетельством в пользу этого? Ответный сигнал, содержащий последовательность абсцисс и ординат такую, что не просто укладывается на график выбранной им функции, но хотя бы приблизительно продолжает его после последнего посланного значения.
Итак, первое, в чем следует удостовериться, очевидно, это "первичный" аргумент в пользу качественности алгоритма экстраполяции: если для известной зависимости у = f(ж) "псевдоэкспериментальные" значения у^ = f(ж^), г = 0,..,Ж заданы на фиксированном отрезке и число N увеличивается, прогнозируемые ординаты ус должны приближаться к f (жс) для всех не слишком больших значений кс. По крайней мере для равномерной сетки это действительно так, и кривые, построенные по отклонению функции, по отклонению производной, по комбинированному отклонению
^е + ^е
шт(1е) тгп(/е1)) С С
практически неразличимы (рис.1, вверху). Минимальные по £ = ус значения отклонений к = 0,1, 2, строго говоря, должны быть функциями кс, но оказывается, что они не зависят от кс,
"Т(1е) =
_ ^з 1 (60^ - 11и2/3 + 11/3) (60^ + 11и2/3 - 11/3) 2 = 3780 77/3 - 132 + 120^2 - 77и4/3 + 264и2 - 132и4 { °} ,
тгп(Л(1)) = - к— х С 90
(-12 + 7/3) (-3 + 3и2 - 4/3^) (-3 + 3и2 + 4/3^) 2
Х (-291 + 168/3 + 3и4 - 56^2/3 + 96^2 - 168и2/3 + 288и2) { °} ,
, N
и ' " '
(-2 + /э)J
с ростом N и уменьшением шага к на фиксированном отрезке изменения минимальные значения резко падают (рис. 1, внизу).
Рис. 1. На фиксированном отрезке [0, 2п] в узлах равномерной сетки Х = hi, i = 0,..,Ж, h = 2п/Ж, заданы значения ординат y^ = sin(xj). Прогнозируется значение yc в точке xc = 2п + hhc, для Ж = 10, 20,.., 100 при меняющемся hc, < hc < h^ax, = 0.01, = п/h. Кри-
вые yc = yc(xc) с ростом Ж все точнее следуют за изгибом синуса.
Список литературы
1. Моги, К. Предсказание землетрясений / Пер. с англ. Б. А. Борисова, под ред. Л. П. Винника. — M.: Мир, 1988. — 382 с.
MOGI, K. (1988) Earthquake prediction. Moscow: Mir.
2. Романовский, Ю.М., Степанова, Н. В., Чернавский, Д. С. Математическая биофизика. — M.: Наука, 1984. — 304 с.
ROMANOVSKY, Yu.M., STEPANOVA, N.V. and CHERNAVSKY, D.S. (1984) Mathematical biophysics. Moscow: Nauka.
3. Завьялов, Ю.С., Квасов, Б. И., Мирошниченко, В. Л. Методы сплайн-функций. — M.: Наука, 1980. — 352 с.
ZAV'YALOV, Yu.S., KVASOV, B.I. and MIROSHNICHENKO, V.L. (1980) Methods of spline functions. Moscow: Nauka.
4. Вершинин, В. В., Завьялов, Ю.С., Павлов, Н.Н. Экстремальные свойства сплайнов и задача сглаживания. — Новосибирск: Наука, 1988. — 102 с.
VERSHININ, V.V., ZAV'YALOV, Yu.S. and PAVLOV, N.N. (1988) Extremal properties of splines and the problem of smoothing. Novosibirsk: Nauka.
Статья поступила в редакцию 14.11.2015
