интерполяции Прежде чем переходить к определенным типам интерполянтов, сделаем два замечания. Во-первых, сами по себе данные таблицы 2.1 не могут определить конкретный вид интерполирующей функции. Например, точки графика на рис. 2.1 можно последовательно соединить отрезками прямых. С другой стороны, всегда можно найти алгебраический полином n -й степени с действительными коэффициентами, график которого точно проходит через n + 1 заданных точек, если все узлы таблицы 2.1 различны. Вообще существует бесконечное количество функций F(x), удовлетворяющих условию (2.1). В качестве интерполирующих функций на практике часто используют алгебраические полиномы, суммы экспонент, Фурье- суммы, сплайны и т.д. На выбор типа интерполянта в конкретной задаче влияет любая дополнительная информация о зависимости между переменными x и y. Второе замечание базируется на свойстве интерполянта (2.1) принимать в узлах точные значения интерполируемой функции. Если значения функции y i таблицы 2.1 имеют существенную погрешность, то точное выполнение требования (2.1) является неосновательным. В таких ситуациях при вычислении значений функции f(x) для аргументов x, не совпадающих с узлами, целесообразно применять методы аппроксимации, рассматриваемые в следующей главе. 14 § 2.2. Интерполяционный полином Лагранжа Алгебраические полиномы с действительными коэффициентами являются хорошо изученными функциями и весьма простыми для вычислений. Их удобно складывать и перемножать, дифференцировать и интегрировать. Алгебраический полином степени n можно записать в стандартной форме (1.4), хотя допустимы и другие представления. Из-за простоты работы с алгебраическими полиномами их часто используют в качестве интерполянтов. Пусть исходной информацией являются ( n + 1) пар чисел xi, yi ( i = 0, 1, ..., n), сведенных в таблицу вида 2.1. Будем строить интерполирующую функцию в виде алгебраического полинома. Из теории функций известно, что через n + 1 произвольных точек на плоскости, заданных координатами xi, yi ( i = 0, 1, .. . , n), всегда можно провести график алгебраического полинома с действительными коэффициентами, если числа xi ( i = 0, 1, .. . , n) различны. Если степень полинома m = n, то, вообще говоря, такой полином является единственным. Согласно требованию (2.1), в узлах интерполяции xi ( i = 0, 1, ... , n) значения этого полинома Pm( x) точно совпадают с данными значениями yi ( i= 0, 1, ... , n) интерполируемой функции f( x): Pm( xi) = yi( i = 0, 1, ... , n). (2.2) Таким образом, всегда можно найти полином Pm( x) степени m = n, который может служить интерполирующей функцией. В специальных случаях, когда точки xi , yi ( i = 0, 1, ..., n) расположены особым образом, степень полинома может быть меньше n. Например, если данные точки расположены точно на одной прямой, то через них пройдет график линейной функции (полинома первой степени) независимо от количества данных точек. При m > n полиномов, графики которых проходят через заданные точки, бесконечно много. Рассмотрим ситуацию, когда расстояния между узлами неодинаковы. В таком случае целесообразно построить интерполяционный полином Лагранжа . Этот полином строится в виде суммы: 15 Ln( x) = iniiyx⋅ ∑ = ) ( 0 l , (2.3) где il ( x) − функции аргумента x, называемые коэффициентами Лагранжа ( i = 0, 1, ..., n). По свойству интерполирующего полинома (2.1) должно выполняться условие: Ln( xi) = yi( i = 0, 1, ..., n). (2.4) Для выполнения (2.4) на коэффициенты Лагранжа налагается следующее требование: il ( xj ) = δ i j ( i, j = 0, 1, ..., n), (2.5) где δ i j – символ Кронекера ( δ i j = 0 при i ≠ j и δ i j = 1 при i = j ). Чтобы обеспечить выполнение равенств (2.5), коэффициенты Лагранжа il (x) представим в виде произведений: il ( x) = с i ( x − xo ) ( x − x1 ) ... ( x − xi-1 ) ( x − xi+1 ) ... ( x − xn), (2.6) где с i − некоторые постоянные величины, пока неизвестные. В выражении для i-го коэффициента Лагранжа отсутствует множитель ( x − xi), поэтому каждое i-е произведение (2.6) содержит n скобок, и каждый коэффициент Лагранжа представляет собой полином степени n. Для определения явного вида постоянных величин с i можно воспользоваться равенством il ( xi) = 1 ( i = 0, 1, ..., n), согласно требованию (2.5). Тогда из (2.6) сразу следует, что 1 / с i = ∏ ≠ = − nijjjixx, 0 ) ( Таким образом, коэффициенты Лагранжа представляются произведениями отношений разностей: il ( x) = ∏ ≠ = − − nijjjijxxxx, 0 ) ( ) ( (2.7) 16 Легко убедиться непосредственной подстановкой, что для каждого коэффициента (2.7) выполняется условие (2.5). Интерполяционный полином в форме Лагранжа получается подстановкой найденных коэффициентов (2.7) в выражение (2.3): Ln( x) = ∑ ∏ = ≠ = − − ⋅ ninijjjijixxxxy0 , 0 ) ( ) ( (2.8) Для большей наглядности запишем полином Лагранжа 3-й степени в развернутом виде: L3 ( x) = y0 ) )( )( ( ) )( )( ( 3 0 2 0 1 0 3 2 1 xх xх xх xх xх xх − − − − − − + y1 ) )( )( ( ) )( )( ( 3 1 2 1 0 1 3 2 0 xх xх xх xх xх xх − − − − − − + + y2 ) )( )( ( ) )( )( ( 3 2 1 2 0 2 3 1 0 xх xх xх xх xх xх − − − − − − + y3 ) )( )( ( ) )( )( ( 2 3 1 3 0 3 2 1 0 xх xх xх xх xх xх − − − − − − Использование полинома Лагранжа для интерполяции иллюстрируется следующим примером. Пример 1. Пусть значения функции заданы таблицей 2.2. Требуется построить по ним полином Лагранжа. Таблица 2.2 Так как таблица 2.2 содержит 5 узлов, интерполяционный полином Лагранжа должен иметь 4-ю степень. Подставляя эти числовые значения таблицы в общую формулу (2.8), находим явный вид полинома Лагранжа для данной задачи: L4 ( x) = −4 ) 5 2 )( 5 , 3 2 )( 2 , 1 2 )( 1 2 ( ) 5 )( 5 , 3 )( 2 , 1 )( 1 ( − − − − − − + − − − − + х х х х + i 0 1 2 3 4 xi −2 −1 1,2 3,5 5 yi −4 0 3,8 −2 7 17 + 3,8 ) 5 2 , 1 )( 5 , 3 2 , 1 )( 1 2 , 1 )( 2 2 , 1 ( ) 5 )( 5 , 3 )( 1 )( 2 ( − − + + − − + + х х х х − −2 ) 5 5 , 3 )( 2 , 1 5 , 3 )( 1 5 , 3 )( 2 5 , 3 ( ) 5 )( 2 , 1 )( 1 )( 2 ( − − + + − − + + xxxx+ + 7 ) 5 , 3 5 )( 2 , 1 5 )( 1 5 )( 2 5 ( ) 5 , 3 )( 2 , 1 )( 1 )( 2 ( − − + + − − + + xxxx Сумма содержит только 4 слагаемых, так как коэффициент Лагранжа 1 l ( x) умножается на нулевое значение y1 = 0. Приведение полинома к стандартному виду даст следующее выражение: L4 ( x) = 0,0819538 x4 − 0,181866 x3 − 1,43424 x2 + 2,19964 x + 3,37007. График построенного полинома приведен на рис. 2.2. -2 0 2 4 6 -4 -2 0 2 4 6 8 yxРис. 2.2. Интерполяционный полином Лагранжа, проведенный через 5 заданных точек таблицы 2.2. Как было указано в начале данного параграфа, полином Лагранжа строится для произвольного набора узлов. Но вычисление значения полинома Лагранжа для заданного аргумента x по формуле (2.8) требует большого количества вычислений. Можно подсчитать, что при этом 18 необходимо выполнить 2n (n + 1) + n сложений и вычитаний, а также 2n (n + 1) умножений и делений. В подробных курсах численных методов доказывается, что полином Лагранжа (2.8) обладает свойством единственности [1]. Полиномы, интерполирующие таблицы данных с постоянным шагом, позволяют проводить более быстрые вычисления, сравнительно с использованием универсального полинома Лагранжа. § 2.3. Интерполяция по равноотстоящим узлам Рассмотрим случай, когда известны значения интерполируемой функции f(x) в равноотстоящих узлах. При этом узлы интерполяции x i выражаются формулой: x i = x 0 + i ⋅ h , где i = 0, 1, ..., n . (2.9) Постоянный параметр h называется шагом интерполяции. В такой задаче исходными данными являются (n + 3) чисел: начальный узел интерполяции x 0 , шаг интерполяции h и (n + 1) значений неизвестной функции y i (i = 0, 1, ..., n) в узлах интерполяции. Будем строить полином P n (x) степени n, обладающий свойством (2.2) P n (x i ) = y i (i = 0, 1, ..., n), (2.10) где значения x i заданы формулой (2.9). Полиномы для интерполяции по равноотстоящим узлам впервые были построены Ньютоном, который был как гениальным физиком, так и гениальным математиком. Первый полином Ньютона строится в форме: P n (x) = a 0 + ∑ ∏ = − = − ⋅ n i i j j i x x a 1 1 0 ) ( (2.11) или в развернутом виде: P n (x) = a 0 + a 1 (x − x 0 ) + a 2 (x − x 0 ) (x − x 1 ) + a 3 (x − x 0 ) (x − x 1 ) (x − x 2 ) +... + a n (x − x 0 )...(x − x n- 1 ).
19 Коэффициенты полинома a i необходимо выразить через известные величины x 0 , h, n, y i (i = 0, 1, . . . n), пользуясь требованием (2.10), т.е. прохождением графика полинома через все точки заданной системы. Сразу заметим, что подстановка x = x 0 в выражение (2.11) дает P n (x 0 )=a 0 . Следовательно, согласно (2.10), мы получаем значение свободного члена искомого полинома: a 0 = y 0 (2.12) Остальные коэффициенты полинома Ньютона выразим через конечные разности, описанные в предыдущей главе. Можно доказать эквивалентность условия (2.10) тому, что ∆ k P n (x 0 ) = ∆ k y 0 (k = 0, 1, ... , n). (2.13) Доказательство основано на том, что согласно (1.15), n-я конечная разность любой функции f(x) выражается через n последовательных значений этой функции в равноотстоящих точках, отличающихся на шаг h. Будем строить последовательно конечные разности ∆ i P n (x), где (i = 0, 1, ... , n), для произвольной точки интерполяции x и заданного шага h. Для нахождения первого коэффициента a 1 составим первую конечную разность полинома (2.11): ∆ 1 P n (x) = P n (x+ h) − P n (x). (2.14) Пользуясь свойством линейности конечных разностей, будем их вычислять для отдельных слагаемых суммы (2.11) и выносить за скобки общие множители a i (i = 0, 1, ... , n): i = 0 a 0 − a 0 = 0 i = 1 a 1 [(x + h − x 0 ) − (x − x 0 )] = a 1 h i = 2 a 2 [(x+ h − x 0 ) (x + h − x 1 ) − (x − x 0 ) (x − x 1 )] = = a 2 [(x + h − x 0 ) (x − x 0 ) − (x − x 0 ) (x − x 0 − h)] = a 2 2h(x − x 0 ). В последнем преобразовании использовалось, что x 1 = x 0 + h. Дальнейшими вычислениями можно показать, что среди сомножителей остальных коэффициентов a i (i > 2) обязательно
20 присутствуют множители ( x − x0 ). В результате получим следующее для первой конечной разности исходного полинома: ∆ 1 Pn( x) = a1 h + 2 a2 ( x − x0) h + 3 a3 ( x − x0 ) ( x − x1 ) h +... + n an ( x − x0 )...( x − xn− 2) h. (2.15) Теперь если в (2.15) положить x = x0 , получим: ∆ 1 Pn( x0 ) = a1 h. С другой стороны, из (2.13) и (2.10) следует: ∆ 1 Pn( x0 ) = y1 − y0 . Сравнивая два выражения для ∆ 1 Pn( x0 ), получаем представление первого коэффициента искомого полинома через известные величины: a1 = ( y1 − y0 ) / h. (2.16) Для вычисления коэффициента a2 рассматривается вторая конечная разность полинома (2.11): ∆ 2 Pn( x) = ∆ 1 Pn( x + h) − ∆ 1 Pn( x) . Ясно, что она не будет содержать слагаемого с коэффициентом a1 , т.к. этот коэффициент является множителем константы в сумме ∆ 1 Pn( x), согласно (2.15). Если провести преобразования, аналогичные проделанным для первой конечной разности исходного полинома, то мы получим, что свободным членом нового полинома ∆ 2 Pn( x) будет 2 a2 h2 , а остальные члены будут содержать сомножители ( x − x0 ). Тогда ясно, что в точке x = x0 получим: ∆ 2 Pn( x0 ) = 2 a2 h2 С другой стороны, по формуле (1.3) для произвольной функции имеем: ∆ 2 Pn( x0 ) = Pn( x2 ) − 2 Pn( x1 ) + Pn( x0 ), а по условию (2.10): ∆ 2 Pn( x0 ) = y2 − 2 y1 − y0 Сравнение дает выражение для второго коэффициента искомого полинома: a2 = ( y2 − 2 y1 − y0 ) / (2 h2 ). (2.17) 21 Для получения остальных коэффициентов искомого полинома (2.11) необходимо повторять аналогичную процедуру с высшими конечными разностями. В результате мы получим следующее общее выражение для коэффициентов интерполянта Ньютона: ai = iihiy! ) ( 0 ∆ (2.18) Следовательно, интерполяционный полином Ньютона выразится через конечные разности порядков 1, ... , n, вычисленные в нулевом узле x0 : Pn( x) = y0 + ∑ ∏ = − = − ∆ niikkiixxhiy0 1 0 0 ) ( ! ) ( (2.19) То, что построенный полином удовлетворяет условию (2.10), проверяется непосредственной подстановкой. Положим x = xj, где xj – произвольный узел, и подставим в (2.19). При этом учтем, что все разности ( xj− xk) кратны целому числу шагов h при k ≠ j и равны нулю при k = j. Следовательно, все слагаемые суммы (2.19) с индексами i > j содержат множитель, равный нулю. В результате получим: Pn( xj) = y0 + j ∆ y0 + 2 ) 1 ( − ⋅ jj ∆ 2 y0 + … + ! 1 ) 1 ( jjj⋅ ⋅ − ⋅ ∆ jy0. Использование (1.14) и (2.13) позволяет представить правую часть последнего уравнения в виде Pn( x0 + jh) = Pn( xj) = yj, что и требовалось доказать. Для практических расчетов полином Ньютона (2.19) удобнее представлять, вводя новую переменную: q = ( x − x0 ) / h. (2.20) Тогда интерполяционный полином Ньютона представится так: 22 Pn( x) = y0 + q⋅∆ 1 y0 + ! 2 ) 1 ( − qq⋅ ∆ 2 y0 + ! 3 ) 2 )( 1 ( − − qqq⋅ ∆ 3 y0 +... + ! ) 1 )...( 1 ( nnqqq+ − − ⋅ ∆ n y0 (2.21) Анализ формулы (2.21) показывает, что погрешность построенного интерполянта Ньютона принимает наименьшие значения для аргументов, близких к нулевому узлу x0 (см. § 2.5). Пример, приведенный в предыдущем параграфе для полинома Лагранжа, показывает, что не всегда целесообразно строить полином по всем узлам, так при этом могут возникать необоснованные «выбросы». В практических расчетах часто заранее выбирается степень интерполирующего полинома. Затем в качестве нулевого x0 берется узел, ближайший слева к значению аргумента x, для которого требуется вычислить значение неизвестной функции f ( x). Интерполянт n-й степени (2.21) строится с помощью значений функции yi в узле x0 и в n узлах, расположенных правее аргумента x. Значения функции в узлах левее x0 не используются. Пример 2. Исходные данные сведены в таблицу 2.3 с постоянным шагом h = 0,1. Таблица 2.3 xi 3,0 3,1 3,2 3,3 3,4 3,5 3,6 yi 0,550 0,450 0,365 0,293 0,234 0,185 0,146 Пусть требуется вычислить значение функции f( x) для аргумента x=3,22. Узлы в таблице 2.3 являются равноотстоящими (шаг h = 0,1), поэтому целесообразно построить интерполянт Ньютона. Ограничимся полиномом 3-го порядка. За нулевой узел, как было рекомендовано ранее, возьмем ближайший слева к значению x = 3,22 , т.е. узел x0 = 3,2. Полином будем строить по четырем узлам xi≥ 3,2. Составим таблицу конечных разностей ∆ kyi ( k =1, 2, 3), необходимых для построения полинома (2.21) при n = 3. Таблицу построим по типу примера 2 предыдущей главы (см. табл. 1.1). Напомним, что разности порядка выше первого вычисляются как разности разностей предыдущего порядка, стоящие в соседнем левом столбце таблицы 2.4. 23 Таблица 2.4 i x i y i ∆y i ∆ 2 y i ∆ 3 y i 0 3,2 0,365 − 0,072 0,013 − 0,003 1 3,3 0,293 − 0,059 0,01 2 3,4 0,234 − 0,049 3 3,5 0,185 Для значений x = 3,22 и x 0 = 3,2 вычислим параметр q = (3,22–3,2)/0,1 = 0,2 и подставим все требуемые числа в формулу (2.21): f (3,22) = P 3 (x = 3,22) = 0,365 + 0,2 ( −0,072) + 2 ) 1 2 , 0 ( 2 , 0 − 0,013 + + ! 3 ) 2 2 , 0 )( 1 2 , 0 ( 2 , 0 − − ( −0,003) ≈ 0,349. Исходные значения функции в узлах интерполирования были заданы с тремя значащими цифрами (см. табл.2.3), поэтому, согласно правилам, результат округляем до трех значащих цифр. Приведем формулы для двух важных частных случаев интерполяции. 1. Линейная интерполяция. Интерполяционный полином (2.21) строится по двум узлам и является линейной функцией: P 1 (x) = y 0 + q ⋅∆ 1 y 0 или P 1 (x) = y 0 + h x x 0 − (y 1 − y 0 ). (2.22) Последнюю формулу легко получить с помощью рисунка 2.1, соединив точки (x 0 , y 0 ) и (x 1 , y 1 ) отрезком прямой линии. Соединяя отрезками прямых точки (x i , y i ) и (x i+ 1 , y i+ 1 ) для соседних узлов, можно получить график кусочно-линейной интерполяции. Практика показывает, что для большинства физических задач такой вид интерполяции неприменим, так как дает большую погрешность. 2. Квадратичная интерполяция. Для построения интерполяционного полинома (2.21) выбираются три узла. Полином становится квадратичной функций следующего вида:
24 P2 ( x) = y0 + q⋅ ∆ 1 y0 + 2 ) 1 ( − qq⋅ ∆ 2 y0 или P2 ( x) = y0 + q⋅ ( y1 − y0 ) + 2 ) 1 ( − qq⋅ ( y2 − 2 y1 + y0 ). (2.23) Выше было показано, что интерполянт (2.21) строится с использованием n значений функции в узлах данной таблицы, которые расположены правее аргумента x, для которого проводится интерполяция. Если же аргумент x находится ближе к концу таблицы, например, лежит в интервале ( xn− 1 , xn), то построить полином вида (2.21) невозможно (за исключением линейного). В подобных случаях следует вместо полинома (2.21) использовать второй Поделитесь с Вашими друзьями: |