Введение
Большинство численных методов решениязадач математического анализа так или иначе связано с аппроксимацией функций.Это и собственно задачи приближения функций (интерполяция, сглаживание,наилучшие приближения) и задачи, вкоторых аппроксимация присутствует как промежуточный этап исследования(численное дифференцирование и интегрирование, численное решениедифференциальных и интегральных уравнений).
Типичной задачей приближения являетсязадача интерполяции: по заданной таблице чисел />, восстановитьфункцию /> с той или иной точностьюна отрезке [а, b] действительной оси.Классический метод ее решения состоит в построенииинтерполяционного многочлена Лагранжа, определяемого равенством
/>
/>
Хотя согласно теореме Вейерштрассавсякая непрерывная функция /> на отрезке /> может быть какугодно хорошо приближена многочленами, практические возможности применениямногочленов Лагранжа ограничены. Прежде всего, используя подобный аппарат, мыдолжны быть уверены, что, выбрав достаточно большое число узлов интерполяции,получим хорошее приближение интерполируемой функции. Однако, как показывает рядпростых примеров, это часто нельзя гарантировать.
С. Н. Бернштейном (1916 г.) былоустановлено, что последовательность интерполяционных многочленов Лагранжапостроенных для непрерывной функции/>на отрезке [—1, 1]по равноотстоящим узлам />, с возрастанием /> не стремится к />. Еще более любопытендругой пример, восходящий к Рунге (1901 г.) и состоящий в том, что указанныйинтерполяционный процесс не сходится на [—1, 1] даже для гладкой сколь угоднораз дифференцируемой функции />(рис. 0.1). В обоихслучаях
/>
Иногда эти трудности удается преодолетьпутем специального выбора узлов интерполяции или за счет перехода к каким-либообобщенным многочленам. Однако такой путь, как правило, весьма усложняетвычисления и к тому же не избавляет нас от второй проблемы — быстрогонакопления ошибок округления с ростом степени многочлена. Поэтому на практикедля того, чтобы достаточно хорошо приблизить функцию, вместо построения интерполяционногомногочлена высокой степени используют интерполяцию кусочными многочленами.
Примером такого рода являетсякусочно-линейная интерполяция. В общем случае отрезок /> точками /> разбивается начасти и на каждом промежутке />, строится свойинтерполяционный многочлен. Полученные таким образом многочлены (обычно одной итой же степени) дают интерполяцию функции /> на всем отрезке />, которая, вообщеговоря, не обеспечивает гладкого перехода от одного звена к другому и может бытьдаже разрывной, если точки /> не включаются вчисло узлов интерполяции. Это допустимо, если не требуется восстанавливатьфункцию с заданной степенью гладкости. В частности, различные таблицысоставляются с таким шагом, чтобы промежуточные значения функции с принятойточностью можно было вычислить с помощью линейной или квадратичнойинтерполяции. Для гладкого восстановления таблично заданной функции нужноувеличить степень составляющих многочленов, а остающиеся свободнымикоэффициенты определять из условий гладкого сопряжения многочленов на соседнихпромежутках. Получающиеся при этом гладкие кусочно-многочленные функции соднородной структурой (составленные из многочленов одной и той же степени)называются сплайн-функциямиили просто сплайнами.Простейший и исторически самый старый пример сплайна — ломаная.
Термин сплайн произошел от английского spline, что в переводеозначает рейка, стержень — название приспособления, которое применяличертежники для проведения гладких кривых через заданные точки. Возьмем гибкуюстальную линейку, поставим ее на ребро и, закрепив один конец в заданной точке />, поместим междуопорами, которые располагаются так, чтобы линейка проходила через заданныеточки (рис, 0.2).
/>
Согласно закону Бернулли—Эйлералинеаризованное дифференциальное уравнение изогнутой оси линейки имеет вид
/>
где /> — втораяпроизводная прогиба, />— изгибающий момент,изменяющийся линейно от одной точки опоры к другой, /> — жесткость.Проинтегрировав это уравнение, получим, что функция />, описывающая профильлинейки, является кубическим многочленом между двумя соседними точками опоры идважды непрерывно дифференцируемой функцией на всем промежутке интегрирования.Для определенности задачи на концах должны быть заданы краевые условия, вчастности, при отсутствии внешних нагрузок па линейку />.
Функция /> представляет собойдругой пример (теперь уже гладкого) сплайна. Она относится к интерполяционнымкубическим сплайнам, обладающим рядом замечательных свойств, которые иобеспечили им успех в приложениях.
В отличие от интерполяционныхмногочленов Лагранжа, последовательность интерполяционных кубических сплайновна равномерной сетке узлов всегда сходится к интерполируемой непрерывнойфункции, причем сходимость повышается с улучшением дифференциальных свойствфункции />. Так, для функции /> из примера Рунгекубический сплайн на сетке с числом узлов /> дает погрешностьтого же порядка, что и многочлен />, но для /> она настолько мала,что в масштабах рис. 0.1 не может быть показана (ср. с многочленом />).
Алгоритмы построения кубическихсплайнов являются весьма простыми и эффективно реализуются на ЭВМ, причемвлияние ошибок округления при вычислениях оказывается незначительным.
Целью данной работыявляется изучение сплайнов, в частности базисных сплайнов. Также мы составимпрограмму для работы со сплайнами.
§1. Определение сплайнов. Пространство сплайнов
Пусть наотрезке [a,b] задано разбиение /> Для целого /> через /> обозначим множество /> раз непрерывнодифференцируемых на /> функций, а через/>— множествокусочно-непрерывных функций с точкамиразрыва первого рода.
Определение.Функция/> называется сплайном степени /> дефекта /> (/> — целое число, />)с узлами на сетке />, если
а) на каждом отрезке /> функция /> являетсямногочленом степени />, т. е.
/> (1)
б) />.
Определение сплайна имеет смыслна всей вещественной оси, еслиположить />
При этом на полуоси /> беретсятолько формула (2), а на полуоси /> только формула (1).
Итак, сплайн /> имеетнепрерывные производные до порядка />. Производныесплайна порядка выше />, вообще говоря,терпят разрывы вточках />. Для определенностибудем считать, что функция />, непрерывна справа,т. е.
/>
/>
/>
Множество сплайнов, удовлетворяющихопределению, обозначим через /> Ясно,что этому множеству принадлежат и сплайны степени n дефекта /> и сплайны степени /> дефекта />, если />, в том числемногочлены степени не выше />. Так как обычныеоперации сложения элементов из />и их умножения надействительные числа не выводят за пределы множества, то оно является линейныммножеством или линейным пространством.
Простейшим примером сплайна являетсяединичная функция Хевисайда
/>
с которой естественным образомсвязана усеченная степенная функция
/>
Функции /> являютсясплайнами соответственно нулевой степени и степени /> дефекта 1 сединственным узлом в нулевой точке (рис. 1.1). Мы будем рассматривать такжеусеченные степенные функции />,связанные с точкамисетки />.При />они принадлежатмножеству />
Теорема1.1. Функции
/>
/> (3)
/>
линейно независимы иобразуют базис в пространстве />размерности />
Доказательство:Предположим противное, т. е. что существуют постоянные />,не все равные нулю и такие, что />
Тогда для /> имеем /> и в силу линейнойнезависимости функций /> находим />Беря /> получаем /> и, по той жепричине, /> Продолжая этотпроцесс, убеждаемся, что все /> Следовательно,функции (3) линейно независимы.
Пусть теперь задан сплайн /> на отрезке /> он являетсямногочленом степени />, /> и может бытьзаписан в виде (1) или (2). При этом, так как первые /> производных сплайнанепрерывны в точках /> т.е.
/>
/>
Покажем, что сплайн />, на отрезке />может быть представленв виде
/> (4)
Где />
Действительно, преобразуя этовыражение при /> получаем
/>
/>
Это доказывает, что всякий сплайн/> может быть представленв виде линейной комбинации функций (3), т. е. эти функции образуют базис в /> и представление (4)единственно. Эта формула называетсяпредставлением сплайна в виде суммы усеченных степенных функций. Итак, множество/> являетсяконечномерным пространством размерности />
/>
§2.Базисные сплайны с конечными носителями
В математическом анализе встречаютсяконструкции, связанные с финитными функциями,т. е. гладкими функциями, которые определяются на всей действительной оси, ноотличны от нуля лишь на некотором конечном интервале (носителе). Ниже мыисследуем финитные сплайны из пространства />. В последующемизложении они играют исключительно важную роль.
Расширим сетку />, добавивдополнительно точки /> (можно положить,например, />).
Возьмем функцию/>ипостроим для нее разделенные разности /> порядкапо значениям аргумента />. В результатеполучаются функции переменной х:
/>
Так как дляразделенной разности /> порядкаот функции /> по точкам/>справедливоравенство
/>
Если использовать тождество /> томожно получить несколько иную форму записи этой функции
/>
Из определения усеченных степенныхфункций следует, что функция/>является сплайномстепени п дефекта 1 на
сетке узлов/>
Лемма 1.1. Справедливотождество
/>
Доказательство.Если/>то разделенная разностьфункции />по точкам />можетбыть вычислена по формуле Лейбница:
/>
Для разности /> порядка путемрассуждений по индукции нетрудно получить
/>
Представим функцию/>ввиде
/>/>
и построим ее разделенную разность /> порядка по формулеЛейбница. Получим
/>
Отсюда, если учестьопределение сплайнов/>, следует тождество(4).
Лемма 1.2. Сплайны /> обладают следующими свойствами:
/>
Доказательство.Функция/>равна нулю при />и являетсямногочленом степени nот х при />.Поэтому ее разделенные разности /> порядкапо значениям аргумента/>тождественно равнынулю при /> и/>т.е./> Внутриинтервала/>
В самом деле, при n = 0 согласно (2) />.Пусть, далее, утверждение а) верно при /> Тогдапри n=lв силу (4) на интервале/>функция/>являетсялинейной комбинацией с положительными весами функций/>причемпо предположению в произвольной точке указанного интервала хотя бы одна из этихфункций больше нуля. Следовательно,/>для />,и утверждение а) установлено.
Докажем утверждение б). Всякую n+1 раз непрерывнодифференцируемую функцию g(t) на промежутке а ≤ t ≤ bможно представить формулой Тейлора с остаточным членом в интегральной форме:
/>
Здесь под знаком интеграла вместо обычногосомножителя/>стоит усеченнаястепенная функция, что позволяет заменить переменный верхний предел tпостоянной величиной b. Из (7) следуетразностное соотношение
/>
то, полагая g(x) = xn+1, поручаем
/>
Поскольку внеинтервала (а, b),то это равенство -совпадает с (6) и лемма доказана.
Лемма 1.3. Функции/>являются сплайнами степени п дефекта 1 с конечныминосителями минимальной длины.
Доказательство.Предположим, что существует сплайн />отличный от нуля наинтервале, меньшем, чем />Такойинтервал, очевидно, не может иметь границей точку, не являющуюся узлом сетки />. Поэтому пусть этобудет интервал (xi, xi+n).
Возьмем представление сплайна дефекта v = 1 через усеченныестепенные функции (1.4). Вследствие того, что />при/> в этомпредставлении/>/>. Так как />при /> тоее производные до порядка n— 1 равны нулю в точке xi+n.Имеем
/>
Последние равенства представляютсобой однородную систему линейных уравнений для определения коэффициентов />/>.Ее определитель пропорционален определителю Вандермонда n-ro порядка, которыйотличен от нуля, и система имеет только нулевое решение. Наконец, из того жеусловия />следует,что />.Значит, />илемма доказана.
Теорема1.2. Функции/> линейно независимы иобразуют базис в пространстве сплайнов />
Доказательство.Покажем сначала линейную независимость функций/>навсей действительной оси. Предположим противное, т. е. что существуют такиепостоянные />,не все равные нулю, что
/>
Выбирая/>получаем, что/>и, значит, />.Беря затем /> находим, что /> ит.д., т.е. /> Следовательно,функции />линейно независимына />
Предположим теперь, что соотношение(8) выполняется только на [а, b].Это значит, что на отрезках /> обращаются в нулисплайны вида
/>
Каждый из них отличен от нуля самоебольшее на интервале />Поэтомуиз предположения />при x/>согласнодоказательству леммы 3 следует, что />0 на интервалах />, а значит, и на всейдействительной оси. В силу линейной независимости функций />на/>должнобыть /> и это для всех i = 0, ..,N-1.
Таким образом, функции/>линейнонезависимы, и так как согласно теореме 1.1 размерность пространства/>равнаn+N, то они образуютбазис в этом пространстве. Теорема доказана.
Функции /> называются базисными сплайнами с конечными носителями минимальной длины(В-сплайнами). В силу теоремы 1.2 всякий сплайн/>может быть единственным образом записан
/>
где — некоторые постоянныекоэффициенты. Эту запись сплайна называют его представлением через В-сплайны.
Из теоремы 1.2 вытекает
Следствие1.1.Всякий сплайн/>,принадлежащий />, сконечным носителем минимальной длины с точностью до постоянного множителясовпадает с В-сплайном.
Доказательство.Минимальным конечным носителем сплайна является один из интервалов /> /> Согласно (9)
/>
Так как /> для /> то, выбираяпоследовательно />, получаем, что />. Аналогично, /> для /> Следовательно, />
Замечание. Представлениесплайнов через B-сплайпы в виде (9) имеетсмысл для конечного отрезка [а, b].Чтобы получить его для всей вещественней оси, нужно положить /> и />. Тогда точки /> оказываютсяузлами кратности /> и при построении B-сплайнов с номерами /> и /> нужно учитыватьправило для разделенных разностей с кратными узлами. Мы не описываем подробноэти конструкции, ибо все практические задачи, где используются B-сплайны, рассматриваются на конечномотрезке.
§3. Нормализованные базисные сплайны ипредставление ими многочленов
При практическихвычислениях удобнее использовать не сами B-сплайны,а функции, получающиеся из них умножением на постоянные множители:
/> (1)
Эти функции называются нормализованными В-сплайнами. Нормирующий множитель равен среднему арифметическому шагов/> наотрезке, где B-сплайн отличен от нуля.
Тождество (2.4) длянормализованных 5-сплайнов имеет вид
/>
С его помощью легко можнопостроить последовательность сплайнов /> Приведем первые четыре функции этой последовательности дляслучая равноудаленных узлов hi= h.
Будем обозначать/>Точка/> — это серединаотрезка-носителя В-сплайна. Тогда имеем
/>
Эти В-сплайны изображенына рис. 1.3, а, б, в,г соответственно.
В § 1 было отмечено, чтомногочленыРп(х) степени не выше n являютсяэлементами пространства сплайнов /> .Следовательно, они представимы через базисы этихпространств, в частности через базис из В-сплайнов в пространстве />. Длявывода формул воспользуемся тождеством (2). После умножения обеих его частей начисло и суммирования по индексу iполучаем
/>
Лемма 1.4. Справедливо тождество
/>
в предположении/>
/>
Доказательство. В формуле (4) положим /> Тогда получаем
/>
Подставляя /> в (3), находим
/>
Повторяя это преобразование n раз, получим справа
/>
Теперь разложим обе части тождества(5) по степеням t.При этом
/>
Здесь /> суть символы элементарных симметрических функций от n аргументов степениа. Это многочлены, состоящие из />слагаемых. Они имеют вид
/>
Подставляяразложения (6) и (7) в (5) и приравнивая коэффициенты при одинаковых степенях t, находимпредставления мономов /> черезнормализованные В-сплайны па отрезке />
/>
Вчастности, при а = 0 получаем соотношение
/>
которое для нормализованных5-сплайпов играет ту же роль, что свойство (2.6) для самих 5-сплайнов.
Полученные формулы (8) решают вопрос опредставлении произвольного многочлена п-й степени черезнормализованные В-сплайны.
Заключение
В данной работе мырассмотрели понятие сплайна и основные определения необходимые для работы с ним.Было изучено понятие базисного сплайна или B-сплайна, а так жеуделено внимание его форме в виде нормализованного сплайна. Так же была созданапрограмма для интерполяции сплайнов при помощи языка программирования высокогоуровня C++.
Список литературы
1. Вержбицкий В.М. Основы численных методов – М.: Высш. шк.,2002.
2. Завьялов Ю.С., Квасов Б.И., МирошниченкоВ.Л. Методы сплайн – функций — М.: Наука, 1980. 352 с
3. Бахвалов Н. С., Жидков Е. П., Кобельков Г. М. Численныеметоды. Учебное пособие. — 4-е издание – М.- СПб.: Физматлит, Невский диалект,Лаборатория базовых знаний, 2003
4. Препарата Ф., Шеймос М. Вычислительнаягеометрия. Введение. — Мир, 1989
5. Колмогоров А. Н., Фомин С. В. Элементы теории функций ифункционального анализа — М.: Наука, 1976
ПриложениеЛистинг программы.
#include
#include
usingnamespace std;
doubleKubichSplain ( int n, double *X, double *Y, double Xp )
{
double*Q, *L, *A, *B, *C, *D, DXp, Yp;
int i,j, k;
Q =new double [n];
L =new double [n];
A =new double [n];
B =new double [n];
C =new double [n];
D =new double [n];
for(i=0;i
if(Xp
{
k=i;
break;
}
Q[1] =(-1) * (X[2]-X[1])/2*(X[1]-X[0] + X[2]-X[1]);
L[1] =(3*(Y[2]-Y[1])/(X[2]-X[1]) — 3*(Y[1]-Y[0])/(X[1]-X[0])) / (2*(X[1]-X[0] +X[2]-X[1]));
C[n-1]=0;
C[0]=0;
for(i=3;i
{
Q[i-1]=(-1)* (X[i]-X[i-1]) / (2*(X[i-1]-X[i-2]+X[i]-X[i-1])+(X[i-1]-X[i-2])*Q[i-2]);
L[i-1]=(3*((Y[i]-Y[i-1])/(X[i]-X[i-1])-(Y[i-1]+Y[i-2])/
(X[i-1]-X[i-2]))-(X[i-1]-X[i-2])*L[i-2])/(2*(X[i-1]-X[i-2]+X[i]-X[i-1])+(X[i-1]-X[i-2])*Q[i-2]);
}
for(i=n-1;i>=2; i--)
C[i-1]=Q[i-1]*C[i]+L[i-1];
A[0]=Y[0];
for(i=1;i
{
A[i]=Y[i];
D[i]=(C[i]-C[i-1])/3.*(X[i]-X[i-1]);
B[i]=(Y[i]-Y[i-1])/(X[i]-X[i-1])+2.*(X[i]-X[i-1])*C[i]/3.+(X[i]-X[i-1])*C[i-1]/3.;
}
DXp=Xp-X[k];
//получениезначения интерполирующей функции
Yp =A[k] + B[k]*DXp + C[k]*DXp*DXp + D[k]*DXp*DXp*DXp;
delete[]A;
delete[]B;
delete[]C;
delete[]D;
delete[]Q;
delete[]L;
returnYp;
}
voidmain ( void )
{
double*X, *Y, Xp;
int n,i;
system(«CLS»);
cout
cin>> n;
X =new double [n];
Y =new double [n];
for (i = 0; i
{
cout
cin>> X[i];
}
for (i = 0; i
{
cout
cin>> Y[i];
}
cout
cin>> Xp;
cout
getch();
delete[]X;
delete[]Y;
}