Методичка (Г.А. Кокотушкин, А.А. Федотов, П.В. Храпов - Численные методы алгебры и приближения функций), страница 5

Интерполяционный многочлен ЛагранжаОпределение. Интерполяцией называется приближенное илиточное нахождение значений какой-либо величины по известнымотдельным значениям этой величины или значениям других величин, связанных с данной.Определение. Интерполяционным многочленом называетсямногочлен Ln(x) степени n, принимающий значение yi в узлах xi,i = 0, 1, 2, …, n (рис. 2.1.1).Рис.

2.1.1. Пример интерполяционного многочленаПример 1. Пусть n = 1. Интерполяционный многочлен проходит через точки (x0, y0), …, (x1,y1) и представляет собой прямуюлинию (рис. 2.1.2):L1 ( x) = y042x − x1x − x0; L1 ( x0 ) = y0 ; L1 ( x1 ) = y1.+ y1x0 − x1x1 − x0В общем случае интерполяционный многочлен n-й степенипроходит через n + 1 точку ( x0 , y0 ),( x1 , y1 ), ( x2 , y2 ),...,( xn , yn ), xi ≠ x jпри i ≠ j, и имеет видnni =0j =0j ≠iLn ( x ) = ∑ yi ∏( x − x0 )( x − x1 ) ( x − xi −1 )( x − xi +1 ) ( x − xn ).( xi − x0 )( xi − x1 ) ( xi − xi −1 )( xi − xi +1 ) ( xi − xn )Записанный в таком виде интерполяционный многочлен называют интерполяционным многочленом Лагранжа.

Интерполяционный многочлен существует и единствен.Рис. 2.1.2. Интерполяционный многочлен степени n = 1Пример 2. Пусть n = 2. Интерполяционный многочлен проходит через точки (x0, y0), (x1, y1), (x2, y2) и представляет собой параболу (рис. 2.1.3):L2 ( x) = y0( x − x1 )( x − x2 )( x − x0 )( x − x2 )( x − x0 )( x − x1 )+ y1+ y2;( x0 − x1 )( x0 − x2 )( x1 − x0 )( x1 − x2 )( x2 − x0 )( x2 − x1 )L2 ( x0 ) = y0 ;L2 ( x1 ) = y1 ;L2 ( x2 ) = y2 .43Рис. 2.1.3.

Интерполяционный многочлен степени n = 2Оценим погрешность интерполяционного многочлена Лагранжа.Пусть f ( x ) непрерывно дифференцируема n + 1 раз. Оценимразность f ( x) − Ln ( x) в фиксированной точке x ≠ xi , i = 0, 1, ..., n(рис. 1.2.4).Рис. 2.1.4. Оценка погрешности интерполяционногомногочлена Лагранжа в фиксированной точке x44Рассмотрим функциюϕ(t ) = f (t ) − Ln (t ) − K ωn (t ),nгде K — некоторая константа; ωn (t ) = ∏ (t − xi ).

Функция ϕ(t ) = 0i =0в точках x0 , x1 , x2 , ..., xn , так как f ( xi ) = Ln ( xi ), а ωn ( xi ) = 0,i = 0, 1, ..., n.Выберем константу K таким образом, чтобы ϕ( x) = 0. Тогдаf ( x ) − Ln ( x ) − K ωn ( x ) = 0 и ϕ( x) = 0 в n + 2 точках. Поэтому потеореме Ролля ϕ′(t ) = 0 в n + 1 точках (рис. 2.1.5).Рис. 2.1.5. Равенство нулю производной ϕ′(t) в n+1 точкеАналогично ϕ′′(t ) = 0 в n точках и т. д.

Наконец, найдется точка ξиз интервала (x0, xn), такая, что ϕn+1 (ξ) = f n+1 (ξ) − K (n + 1)! = 0. ОтсюдаK=f ( n+1) (ξ),( n + 1)!поэтому45f ( x ) − Ln ( x) =f ( n+1) (ξ)ωn ( x ).(n + 1)!Получили оценку погрешности интерполяционного многочлена Лагранжа. ПустьM n = max | f ( n ) (ξ) |,ξ∈[ x0 , xn ]тогда| f ( x ) − Ln ( x) | ≤Mn| ωn ( x ) | .(n + 1)!Недостатком интерполяции с помощью интерполяционногомногочлена Лагранжа является то, что при небольшом значении nинтерполяционный многочлен Лагранжа достаточно хорошо приближает гладкую функцию, а при большом значении n наблюдаются значительные колебания интерполяционного многочленамежду узлами интерполяции.2.2. Сплайн-интерполяцияОпределение.

Сплайном порядка m называется функция, являющаяся многочленом степени m на каждом из частичных интервалов [xi–1, xi], i = 0, 1, 2, …, n, принимающая заданные значения yiв узлах xi, i = 0, 1, 2, …, n и имеющая непрерывные производныедо (m – 1)-го порядка включительно (рис. 2.1.1).Рассмотрим более подробно случай, когда m = 3, — кубический сплайн.На каждом из частичных интервалов [xi–1, xi] сплайн будем искать в видеf ( x ) = ai + bi ( x − xi −1 ) + ci ( x − xi −1 ) 2 + di ( x − xi −1 )3 ,x ∈ [ xi−1 , xi ], i = 1, 2,..., n. (2.2.1)При этом из условия непрерывности в узлах сплайна его первых и вторых производных получаемf ( xi − 0) = f ( xi + 0);46(2.2.2)f ′( xi − 0) = f ′( xi + 0);(2.2.3)f ′′( xi − 0) = f ′′( xi + 0);(2.2.4)f ′′( x0 ) = f ′′( xn ).(2.2.5)Из (2.2.1) имеем f ( xi−1 ) = yi −1 = ai , отсюдаf ( xi ) = yi −1 + bi hi + ci hi2 + di hi3 , i = 1, 2, ..., n,где hi = xi − xi −1 , i = 1, 2, ..., n.В то же время, если рассматривать сплайн (2.2.1) на интервале[ xi , xi +1 ] , то f ( xi ) = yi = ai +1 .Отсюда и из (2.2.2) получаемyi = yi −1 + bi hi + ci hi2 + di hi3 , i = 1, 2, ..., n.(2.2.6)Найдем первую и вторую производные (2.2.1):f ′( x ) = bi + 2ci ( x − xi −1 ) + 3di ( x − xi −1 ) 2 ;f ′′( x ) = 2ci + 6di ( x − xi−1 ).Отсюда и из условий (2.2.3) и (2.2.4) следуетf ′( xi − 0) = bi + 2ci hi + 3di hi2 ;f ′( xi + 0) = bi +1;bi + 2ci hi + 3di hi2 = bi +1 ;(2.2.7)f ′′( xi − 0) = 2ci + 6di hi ;f ′′( xi + 0) = 2ci +1 ;2ci + 6di hi = 2ci +1.(2.2.8)Из (2.2.8) получаемdi =ci +1 − ci.3hi(2.2.9)47Из (2.2.6) и (2.2.9) следуетbi =yi − yi −1− ci hi − di hi2 .hiПосле подстановки di в последнее выражение получаемbi =yi − yi −1 hi− (2ci + ci +1 ).3hi(2.2.10)Подставим (2.2.10) в (2.2.7):yi − yi −1 hiy −y h− (2ci + ci +1 ) + 2ci hi + 3di hi2 = i+1 i − i +1 (2ci+1 + ci+ 2 ).33hihi+1Используя (2.2.9), находимyi − yi −1 hi3(c − c ) y − y h− (2ci + ci +1 ) + 2ci hi + i +1 2 i = i+1 i − i +1 (2ci +1 + ci+ 2 ).hi3(3hi )hihi+13После упрощений получим СЛАУ с трехдиагональной матрицей; эту СЛАУ можно решить методом прогонки:hi2hy − y y − yi −1⎛2⎞ci + ⎜ hi + hi +1 ⎟ ci +1 + i +1 ci+ 2 = i +1 i − i, i = 1, 2, ..., n − 1.333hi +1hi⎝3⎠(2.2.11)Из условия (2.2.5) и (2.2.8) имеем c1 = 0, cn+1 = 0.Условие (2.2.5) нужно было для единственности решения возникающей из (2.2.2) – (2.2.4) системы линейных уравнений.Если узлы равноотстоящие, то шаг h = (b − a ) / n , где [ a, b] —рассматриваемый отрезок, система (2.2.11) упрощается:h4hy − 2 yi + yi +1ci + h ci +1 + ci + 2 = i −1, i = 1, 2, ..., n − 1;333hc1 = 0, cn+1 = 0.Остальные коэффициенты сплайнов находят по формулам48ai = yi −1 ;bi =yi − yi −1 h− (2ci + ci+1 );3hdi =ci +1 − ci, i = 1, 2, ..., n.3hПусть sh ( x) — кубический сплайн, построенный для функцииf ( x ) на интервале [ a, b] с равноотстоящими узлами, т.

е. sh ( a + ih) == f ( a + ih), i = 0,1, ..., n. Тогда имеет место следующая теорема.Теорема. Для функции f ∈ C ( 4) [a, b] справедливы оценкиf ( x) − sh ( x)C [ a ,b ]≤ M 4h4 ;f ′′( x) − sh′′ ( x)f ′( x) − sh′ ( x )C [ a ,b ]C [ a ,b ]≤ M 4 h3 ;≤ M 4h2 .Из этих оценок следует, что при шаге h → 0 (т. е. при n → ∞)последовательности sh( i ) ( x ), i = 0, 1, 2, сходятся соответственно кфункциям f (i ) ( x), i = 0, 1, 2.Задание к лабораторной работе«Сплайн-интерполяция»1. Построить таблицу значений yi = f (a + ih) (табл. 2.2.1) наотрезке [ a, b] с шагом h = (b − a) / n.2. По полученной таблице вычислить коэффициенты сплайна,используя метод прогонки.3. Вычислить значения сплайна и заданной функции в серединах получившихся интервалов, т.

е. в точках xi = a + (i − 0,5) h,i = 1, 2, ..., n.4. Вычисления произвести при числе разбиений n = 5, 25, 125.Оформить таблицу, столбцами которой являются:1) значения функции xi = a + (i − 0,5) h, i = 1, 2, 3, 4, 5;492) значения заданной функции yi = f ( xi ), i = 1, 2, 3, 4, 5;3) значения сплайна при n = 5yiSpline5 = yi−1 +bi h ci h 2 di h3++, i = 1, 2, ..., n,248в серединах получившихся интервалов;4) значения сплайна при n = 25yiSpline 25 = yi−1 +bi h ci h 2 di h3++, i = 3 + 5 j, j = 0, 1, 2, 3, 4,248т.

е. в тех же точках, что и при n = 5;5) значения сплайна при n = 125yiSpline125 = yi −1 +bi h ci h2 di h3++, i = 13 + 25 j, j = 0, 1, 2, ..., 24,248т. е. в тех же точках, что и при n = 5.Убедиться, что при увеличении числа разбиений n качествосплайн-интерполяции повышается.Таблица 2.2.1Варианты функцийВариант123456789101150Функцияe + sin xln(2 x − 1) − sin x 2arctg(2 x + 3)x3x + 2 + tg xπ⎞⎛sin ⎜ 2 x − ⎟3⎠⎝π⎞⎛cos ⎜ 2 x + ⎟3⎠⎝arcsin(2 x − 1)sh xch xth xe x − e2 xИнтервал[0, 2][1, 3][–1, 3][–1, 1][0, π][0, π][0, 1][0, 2][0, 2][0, 2][–1, 1]Окончание табл.

2.2.1Вариант121314ФункцияИнтервалπ⎞⎛sin ⎜ x + ⎟6⎠⎝ln(2 x + 1) + 2 sin 3 xe − cos 2 x2x[0, π][0, 3][0, 3]15π⎞⎛ln(2 x + 1) + sin ⎜ x + ⎟3⎠⎝[0, 5]162x − 1π⎞⎛+ sin ⎜ 4 x − ⎟3x + 13⎠⎝[0, 3]17181920π⎞⎛2 x + 1 − sin ⎜ 5 x − ⎟6⎠⎝π⎞⎛sh x − cos ⎜ 5 x − ⎟3⎠⎝π⎞⎛9 x − 2 + sin ⎜ 3 x + ⎟3⎠⎝arctg(2 x + 3) + cos x 2[0, 5][–1, 4][1, 6][–2, 3]21π⎞⎛e 2 x + sin ⎜ 5 x − ⎟4⎝⎠[–3, 3]22sin(cos x + 3)[–3, 3]2324252627282930[0, 1][0, 5]5xsin 2x sin x 2sin 5xxcos(sin x 2 )[0, 2][0, 3]x + 2 + cos(sin x )2π⎞⎛2 x − 1 + cos ⎜ 5 x + ⎟6⎠⎝arctg(2 x − 1) − sin x 22x + 1π⎞⎛+ cos ⎜ 4 x + ⎟3x − 13⎠⎝[0, 3][1, 4][0, 3][1, 3]512.3. Метод наименьших квадратовПусть известны значения yi в узлах xi, i = 0, 1, 2, ..., n;ϕ(a0, a1, ..., am, x) — функция, зависящая от параметров a0, a1, ...,am. Рассмотрим функцию S:nni =0i =0S = ∑ (ϕ( a0 , a1 , ..., am , xi ) − yi ) 2 = ∑ εi2 .Выберем параметры a0 , a1 , a2 , ..., am так, чтобы минимизироватьзначение S, т.

е. сумму квадратов невязок εi2 (рис. 2.3.1).Рис. 2.3.1. Минимизация суммы квадратов невязокПолучим систему уравнений∂S= 0, i = 0,1, 2, ..., m.∂aiЭту систему уравнений (часто нелинейную) можно решить методом Ньютона.52Рассмотрим подробнее случай, когда функция ϕ(x) = am xm ++ am–1xm–1 + ... + a1x1 + a0. Условие∂S= 0, i = 0,1, 2, ..., m .∂aiприводит к следующей СЛАУ:n∂S= 2∑ ( am xim + am−1 xim−1 + ... +a1 xi1 + a0 − yi )1 = 0;∂a0i =0n∂S= 2∑ ( am xim + am−1 xim−1 + ... +a1 xi1 + a0 − yi ) xi = 0;∂a1i =0n∂S= 2∑ ( am xim + am−1 xim−1 + ...

+a1 xi1 + a0 − yi ) xi2 = 0;∂a2i =0...n∂S= 2∑ (am xim + am−1 xim−1 + ... +a1 x1i + a0 − yi ) xim = 0.∂ami= 0Преобразуем ее к видуn⎛ n ⎞⎛ n⎞⎛ n⎞(n + 1) a0 + ⎜ ∑ xi ⎟ a1 + ⎜ ∑ xi2 ⎟ a2 + ... + ⎜ ∑ xim ⎟ am = ∑ yi ;i =0⎝ i =0 ⎠⎝ i =0 ⎠⎝ i =0 ⎠n⎛ n ⎞⎛ n 2⎞⎛ n 3⎞⎛ n m+1 ⎞⎜ ∑ xi ⎟ a0 + ⎜ ∑ xi ⎟ a1 + ⎜ ∑ xi ⎟ a2 + ... + ⎜ ∑ xi ⎟ am = ∑ xi yi ;i =0⎝ i =0 ⎠⎝ i =0 ⎠⎝ i =0 ⎠⎝ i =0⎠…n⎛ n m⎞⎛ n m+1 ⎞⎛ n m+2 ⎞⎛ n 2m ⎞m⎜ ∑ xi ⎟ a0 + ⎜ ∑ xi ⎟ a1 + ⎜ ∑ xi ⎟ a2 + ... + ⎜ ∑ xi ⎟ am = ∑ xi yi .i =0⎝ i =0 ⎠⎝ i =0⎠⎝ i =0⎠⎝ i =0⎠nni =0i =0Введем коэффициенты b pq = ∑ xip+ q , c p = ∑ xip yi .Получим СЛАУ:b00 a0 + b01a1 + ...

+ b0 m am = c0 ;b10 a0 + b11a1 + ... + b1m am = c1;53…bm 0 a0 + bm1a1 + ... + bmm am = cm .Эту СЛАУ решаем методом Гаусса.Пример. Пусть даны точки:x–1012y1–1141. Найдем методом наименьших квадратов прямую ϕ( x) == a0 + a1 x, на которой минимизируется сумма квадратов невязок.Получим систему уравнений4a0 + 2a1 = 5;2a0 + 6a1 = 8.Отсюда ϕ( x) = 0,7 + 1,1x (см. рис. 2.3.2).2. Найдем методом наименьших квадратов параболу ψ ( x) == a0 + a1 x + a2 x 2 , на которой минимизируется сумма квадратов невязок. Получаем систему уравнений4a0 + 2a1 + 6a2 = 5;2a0 + 6a1 + 8a2 = 8;6a0 + 8a1 + 18a2 = 18.Решив ее, найдем, что (рис. 2.3.3) ψ ( x) = 1, 25 x 2 − 0,15 x − 0,55.54Рис. 2.3.2.