Нелинейная интерполяция

1.3. Нелинейная интерполяция Для табулирования быстроменяющихся функций требуется весьма малый шаг, т.е. возникает необходимость создавать таблицы очень больших объемов, что в ряде случаев неприемлемо. Оказывается, что преобразованием переменных   ( y ) и добиться того, чтобы в новых переменных график    ( x ) можно  (  ) был близок к прямой хотя бы на отдельных участках. В этом случае интерполяцию проводят в переменных ( , ) , а затем обратным интерполированием находят y i  y(  i ) . Преобразования  ( y ) и  ( x ) должны быть достаточно простыми (логарифмиче- ская, экспоненциальная, тригонометрические и некоторые другие функции). При этом надо заботиться о том, чтобы и обратное преобразование y(  ) оказалось несложным. Во многих задачах теплофизики, гидродинамики, оптики, переноса излучения) и других областей науки и техники часто встречается степенная зависимость функции от своих аргументов. В этом случае удобны преобразования типа логарифмирования. Вообще говоря, в каждом конкретном случае приходится специально подбирать вид функций  ( y )и   ( x ) . Например. функция y  a может быть представb  c  e dx лена следующим образом ln( a a  b )  d  x  ln c , т.е.   ln(  b ),   x . y y Для функциональной зависимости y  переменных   x линеаризация достигается заменой ax  b y 1 1 ,   , а для функции y  axe  bx -   ln( ),   x . x y x Пример. Получить формулу для нелинейной двухточечной интерполяции функции   ln( y ) и y( x ) , если переменные можно преобразовать по формулам   x. Составим интерполяционный полином Ньютона на двухточечном шаблоне:   0  1  0 (    0 ). 1  0 В исходных переменных имеем ln( y )  ln( y 0 )  ln( y 1 )  ln( y 0 ) ( x  x 0 ), x1  x0 и окончательно y  y 0 ( y 1 / y 0 )( x  x ) /( x1  x0 ) . 1.4. Интерполяция сплайнами Слово “сплайн” переводится как “гибкая линейка”. Такую линейку можно использовать для проведения кривых через заданную совокупность точек, изгибая и придерживая ее так, чтобы ребро проходило через все точки на плоскости. Равновесие гибкой линейки описывается уравнением  IV ( x )  0 , т.е. интерполяционный полином на участке между каж- дой парой соседних точек имеет третью степень:  ( x )  a i  bi ( x  xi 1 )  ci ( x  xi 1 ) 2  d i ( x  xi 1 ) 3 , (3.1) xi 1  x  xi ,0  i  N . В узлах значения многочлена и интерполируемой функции совпадают:  ( x i 1 )  y i 1  a i , (3.2)  ( xi )  y i  ai  bi hi  ci hi2  d i hi3 , (3.3) hi  xi  xi 1 ,1  i  N . Число таких уравнений меньше числа неизвестных в два раза. Недостающие уравнения получают, приравнивая во внутренних узлах первые и вторые производные, вычисляемые по коэффициентам на соседних участках:  ( x )  bi  2ci ( x  xi 1 )  3d i ( x  xi 1 ) 2 ,  ( x )  2ci  6 d i ( x  xi 1 ), xi 1  x  xi , bi 1  bi  2ci hi  3d i hi2 , (3.4) ci 1  ci  3d i hi , 1  i  N  1. (3.5) Недостающие условия можно получить, полагая, например, что вторая производная равна нулю на концах участка интерполирования:  ( x0 )  0,c1  0, (3.6)  ( x N )  0,c N  3d N hN  0, (3.7) Уравнения (3.2)-(3.7) позволяют определить все 4 N неизвестных коэффициентов: ai ,bi , ci , d i ( 1  i  N ). Решение полученной системы уравнений можно сильно упростить, если привести ее к специальному виду. Используя уравнение (3.2), можно получить все коэффициенты ai . Из (3.5) и (3.7) следует di  c i 1  c i , 3hi dN   1  i  N  1. (3.8) cN 3h N (3.9) Из (3.3) и (3.8): bi  y i  y i 1 c  2c i  hi i 1 , hi 3 Из (3.3) и (3.9): 1 i  N 1 . (3.10) bN  y N  y N 1 2c  hN N hN 3 (3.11) Исключим теперь из (3.4) величины bi и bi 1 с учетом (3.10), наращивая во втором случае индекс на 1, а величину d i - с учетом (3.8). В результате получим систему уравнений для определения коэффициентов c i : c1  0 , hi 1 c i 1  2( hi 1  hi )c i  hi c i 1  3( y i  y i 1 y i 1  y i  2  ) , 2  i  N  1 (3.12) hi hi 1 c N 1  0 . ) После нахождения коэффициентов c i остальные коэффициенты определяют по следующим формулам: a i  y i 1 , 1  i  N , bi  y i  y i 1 c  2c i  hi i 1 , 1 i  N 1 , hi 3 bN  y N  y N 1 2c  hN N , hN 3 di  ci 1  ci , 1 i  N 1 , 3hi dN   cN . 3h N Осталось выяснить, как решать систему (3.12). Матрица этой системы трехдиагональна, т.е. все ее элементы равны нулю, кроме тех, которые находятся на главной и двух соседних диагоналях. Такие системы удобно решать методом прогонки. Суть метода в следующем. Применение метода исключения Гаусса для решения системы уравнений с трехдиагональной матрицей приводит к тому, что система уравнений преобразуется к виду, когда в каждом уравнении содержится только два неизвестных, и при обратном ходе одно из этих неизвестных выражается через другое. Поэтому применительно к (3.12) можно записать: ci   i  1 ci  1   i  1 , где (3.13)  i 1 , i 1 - некоторые не известные пока прогоночные коэффициенты; c i 1   i c i   i . Подставляя последнее выражение в (3.12) и преобразуя, получим ci   hi f i  hi 1 i ci 1  . hi 1 i  2( hi 1  hi ) hi 1 i  2( hi 1  hi ) (3.14) Сравнивая (3.13) и (3.14), получим  i 1   hi f i  hi 1 i , i 1  . hi 1 i  2( hi 1  hi ) hi 1 i  2( hi 1  hi ) (3.15) В этих формулах введено обозначение f i  3( y i  y i 1 y i 1  y i  2  ). hi hi 1 Из условия c1  0 , следует  2  0 , 2  0 . Теперь алгоритм решения (3.12) выглядит следующим образом. По формулам (3.15) при известных  2 , 2 , равных нулю, вычисляют прогоночные коэффициенты  i 1 , i 1 ( 2  i  N ) (прямой ход). Затем по формулам (3.13) при условии c N 1  0 определяют все c i (обратный ход).