Интерполяция функций одной и нескольких переменных

1. ИНТЕРПОЛЯЦИЯ ФУНКЦИЙ ОДНОЙ И НЕСКОЛЬКИХ ПЕРЕМЕННЫХ В задачах математического моделирования, которые для своей реализации требуют, как правило, применения ЭВМ, из различных способов задания функций наиболее часто используется табличный способ. При этом возникает необходимость определения значений функции для значений аргумента, отличающихся от тех, которые содержатся в таблице. В этих целях применяют процедуру интерполяции, а также экстраполяции. По сути, речь идет о восстановлении непрерывной функции по заданному набору ее значений для счетного множества значений аргумента. Задача восстановления функции может ставиться и иначе, когда ищется функция, близкая в определенном смысле к таблично заданной функции (например, наилучшее среднеквадратичное приближение). В обоих случаях исходную функциональную зависимость y( x ) , по которой сопоставлена таблица, заменяют приближенной функцией  ( x , a ) , которую и используют в дальнейших вычисления. При этом близость y( x ) и  ( x , a ) обеспечивается подбором свободных параметров a  a0 , a1 ,..., a n . В простейшем варианте лагранжевой интерполяции, когда функция y( x ) задана в узлах некоторой сетки x i , параметры a определяют из условия совпадения  ( x , a ) со зна- чениями функции в фиксированном числе узлов. Получаемая при этом система уравнений имеет вид ( xi , a0 , a1 ,..., a n )  y( xi )  yi ,0  i  n. (2.1) Из этой системы можно определить все ai . Интерполяция может быть линейной или нелинейной, в соответствии с характером зависимости функции  ( x , a ) от параметров a . 1.1. Линейная интерполяция При линейной интерполяции функция  ( x , a ) имеет вид n ( x , a )   a j  j ( x ), j 0 (2.2) где все функции  j ( x ) линейно независимы. Подставляя (2.2) в (2.1), получим систему линейных уравнений для определения ai : n a  j 0 j j ( xi )  y i , 0  i  n. (2.3) Единственность решения обеспечивается требованием неравенства нулю определителя системы (1.3):  0 ( x0 )  1 ( x0 )...  n ( x0 )  . ... . 0 (2.4)  0 ( x n )  1 ( x n )...  n ( x n ) при x k  xl (т.е. при любых несовпадающих узлах). Систему функций, удовлетворяющую условию (2.4), называют чебышевской. Из различных систем функций наиболее распространены многочлены, хотя применяют также тригонометрические и экспоненциальные функции. Если в качестве системы функций выбрать степенные функции аргумента, т.е.  j ( x )  x j , то определитель (2.4) окажется определителем Вандермонда, который не равен нулю при условии x k  xl . Следовательно, интерполяционный полином всегда существует и он единствен. Для практических вычислений удобно использовать многочлен Pn ( x ) в форме интерполяционного полинома Ньютона. Введем понятие разделенных разностей функции y( x ) , заданной в узлах x i : y( x i , x j )  [ y( x i )  y( x j )] /( x i  x j ), y( x i , x j , x k )  [ y( x i , x j )  y( x j , x k )] /( x i  x k ), (2.5) y( x i , x j , x k , x l )  [ y( x i , x j , x k )  y( x j , x k , x l )] /( x i  x l ), Правило образования таких конструкций понятно из приведенной записи. Разделенные разности первого, второго и более высоких порядков связаны с производными соответствующих порядков. 2 Полином n-ой степени может быть представлен через разделенные разности в виде: Pn ( x )  Pn ( x0 )  ( x  x0 )P( x0 , x1 )  ( x  x0 )( x  x1 )P( x0 , x1 , x 2 ) ...  ( x  x0 )( x  x1 )...( x  x n 1 )P( x0 , x1 , x 2 ,..., x n ) (2.6) Ввиду того, что значения интерполяционного полинома в этих узлах совпадают со значениями интерполируемой функции, разделенные разности, входящие в (2.6), выражаются через узловые начения функции. В результате получается полином, называемый интерполяционным многочленом Ньютона: n y( x )  y( x0 )   ( x  x0 )( x  x1 )...( x  x k 1 ) y( x0 ,..., x k ). (2.7) k 1 При вычислениях по этой формуле точность расчетов удобно оценивать, наблюдая за тем, насколько быстро убывают члены ряда. Если это происходит достаточно быстро, можно оставлять только те члены, которые больше заданной погрешности расчетов. Погрешность многочлена Ньютона можно оценить по формуле y( x )  Pn ( x )  где M n 1  max y  n 1 M n 1  n( x ), ( n  1 )! (  ) - максимальное значение производной интерполируемой функ- ции на отрезке между наименьшим и наибольшим из значений x0 , x1 , x 2 ,..., x n , а полином n  n ( x )   ( x  xi ) . i 0 На равномерной сетке узлов приведенная формула трансформируется к виду y( x )  Pn ( x )  h M n 1    2 2n n 1 , (2.8) ,где h  xi 1  xi  xi  xi 1     (предполагается, что шаг сетки постоянен). Выражение (2.8) свидетельствует, что с уменьшением расстояния между узлами (шага) h погрешность представления функций полиномом Ньютона убывает как O( h n 1 ). Трудность использования указанных теоретических оценок на практике состоит в том, что производные интерполируемой функции обычно неизвестны, поэтому для определения погрешности удобнее воспользоваться оценкой первого отброшенного члена. Пример. Построить интерполяционный полином Ньютона четвертой степени для функции y( x )  cos( 2 x ) в области значений аргумента 0  x  / 4 . Заполним таблицу разделенных разностей, вычисляемых по пяти узлам, представив для удобства вычислений y( z )  cos[(  / 2 ) z ] , 0  z  1. z0 y( z 0 ) y( z 0 , z 1 ) y( z 0 , z 1 , z 2 ) y( z 0 ,...,z 3 ) y( z 0 ,...,z 4 ) 1 -0,304 0,25 0,924 -1,128 -0,868 0,5 0,707 0,363 -0,856 -1,296 0,75 0,383 0,149 0,512 -0,472 -1,532 1 Полином Ньютона y ( z )  1  0,304 z  1,128 z ( z  0,25)  0,363z ( z  0,25)( z  0,5)   0,149 z ( z  0,25)( z  0,5)( z  0,75) Вычислим y( 0 ,6 )  1  0 ,182  0 ,237  7 ,623  10 3  4 ,694  10 3  0 ,589. Точное значение y( 0,6 )  0,588 , т.е. точность вычислений по приближенной формуле оказалась весьма высокой. Опишем алгоритм построения полинома Ньютона. В качестве исходных данных задаются: таблично интерполируемая функция, значение аргумента x, для которого выполняется интерполяция, и степень полинома. 1. Формируется конфигурация из (n+1)-го узла, по возможности симметрично относительно значения x (понятно, что вблизи краев таблицы это невозможно). 2. На сформированной конфигурации узлов строится таблица, аналогичная той, которая представлены в вышеприведенном примере. 3. Строится полином Ньютона с разделенными разностями из верхней строки таблицы. Обратная интерполяция. Применяется для приближенного нахождения корня монотонных функций. Суть ее состоит в том, что меняются местами столбцы табличной функции и выполняется обычная интерполяция при задании x=0. Факт наличия корня у функции устанавливается по наличию смены знака у функции при продвижении по строкам таблицы. Помимо полинома Ньютона существует еще одна важная форма интерполяционного многочлена- полином в форме Лагранжа n Ln ( x )   L(i n ) ( x ) y( xi ) , i 0 (n) где лагранжевы коэффициенты Li ( x ) определяются по формуле ( x  xk ) . k 0 ( x i  x k ) n L(i n ) ( x )   k i Полином Лагранжа в отличие от полинома Ньютона содержит в явном виде значения интерполируемой функции в узлах y( x i ) . В заключение отметим, что в практике вычислений для интерполяции полиномы степени выше пятой обычно не используют, т.е. число узлов интерполяции не превышает шести. Если при таком числе узлов не обеспечивается заданная погрешность, следует уменьшать расстояние между узлами. 1.2. Многомерная интерполяция В различных приложениях широко используют двумерные и трехмерные таблицы. Например, теплофизические свойства различных веществ зависят от температуры и давления, а оптические характеристики - еще и от длины волны излучения. Осуществляя многомерную интерполяцию, следует помнить, что расположение узлов не может быть произвольным. Например, при интерполяции полиномом первой степени P1 ( x , y ) узлы не должны лежать на одной прямой в плоскости.. Проверять условия подобного типа достаточно сложно, поэтому на практике целесообразно строить регулярные сетки, как правило, прямоугольные и равномерные, когда узлы являются точками пересечения двух взаимно перпендикулярных систем параллельных прямых. На этой сетке проводят простую последовательную интерполяцию: сначала по строкам, а затем по столбцам.