Учебник - Практические занятия по вычислительной математике - Аристова (1238839), страница 23

Просмтор этого файла доступен только зарегистрированным пользователям. Но у нас супер быстрая регистрация: достаточно только электронной почты!

Текст из файла (страница 23)

задачи на доказательство), для которой алгебраический интерполяционный процесс наравномерной сетке не сходится, являются причинами того, что дляфункций, заданных таблично на большом числе точек, не используютинтерполянты высоких степеней. Обычно в таких случаях переходят кпостроению сплайнов. Рассмотрим вначале более простую задачу кусочно-многочленной интерполяции.Пусть функция f (x) задана таблицей. Для восстановления функциимежду узлами можно воспользоваться функцией, которая между кажды150VI.4.1. Оценка неустранимой погрешности при интерполяции. Выборстепени кусочно-многочленной интерполяциими двумя соседними узлами является многочленом заданной невысокойстепени, например, первой, второй, третьей и т. д.Соответствующая интерполяция называется кусочно-линейной,кусочно-квадратичной и т.

п. Такую интерполяцию можно строитьлокально, не привлекая данные с отдаленных участков задания функции,а можно строить интерполяцию, пытаясь сохранить максимальную гладкость построенной функции. В последнем случае все коэффициенты локальной интерполяции (кроме кусочно-линейной) оказываются связанными между собой.VI.4.1.

Оценка неустранимой погрешности при интерполяции.Выбор степени кусочно-многочленной интерполяцииПусть функция f (x) определена на отрезке [0, π] и пусть заданы еезначения в узлах равномерной сетки xk = kπ/n, k = 0, 1, …, n. По таблицеf (x0), f (x1), …, f (xn) в принципе нельзя восстановить функцию f (x) точно, потому что значения различных функций могут совпадать в точкахxk, k = 1, …, n, т. е.

различные функции могут иметь одинаковую таблицу.Если о функции известно лишь то, что она непрерывна, то ее нельзявосстановить в точке x ≠ xk, k = 0, 1, …, n, ни с какой гарантированнойточностью.Пусть о функции f (x) известно, что она имеет производные порядкаs + 1, причемmax | f ( s 1) ( x) |  M s  const.(4.1)xУкажем две функцииf ( I) ( x ) sin nxn s 1,f (II) ( x)  sin nxn s 1(4.2),для которых таблицы f(I)(xk) = f(II)(xk), k = 0, 1, …, n, совпадают (обе таблицы содержат лишь нули).

Эти функции уклоняются друг от друга навеличину порядка hs + 1:max | f (I)( x)  f (II)( x) |  2 maxxsin nxxns 1 2h s 1 .(4.3)Таким образом, зная лишь оценку s + 1 производной, в принципенельзя восстановить табличную функцию с точностью, большей, чемO(hs + 1). Данная погрешность неустранима.151VI. ТАБЛИЧНОЕ ЗАДАНИЕИ ИНТЕРПОЛИРОВАНИЕ ФУНКЦИЙVI.4.2.

Насыщаемость (гладкостью)кусочно-многочленной интерполяцииПусть функция f(x) определена на отрезке [a, b], и задана ее таблицаf(xk) в равноотстоящих узлах xk, k = 0, 1, …, n; с шагом h = (b – a)/n.Погрешность кусочно-многочленной интерполяции степени s (с помощью интерполяционных многочленов Ps(x, fkj) на отрезке xk ≤ x ≤ xk+1)в случае, если на [a, b] существует и ограничена f(s+1)(x), имеет порядокO(hs+1).Если о функции f (x) известно лишь, что она имеет ограниченнуюпроизводную до некоторого порядка q, q < s, то неустранимая погрешность при ее восстановлении по таблице есть O(hq+1). Можно показать,что при интерполяции с помощью Ps(x, fkj) порядок O(hq+1) достигается.Если f (x) имеет ограниченную производную порядка q + 1, q > s, топогрешность интерполяции с помощью Ps(x, fkj) остается O(hs+1), т. е.порядок погрешности не реагирует на дополнительную, сверх s + 1 производной, гладкость функции f (x).

Это свойство кусочно-многочленнойинтерполяции называют свойством насыщаемости (гладкостью).VI.4.3. Нелокальная гладкая кусочно-многочленная интерполяцияПусть задана таблица функции на некотором отрезке. Поставим задачунайти на каждом отрезке разбиения xk ≤ x ≤xk+1, k = 0, 1,…, n – 1, кубическиймногочлен P3(x, k) так, чтобы возникающая при этом на отрезке a ≤ x ≤ b кусочно-многочленная функция совпадала с заданной сеточной функцией в узлах иимела непрерывные производные до порядка s = 2. Разность между степеньюинтерполяционного полинома и глобальной гладкостью сплайна называетсядефектом сплайна.

Точки, в которых заданы значения функции, называютсяузлами интерполяции, а точки, в которых сшиваются интерполяционные многочлены, – узлами сплайна. Вообще говоря, узлы интерполяции и узлы сплайна могут не совпадать.Будем строить кубический сплайн дефекта 1 при условии совпадения узлов интерполяции и узлов сплайна. Кубический многочлен задается четырьмякоэффициентами на каждом интервале, всего 4n коэффициентов. Условие интерполяции на каждом из отрезков разбиения слева и справа даст 2n условий,требование непрерывности первой и второй производной во внутренних узлахприведет еще к 2(n – 1) условию, всего 4n – 2 условия для нахождения 4n коэффициентов.Недостающие условия можно задавать различными способами. Наиболееупотребляемыми являются следующие:152VI.4.3.

Нелокальная гладкая кусочно-многочленная интерполяция1) «свободный» сплайн, соответствующий минимуму потенциальнойэнергии упругой линейки, поставленной на ребро и закрепленной так, чтобыона проходила через заданные точки:d 2 P3 ( x, 0)dx2d 2 P3 ( x, n)dx20;(4.4)2) кубический сплайн Шонбергаd 3P3 ( x, 0)dx3d 3c0,dx3d 3P3 ( x, n)dx3d 3cn.dx3(4.5)Здесь c0(x), cn(x) – единственные кубические кривые, которые проходятсоответственно через четыре первые и четыре последние из заданныхточек;3) эмпирический сплайн, требующий непрерывности третьей производной в первой и предпоследней точке интервала;4) естественный сплайн, обеспечивающий минимизацию разрывовпоследней существующей производной.Теорема 12.

Интерполяционный кубический сплайн S(x) дефекта1, удовлетворяющий одному из краевых условий (4.4–4.5), существует иединственен.Рассмотрим построение кубического сплайна в общем случае на неравномерной сетке xn – xn - 1 = hn – 1, xn + 1 – xn = hn. Пусть mn – значениевторой производной в точке xn (пока неизвестное!). На каждом отрезке[xn, xn + 1] для второй производной кусочно-кубического сплайна имеем S xx1hn mn( xn1  x)  mn1 ( x xn )  .Так как сплайн – полином третьей степени, то его вторая производная – линейная функция. Интегрируем выражение для второй производной сплайна, получаем на отрезке [ xn , xn 1 ]S x 1hn( xn1  x)2( x  xn )2 mn mn 122  An .An – константа интегрирования.

Интегрируя последнее соотношение ещераз, получаем153VI. ТАБЛИЧНОЕ ЗАДАНИЕИ ИНТЕРПОЛИРОВАНИЕ ФУНКЦИЙS ( x) 16hn mn ( xn1  x)3  mn1( x  xn )3   An x  Bn .ПослевторогоинтегрированияположимAn + Bn = (βn - αn)x + αnxn+1 - βnxn, вместо двух констант интегрированияAn, Bn введем две новые константы, более удобные для дальнейших выкладок. Тогда последнее равенство перепишем в видеS ( x) 16hn mn ( xn1  x)3  mn1( x  xn )3  n ( xn 1  x)  n ( x  xn ).Из условий интерполяции S(xn) = fn, S(xn + 1) = fn + 1 получаемfm hm h2f n  n n   n hn   n  n  n n .6hn6f n 1 fm hmn 1hn2  n hn  n  n 1  n 1 n .hn66ff(m mn )hnAn  n1 n  n1.hn6Приравняем первые производные в каждом узле сетки, кроме граничных, справа и слева Sx'(xn + 0) = Sx'(xn - 0), получим систему уравнений для определения коэффициентов сплайна:mn hn 1 mn 1hn 1 f n  f n 1 (mn  mn 1 )hn 122hn 16mn 1hn mn hn f n 1  f n (mn 1  mn )hn.22hn6(4.6)Эта система дополняется соответствующими граничными условиями.

В случае свободного сплайна m0  mN  0.Для определения коэффициентов mn получена система линейныхуравнений с трехдиагональной матрицей. Матрица этой системы симметрична, имеет свойство диагонального преобладания и, как можнопоказать, положительно определена, а следовательно, неособенная. Значит, решение рассматриваемой СЛАУ существует и единственно. Следовательно, и задача о построении кубического сплайна имеет един154VI.5. Дробно-полиномиальные аппроксимацииственное решение. Для других типов краевых условий доказательствопроводится аналогично.Коэффициенты mi называются моментами кубического сплайна.Теорема 13.

Пусть S3(x) – кубический сплайн дефекта 1, интерполирующий на системе узлов a = < x0 < x1 < …< xn = b четыреждынепрерывно дифференцируемую на [a, b] функцию f(x). Тогда дляn c  0 : x  [a, b] справедливо неравенствоf ( x)  S3 ( x)  c4 , где  max hi .i 0,...,n 1VI.5. Дробно-полиномиальные аппроксимацииVI.5.1. Рациональная аппроксимацияВ ряде случаев большую точность приближения можно получить,используя рациональную интерполяцию. Это важно в тех случаях, когдаинтервал аппроксимации бесконечен и поведение функции на бесконечности не описывается степенной зависимостью или при наличии у функции полюсов.При заданных значениях функции в узлах f(x1), …, f(xn) приближение к f(x) ищется в видеR( x) a0  a1 x  ...

 a p x pb0  b1 x  ...  bq x q, n  p  q  1.(5.1)Коэффициенты ai, bi определены с точностью до общего множителяи находятся из условий R(xj) = f(xj), j = 1,...,n, илиpqk 0k 0 ak xkj  f ( x j )  bk xkj  0,j  1,..., n.(5.2)Система (5.2) представляет собой систему n уравнений относительноn + 1 неизвестного. Превышение количества неизвестных над количеством уравнений является естественным, т.к.

числитель и знаменательопределены с точностью до общего множителя.Существует рекурсивный алгоритм вычисления функции R(x) в случаях, когда n – нечетно и p = q и когда n – четно и p = q + 1.Для начала определим два набора величин:n величин R1,0 = 0 R2,1 = 0…Rn,n – 1 = 0,155VI. ТАБЛИЧНОЕ ЗАДАНИЕИ ИНТЕРПОЛИРОВАНИЕ ФУНКЦИЙn величин R1,1 = f1 R2,2 = f2…Rn,n = fn.После этого с помощью формулыRi,k  Ri 1,k Ri 1,k  Ri ,k 1 x  xi x  xkRi 1,k  Ri ,k 1   1 1   Ri 1,k  Ri 1,k 1 (5.3)последовательно вычисляютсяn величинR1,2, R2,3 …n – 1 величин R1,2, R2,3Rn – 1,n,…Rn – 1,nи так далее до вычислениядвух величинодной величиныR1,n–1и R2,n,R1,n,которая и является ответом.Все деления в (5.3) должны сопровождаться проверками на нуль взнаменателе, например, при x = xk вся большая дробь в (5.3) должна бытьположена равной нулю.Как и для полиномиальной интерполяции, алгоритм позволяет оценить погрешность по последней поправке – разнице междуR1,n и R1,n – 1 или R2,n .Пример.

Рассмотрим построение рациональной интерполяции дляфункции f(x) = arctg x на интервале 0  x < ∞. Существует пределf (+∞) = π / 2. Это означает, что в формуле (5.10) нужно выбрать p = q.Кроме того, при x  0 функция f(x) ≈ x и раскладывается в ряд по нечетным степеням x. Для того чтобы удовлетворить обоим требованиям прине слишком большой степени полиномов, удобнее аппроксимировать неf (x), а ее квадрат:f 2(x) = arctg2 x ≈ Pp(x2) / Qp(x2).Для удовлетворения граничным условиям надо полагатьa0 = 0, a1/b0 = 1,ap = (π/2)2, bp = 1.В этом случае свободными остаются 2p –2 параметров. Столько же нужно брать узлов интерполяции на [0, +∞].

Характеристики

Список файлов книги