CC BY
90
ВЕСТНИК ТОМСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА
2008 Математика и механика № 2(3)
УДК 519.6
О.П. Федорова ОБ ОДНОМ ПОДХОДЕ К ПРИБЛИЖЕНИЮ ФУНКЦИИ СПЛАЙНАМИ
В настоящей работе построен сплайн, который аппроксимирует функцию одной переменной. Коэффициенты сплайна выбраны так, что интегралы по области задания функции и сплайна совпадают. Оценена погрешность сплайна. Приведены результаты численного эксперимента.
Ключевые слова: Аппроксимация сплайнами, приближение функции, оценивание погрешности приближения.
Пусть задана функция f (х), х е Б , обладающая свойством
| f (х) dx = 1.
э
Приблизим функцию сплайнами так, чтобы значения определенного интеграла по области задания приближаемой функции и сплайна совпадали. Будем использовать аппроксимацию сплайнами дефекта один.
На сетке А : -да < х_т <... < х0 < х1 <... < хп < хп+1 <... < хп+т < да , сплайн 5(х) степени т на отрезке [х0; хп ] выражается линейной комбинацией нормализованных В-сплайнов:
5 Х ( /) В1 ( Х ) , (1)
I
где (f) - последовательность линейных функционалов, вид которых определяет вид приближения. Выберем (f), так чтобы значения определенного интеграла по области задания приближаемой функции и сплайна были равными:
15 (х) с1х = | f (х) с1х .
э э
Пусть область Б совпадает с отрезком [х0; хп ], т.е. Б = [х0; хп ]. Тогда
П Х/+1 П Х1+1
15(х)dx = ^ | 5(х)dx, а | f (х)сХ = £ | f (х)Сх.
О ¡=0 XI д г=0 х
Исходя из целей построения сплайна, потребуем, чтобы
Х+1 Х+1
| 5(х)сХ = | f (х)Сх . (2)
На промежутке [х;; х;+1 ] кубический сплайн £(х) дефекта 1 (1) запишется в виде
1+2
5 (х) = Т грвр (х). (3)
р=/-1
Подставим ВД в виде (3) в формулу (2) и получим для нахождения коэффициентов систему линейных уравнений
а2г-1 + Ь121 + С2г+1 + ¿121+ 2 = I / (Х) , г = 0 П - 1 . (5)
X
Здесь
-*7+1 -*7+1 Х1+1 -*7+1
« = I Вг'-1 (ХМХ , Ъ1 = I В (Х)<& , С = I В1+1 (Х)ЛХ , ¿1 = I В+2 (х)^х.
X X X X
Для определения п+3 коэффициентов сплайна имеется п условий вида (2). Дополнительные условия можно задать различным способом, исходя из знаний о свойствах приближаемой функции Д (х). Например, потребуем, чтобы значение сплайна или его г-й производной и соответственно значения функции или ее производной совпадали в трех точках из четырех: х0, хь х„_ь х„.
£ 2рв{р)(хк) = /(г)(хк), к = 0, 1, п - 1, п. (6)
р=- 1
Система (5), (6) имеет четырехдиагональную матрицу коэффициентов, решение которой будем искать методом прогонки в виде
Щ 2М + & +2 + Я; .
Здесь ш = -с' + а, р, =----------------±------,Л, = р< -Л.
' Ь + Ш,-а, Ь, + Ш,• а, ' 6, + У,-о,
Рассмотрим функцию Л(х) = £(х) -f (х) - погрешность приближения. Если функцияf (х) е С2 \Б] и сплайн - кубический дефекта 1, то погрешность Л(х) есть дважды непрерывно дифференцируемая функция. По условию (2)
| 5 (х)с1х - | f (х)с1х = | Я (х)с1х = 0.
Или | Я (х) dx = Я (х) • ^ = 0.
х
Из этого равенства следует, что 3 х; е [х;; х;+1 ] такая, что Я (х; )=0, т.е. на каждом промежутке существует хотя бы одна точка, в которой значение сплайна £(х) и функции Дх) совпадают. Запишем для Ух е [хг-; х1+1 ] формулу Тейлора для
погрешности Л(х)= Я (х; +5;):
Я(х) = Я (X;) + Я' (е;)5г-, (7)
где 0; е [х.; х.+1 ]. Формулу (7), учитывая что Я (х;) =0, каждого промежутка можно записать в виде Л(х)= Я'(е.)5г-. Таким образом, для остаточного члена аппроксимации в области Б выполняется оценка
|Я (х)| < М 5, (8)
здесь М=тах|0;| , 5 = тах|8;|. Для построения оценки погрешности приближения
I I
можно использовать следующий подход. Построим кубический сплайн Б (х) по точкам хг-, этот сплайн будет интерполирующим как для функции Д(х),
X
так и для сплайна S(x). При выбранных граничных условиях выполняются следующие неравенства [1]:
f (x) - S (x )|| < 96 h2« (f ”). (9a)
Здесь ю(/") = maxю; (/"), ю (/") = max |/"(£)-/"(n),
i k> x; x+il
h = max hi ; i = 0, 1, ... , n — 1;
i
||S?(x) - S(x)|| <96h2«(S"). (96)
Для R(x), используя неравенства (9а) и (96), построим оценку
II/ (x) - S (x)|| <19 h 2Q ,
96
где Q = max {(/"), ю(S ")}.
Если сетка является равномерной, то на промежутке [х;; х;+1 ] отличные от нуля В-сплайны имеют вид
Bi_1(x) = (1 -t)3 /6 , в (x) = [1 + 3 • (1 -1) + 3 • t • (1 -1 )2 ]/6,
Bi+1 (x) = [1 + 3 • t + 3• t2 • (1 -t)]/6,
Bi+2 (x) = t3 /6 . (10)
Здесь t = (x - x;)/(x;+1 - x;) и, с учетом выражений для B-сплайнов (10), получим
h
| /(х)Лх = —(( + Иг; + 11гг+1 + г;+2), / = 0,п-1, (11)
Х1
где Н = х;+1 - х; - шаг равномерной сетки. Дополним систему (10) уравнениями:
1 2 1
6z- +— zr 3 0 і + T zi = 6
1 2 1
¿Z Zn-1 6 + з Zn + 7 zn+i 6
1 2 1
7z«-2 + 6 зZn-1 + 7 zn+1 6
(*я-1). (12)
Выразим из первого уравнения системы (12) г_1 = -4г0 - г1 + /(х0), подставим в уравнение (11) при і = 0 и получим следующие выражения для Ж0 = -10/7, 24 х‘г
Q0 = -1/7, Л0=— |/(х)(1х-Дх0)/7. Для получения решения на правой границе
1Н „
х0
введем обозначения а = 1/6, р = 2/3, у = 1/6 и запишем систему уравнений на правой границе в виде
а'2«-і +Р'2« +У'2«+і = /(хп), а • 2«-2 + Р • 2«-1 + У • 2«+1 = / (Хп-1 ) ,
2п-2 = ^п-22п-1 + Qn-22п + Я«-2 ■
2п-1 = Щп-12п + бп-12п+1 + Кп-1 . (13)
Первое уравнение системы (13) добавим к (5), зададим значения ап = а,
Ьп = в, сп = у, Лп = 0 и получим (п+1)-е уравнение системы (5) в виде
ап2п-1 + Ьп 2п + Сп 2п+1 + Лпг«+2 = / ( хп ) .
Вычислим коэффициенты Щ, Qi, Яі для і = 1,2,...,п . Из системы (13) выразим
2я+1 = и/У:
и = /(хп_!)-а2Яп + Жп_2Яп_х + е„_2*п) -в№_А + Ля_1) -у*я,
V = а (Шп_2шп_11¥п + шп_1дп_1 + шпдп_2) + р (жп _1шп + дп_1) +у .
Будем рассматривать сплайн как плотность некоторого распределения, тогда для первого начального момента Ц1= | х ■ Б (х) dx, интегрируя (1) для случая рав-
ноотстоящих узлов, получаем выражение
Ці= Е
і( (к2
1=0
х, к
\
120 24
(
11к2 11хк^ п-1( (11к2 11хкЛ ('-2
60
24
-Е
¡=0
40
24
к х, к
30 + "24"
Применим построенный сплайн для аппроксимации функции плотности вероятностей, например, стандартного нормального распределения
(14)
Здесь х є (-да; да). Однако, | f (х)СХ = 0,999999998026825, если х0 = -6 ; хп = 6 ,
-6
поэтому полагаем область Б = [— 6; 6].
В табл. 1 даны максимальные отклонения функции Дх) от сплайна 5(х) для различных значений числа разбиений п.
Таблица 1
Максимальное отклонение значений функции (13) от значений сплайна
для различных п
Число разбиений 10 20 30 50
Максимальное отклонение 0,010041 0,000353 0,000045 0,000002
Проведенные расчеты показали, что значения интеграла от 5(х) по области Б = [-6; 6] отклоняются от единицы менее чем на 10 - 7, при этом интегралы
X+1 - 8
| /(х)йх вычислялись с точностью порядка 10 . Значение Ць вычисляемое по
формуле (12), равно нулю также с точностью порядка 10 . Заметим, что функция
о
6
X
(14) при |х|—— да асимптотически приближается к нулю, поэтому для верификации применяемого способа задания дополнительных условий приведем результаты, построения сплайна для функции
Я(х) = 0,75х°’25 - 2,5х3 - 0,5х2+1 на промежутке Б = [0; 4]. Значения функции g(x) на границах не равны нулю:
g(0) = 1 и g(4) = 25, интеграл по области Б = [0; 4]: |g(х)dx = -13,06666(6),
э
| х • g (х) сіх = - 24. В табл. 2 приведены величины максимального отклонения сплай-
э
на от функции тах |5(х) - g(х)|, интеграла от сплайна 15(х) dx и (!1 = | х • 5(х) dx , хє° э э
вычисленные при нескольких п.
Таблица 2
Результаты расчетов для функции g(x) при различных n
Число разбиений 10 20 30 50
Максимальное отклонение 0,005066 0,001278 0,000589 0,000220
Интеграл от сплайна по промежутку [0;4] - 13,064400 - 13,066100 - 13,66415 - 13,066576
Й1 - 23,991763 - 23,997935 - 23,999082 - 23,999669
Общее решение системы разностных уравнений четвертого порядка с постоянными коэффициентами (11)
z; = с (-1) + с2 (-5 + 2^6) + с3 (-5 - 2л/б) + z*.
Здесь z* - частное решение системы (11).
Рассмотрим применение предлагаемого сплайна к приближению эмпирической плотности вероятности, построенной по модельным данным. В работе [2] предложено семейство линейных конгруэнтных методов генерирования базовой случайной величины.
В общем случае линейный конгруэнтный метод определяется рекуррентной последовательностью вида
an= an / m, где an = (вan-1 +c) mod (m), n = 1, 2,... (15)
Эта последовательность является циклической. Согласно (15), an*e{0,1,...,m}, поэтому период T < т. Максимальная длина периода Т = m достижима, об этом говорит следующая теорема [2].
Теорема 1. Максимально возможная величина периода последовательности
(15) равна m и достигается тогда и только тогда, когда:
1) с и m взаимно просты,
2) в-1 кратно р для всех простых р, являющихся делителем m,
3) в-1 кратно четырем, если m кратно четырем.
Для проведения численного эксперимента будем использовать генераторы базовой случайной величины с параметрами с = 51, х0 = 5, в = 101 и с = 7, х0 = 2, в = 81 и длиной периода m = 1 000 000. Методом полярных координат построим последовательность случайных чисел, имеющих стандартное нормальное распределение. Принадлежность используемых для численного эксперимента выборок к нормальному распределению проверялась при помощи критерия согласия х2.
Анализ результатов численного эксперимента показал, что сплайн, приближающий эмпирическое распределение, демонстрирует сходимость по норме пространства С [-6; 6] к плотности вероятности стандартного нормального распределения (рис. 1).
Рис. 1. Сплайн-аппроксимация S(x) эмпирической плотности вероятности и график плотности стандартного нормального распределения
В табл. 3 показаны значения тах |£(х) - /(х)|, вычисляемые при различных п.
хе[-6;6]
Таблица 3
Объем выборки п = 10 п = 30 п = 100
5000 0,008653 0,033278 0,078358
50000 0,004839 0,009868 0,022497
500000 0,000846 0,003912 0,008984
1000000 0,000862 0,003321 0,004401
Здесь f (х) - функция плотности вероятностей стандартного нормального распределения. По построению, 5(х) для каждого разбиения есть плотность некоторого распределения. Задавая увеличивающееся число разбиений к, получаем последовательность функций 5к(х), которая отклоняется от эмпирической плотности на величину порядка 0(к2), т.е. при увеличении объема выборки и увеличении числа разбиений так, что к ^ 0, сплайн-функция сходится к плотности стандартного нормального распределения со вторым порядком малости относительно к .
ЛИТЕРАТУРА
1. Завьялов Ю.С., Квасов Б.И., Мирошниченко В.Л. Методы сплайн-функций. М.: Наука, 1980. 352 с.
2. Кнут Д. Искусство программирования для ЭВМ: В 3 т. М.: Мир, 1977. Т. 2: Получис-ленные алгоритмы. 728 с.
Статья принята в печать 05.07.2008 г.
