CC BY
29
УДК 004.021 004.043
Обобщение методов аппроксимации наборов дискретных
данных
И. А. Маркова
Кафедра прикладной математики и информатики Международный университет природы, общества и человека «Дубна» ул. Университетская, д. 19, Дубна, Московская область, Россия, 141980
В процессе математического моделирования нередко возникает необходимость в гладком приближении различных зависимостей, заданных дискретно или графически, особенно, если значения таких зависимостей приходится получать в результате проведения сложных экспериментов или из громоздких расчётов. Обратное преобразование непрерывных моделируемых объектов в дискретный цифровой формат, который используется для их хранения и компьютерной обработки, тоже требует определённого упорядочения.
Предварительно принимается, что гладкое приближение заданных дискретных точек на плоскости осуществляется линейной аналитической моделью. Условия интерполяции приводят к системе линейных уравнений с квадратной матрицей.
При интерполировании полиномами путём разбивки и перегруппировки членов степенного ряда можно получить такие базисные функции, как полиномы Лагранжа, или полиномы Бернштейна. Другими способами интерполяции являются полиномы Ньютона, итерационный процесс Эйткена и т.д.
Однако перечисленные методы реализуют лишь некоторые частные случаи из всех возможных вариантов приближения дискретных данных произвольными базисными функциями и в основном ориентированы на вычисления вручную. В компьютерных вычислениях желательно найти по возможности общий алгоритм решения, чтобы избежать программирования многочисленных частных случаев.
Рассмотрена задача обобщения существующих методов аппроксимации наборов дискретных данных (обобщённый алгоритм) и приведения этих дискретных данных к единому виду (дискретная унифицированная структура).
Ключевые слова: дискретная структура, линейная аналитическая модель, условия интерполяции, аппроксимация, метод коллокаций, метод Галёркина.
1. Введение
Предварительно принимается, что гладкое приближение заданных дискретных точек {xi,yi} на плоскости осуществляется линейной аналитической моделью (например, степенным рядом):
/ (ж, Со, ...,Сп) = Со^о(х) + С1^1(х) + ... + Сп(рп(х),
(1)
где Сз — некоторые коэффициенты, (х) — базисные функции (линейно независимые).
Условия интерполяции, то есть условия точного выполнения равенств /(хг) = Уг в точках XI, (г = 0,1,... ,т), приводят к системе линейных уравнений с квадратной матрицей плана Р (причём п = т)
( <Ро(хо) <Р!(хо) <ро(хг) (¿>1(хг)
т ) <Р1(хт)
<Рп(х о) \
) /
( Со \
С1
\сп)
( Уо \
У1
\ Уш )
или Р ■ С = У. (2)
Статья поступила в редакцию 30 декабря 2013 г.
При аппроксимации минимизируется сумма квадратов отклонений значений f (Хк) от значений ук на заданной системе точек Хк(к = 0,1,... ,т) (метод наименьших квадратов). Необходимое для этого условие экстремума приводит к системе нормальных уравнений с квадратной матрицей N размера (п + 1) х (п + 1) (причём должно соблюдаться условие п ^ т):
N -С = В.
(3)
Можно показать, что N = РТ ■ Р; В = Рт ■ У или в развёрнутом матричном виде:
( Ро(хо) ф0 (жх) ^х(жо) (хх)
<Ро(хт) \
\ <£п(хо) ^п(хх) ... <£п(хт) )
( ^о(хо) ^1(жо) ^о(Ж1) ^1(^1)
V <£о(Х ) <Р1(хт)
Рп(хо) \ / Со \ <рп(х 1) _ С1
'■Рп(хт) / \ Ст )
( Уо \
У1
\ Уш )
В вычислительной практике прямого решения систем линейных уравнений (2) и (3) стараются избегать из-за того, что они часто получаются плохо обусловленными [1] и существуют более простые способы, основанные на специально сконструированных базисных функциях [2,3]. Например, при интерполировании полиномами путём разбивки и перегруппировки членов степенного ряда можно получить такие базисные функции, как полиномы Лагранжа, или полиномы Бернштейна, в которых коэффициенты С^ принимают значения, равные заданным значениям у^. Другими способами интерполяции являются полиномы Ньютона, итерационный процесс Эйткена и т.д.
Однако перечисленные методы реализуют лишь некоторые частные случаи из всех возможных вариантов приближения дискретных данных произвольными базисными функциями и в основном ориентированы на вычисления вручную. В компьютерных вычислениях желательно найти по возможности общий алгоритм решения, чтобы избежать программирования многочисленных частных случаев.
х
2. Обобщённый алгоритм
Предлагается алгоритм, в котором вычисляемое значение f в точке х : / = /(х,Ск), определяется из условия разрешимости расширенной однородной системы уравнений. Эта система получается из системы (2) или (3) переносом из её правой в левую часть столбца У или В и добавлением к ней строки: Со^о(х) + С1^1(ж) + ... + Сп^п(х) - / = 0
/
Я
\ Фо(х) ^1(ж)
<Рп(х)
/о \
/1
/п
( Со \ С1
Сп 1
0
0
0 0
или Н ■ С = 0, (4)
где через Ц обозначена матрица Р или N, через ^ значения -у^ или —bi.
Линейные преобразования Гаусса приведут квадратную матрицу Н размера г х г ((п + 2) х (п + 2)) к треугольному виду, и (Н) будет равен произведению
её диагональных элементов. Если в этом процессе не подвергать перестановкам нижнюю строку, то условие разрешимости однородной системы (4) примет вид:
(Н) = ёе! ■ (-/ + Нгг) = 0, (5)
где Нгг — линейная комбинация значений /¿, добавленная в нижний правый элемент матрицы Н после преобразований Гаусса.
Поскольку здесь ёе! = 0, что необходимо для существования единственного решения системы (2) или (3), то должно выполняться соотношение (—/ + Нгг) = 0 или Нгг = f. Таким образом, алгоритм вычисления f состоит в обнулении правого нижнего элемента матрицы Н и последующего прямого хода алгоритма Гаусса без перестановок нижней строки. Накопленная в правом нижнем элементе матрицы Н линейная комбинация значений fi даст искомое значение f.
Следует добавить, что полученный алгоритм остаётся справедливым не только при вычислении в точке х значения f, но и при вычислении результата любой операции над рядом (1), сохраняющей его линейность относительно коэффициентов С к. В частности, аналогичным путём будет вычисляться производная f' в точке х или интеграл Р на интервале [х-\,х\ после замены нижней строки системы (4):
Со^о(х) + С\ф1 (х) + ... + Спрпх) — / = 0 на продифференцированную (6), или проинтегрированную (7) строку
С0^'0(х) + СМх) + ... + Сп^'п(х) — /' = 0, (6)
X X
С-1 + Со у уо(х)ёх + ... + Сп ! (х)ёх — Р = 0. (7)
3. Дискретная структура
В случаях, когда серия вычислений производится при неизменных исходных дискретных данных (точки {хг,Уг}), можно не повторять алгоритм Гаусса для всей матрицы Н. Поскольку в этих случаях будет меняться только нижняя строка, то приведённые к треугольному виду строки матрицы Н, кроме нижней, следует сохранять в памяти и применять алгоритм Гаусса только к нижней строке
( Лц 0
0
¿12 ¿22
0
\ фо(х) ф1(х)
¿1,п+1 ¿1,п+2 \ ¿2,п+1 ¿2,п+2
Лп+1,п+1 Лп+1,п+2
фп(х)
0
(8)
7
Здесь через обозначены элементы дискретной структуры. Полученная треугольная дискретная структура размера
(п + 2)(п + 3)
ис-
пользуется в данном случае вместо набора коэффициентов Cj размера п +1 из традиционного подхода. Такая реструктуризация сохраняемых данных, увеличивающая с одной стороны их объём, с другой стороны окупается тем, что предоставляет следующие возможности:
— получать результаты, в том числе производные и интегралы, путём применения единого алгоритма к одинаковым структурам данных, практически исходя из любого набора базисных функций;
— уменьшить потери точности в вычислениях, поскольку выполняется только прямой ход алгоритма Гаусса, а обратный ход исключается.
Известно, что неточность в задании элементов матрицы в случае плохо обусловленной системы уравнений приводит к значительным погрешностям при вычислении решения (вектора коэффициентов С). В предложенном алгоритме процесс решения выполняется частично (только прямой ход алгоритма Гаусса), что, однако, не отменяет чувствительности результата к погрешностям в исходных данных, которые определяется обусловленностью матрицы Н.
Помимо исходных данных {xi, у^}, на обусловленность системы (8) можно влиять путём подбора базисных функций (pj (х), также переносом начала системы координат в другую точку {хо,уо} и масштабированием координатных осей, умножая их на коэффициенты зх, 8У. Итоговые линейные преобразования зх • (х — Хо) и 8У • (у — У о) не изменят линейности исходной системы уравнений (4).
Для оценки изменённой обусловленности системы уравнений можно применить критерий, изложенный в работе [4].
Универсальность данного алгоритма позволяет распространять область его применения на новые задачи, что приводит к следующим возможным обобщениям в процессе формирования расширенной матрицы Н: 1. В соответствии с соотношениями (6) и (7) можно включать в матрицу плана Р строки с желаемыми или измеренными значениями наклонов касательных у' в точках Хг, т.е. строки:
Со^о(хн) + С1<р[(хг) + ... + Сп<р'п(хг) — у\ = 0 (задача Эрмита),
и/или строки со значениями интегральных площадей Уг на некоторых подынтервалах [х-1,хг], т. е. строки:
В принципе, если возникнет такая практическая необходимость, ничего не мешает повторять дифференцирование и интегрирование и включать в нижнюю строку системы (4) и/или в матрицу плана Р, соотношения для производных различных порядков (р = 1,2,3,...) и/или для повторных интегралов различной кратности (р = —1, —2, —3,..., добавляя слева новые коэффициенты Ср).
Алгоритм далее можно обобщить на двумерные зависимости, для вычисления в некоторой точке {х, у} их значений и частных производных, а по плоской подобласти Б - двойных интегралов.
Двумерную аналитическую линейную модель следует представить в виде, аналогичном ряду (1):
где С к — неизвестные коэффициенты ^ (х,у) = иг(х) • (у), где иг(х) и Vj (у) — одномерные базисные функции.
Дополнительно к дискретным значениям в точках {хг, yj} в матрицу плана Р можно тоже включать строки со значениями наклонов касательной плоскости по направлению х и/или по направлению у, а также для объёмов по некоторым подобластям в плоскости Оху.
Можно продолжить распространять данный подход на многомерные зависимости с целью вычисления (нахождения) в том числе градиента (от скалярной функции нескольких переменных), матрицы Якоби (от вектор-функции
4. Последующие обобщения
/ (x, у) = Cо>Mx, у) + у) + ... + Сп^р^х^ y),
(9)
нескольких переменных), многомерного интеграла по некоторой подобласти независимых переменных.
2. Обширной неосвоенной областью приложений предложенного алгоритма могут служить задачи построения приближённых решений уравнений (функциональных, дифференциальных, интегральных), т. е. уравнений, решениями которых являются функции одной переменной. Существующие численные методы построения приближённых решений перечисленных уравнений основаны, в основном, на их интерполяции полиномами (частный случай линейной модели (1)).
В данном случае требуется решать обратную задачу - исходя из уравнений F (t, y(t), y'(t)J y(t)dt,...) =0, связывающих зависимые (у) и независимые (t) переменные, подобрать удовлетворяющие уравнению элементы дискретной структуры (8). Поскольку зависимая функция является неизвестной, то критериями близости вычисляемой функции f (t) к точному решению y(t) могут служить используемые в вычислительной практике:
— метод коллокаций, состоящий в точном выполнении равенства ну-
— метод наименьших квадратов, вытекающий из условия минимума суммы квадратов всех невязок на заданной системе точек ¿¿;
— метод Галёркина, основанный на требовании ортогональности базисных функций (!) и невязки:
После выбора критерия можно приступать к разработке алгоритма, позволяющего построить последовательное приближение элементов дискретной структуры (8) с нужной точностью.
В перспективе ничто не препятствует переходу к уравнениям, решениями которых являются функции двух и более переменных.
1. Райс Д. Матричные вычисления и математическое обеспечение. — М.: Мир, 1984. [Rice JohnR. Matrix Computation and Mathematical Software. — Moscow: Mir, 1984. — (in russian). ]
2. Бахвалов Н. С., Жидков Н. П., Кобельков Г. М. Численные методы. — М.: БИНОМ. Лаборатория знаний, 2006. [Bakhvalov N. S, Zhidkov N. P, Kobelkov G. M. Numerical Methods. — Moscow: BINOM, 2006. — (in russian). ]
3. Форсайт Д., Малькольм М., Моулер К. Машинные методы математических вычислений. — М.: Мир, 1980. [Forsythe G.E., Malcolm M.A., Moler C.B. Computer Methods for Mathematical Computations. — Moscow: Mir, 1980. — (in russian). ]
4. Митюков В. В. Наглядная геометрическая оценка обусловленности линейных преобразований // «Решетневские чтения»: Материалы Х-й Международной научной конференции (8-10 ноября 2006 г.). — 2006. — С. 253254. [Mitjukov V. V. Evident Geometrical Estimation of Conditionality of Linear Transformations // "Reshetnev Readings": Materials of 10th International scientific conference (November 8-10, 2006). — 2006. — Pp. 253-254. — (in russian). ]
для j = 1, 2,... ,n.
0
Литература
UDC 004.021 004.043
The Generalization Known Methods to Approximate Various
Sets of Discreet Data
I. A. Markova
Department of Applied Mathematics and Informatics Dubna International University for Nature, Society and Man 19, University str., Dubna, Moscow region, Russia, 141980
In the process of mathematical modeling a necessity often arises to smoothly approximate various dependencies which are defined discretely or graphically. Especially if the value of such dependencies are obtained as a result of complex experiments or cumbersome calculations. The inverse transform of continuous simulated objects in discrete digital format that is used for storage and computer processing also requires a certain ordering.
It is assumed beforehand that the smooth approximation of discrete set of points on the plane is performed by linear analytical model. Interpolation conditions lead to a system of linear equations with a square matrix.
When interpolating polynomials by breaking and rearranging the terms of a power series one can get such basic functions as Lagrange polynomials or Bernstein polynomials. Other methods of interpolation are Newton polynomials, Aitken iterative process, etc.
However, these methods realize only some particular cases of all possible approximations of discrete data by arbitrary basis functions and are mainly focused on manual calculations. In computer calculations, it is desirable to find a general algorithm for solutions in order to avoid programming many particular cases.
The problem of generalization of existing methods for approximation of discrete data sets (generalized algorithm) and bringing these discrete data to a common form (a discrete unified structure) is considered.
Key words and phrases: discrete structure, the linear analytical model interpolation conditions, approximation, collocation method, Galerkin method.
