Интерполирование

Интерполирование. Интерполирование и приближение функций   7.1. Задача интерполирования и аппроксимации функций Задача интерполирования состоит в том, чтобы по значениям функции f(x) в нескольких  точках отрезка восстановить ее значения в остальных точках данного отрезка. Разумеется, такая постановка задачи допускает сколь угодно много решений. Задача интерполирования возникает, например, в том случае, когда известны  результаты измерений yk = f(xk) некоторой физической величины f(x) в точках xk, k = 0, 1, …, n и требуется определить ее значение в других точках. Интерполирование  используется также при необходимости сгущения таблиц, когда вычисление  значений f(x) по точным формулам трудоемко. Иногда возникает необходимость приближенной замены (аппроксимации) данной  функции (обычно заданной таблицей) другими функциями, которые легче вычислить. При обработке эмпирических (экспериментальных) зависимостей, результаты обычно  представлены в табличном или графическом виде. Задача заключается в аналитическом  представлении искомой функциональной зависимости, то есть в подборе формулы,  корректно описывающей экспериментальные данные.   7.2. Интерполирование алгебраическими многочленами Пусть функциональная зависимость задана таблицей y0 = f(x0);…, y1= f(x1);…,yn  = f(xn). Обычно задача интерполирования формулируется так: найти многочлен P(x) =  Pn(x)степени не выше n, значения которого в точках xi (i = 0, 1 2,…, n) совпадают со  значениями данной функции, то есть P(xi) = yi. Геометрически это означает, что нужно найти алгебраическую кривую вида  (7.1) проходящую через заданную систему точек Мi(xi,yi) (см. рис. 7.1). Многочлен Р(х)  называется интерполяционным многочленом. Точки xi (i = 0, 1, 2, …, n) называютсяузлами интерполяции  Рис. 7.1. Интерполирование алгебраическим многочленом Для любой непрерывной функции f(x) сформулированная задача имеет единственное  решение. Действительно, для отыскания коэффициентов а0, а1, а2 ,…, аn получаем систему линейных уравнений  (7.2) определитель которой (определитель Вандермонда) отличен от нуля, если среди  точек xi (i = 0, 1, 2,…, n) нет совпадающих. Решение системы (7.2) можно записать различным образом. Однако наиболее  употребительна запись интерполяционного многочлена в форме Лагранжа и в форме  Ньютона. Запишем без вывода интерполяционный многочлен Лагранжа:  (7.3) Нетрудно заметить, что старшая степень аргумента х в многочлене Лагранжа равна n.  Кроме этого, несложно показать, что в узловых точках значение интерполяционного  многочлена Лагранжа соответствует заданным значениям f(xi).   7.3. Интерполяционная формула Ньютона Интерполяционная формула Ньютона позволяет выразить интерполяционный  многочлен Pn(x) через значение f(x) в одном из узлов и через разделенные разности  функцииf(x), построенные по узлам x0, x1,…, xn. Эта формула является разностным  аналогом формулы Тейлора:  (7.4) Прежде чем приводить формулу Ньютона, рассмотрим сведения о разделенных разностях. Пусть в узлах  что среди точек xk, k = 0, 1,…, n нет совпадающих. Тогда разделенными разностями  первого порядка называются отношения  известны значения функции f(x).Предполагаем,  Будем рассматривать разделенные разности, составленные по соседним узлам, то есть  выражения  первого порядка можно построить разделенные разности второго порядка: . По этим разделенным разностям   (7.5)  (7.6) Аналогично определяются разности более высокого порядка. То есть пусть известны  разделенные разности k­го порядка  разделенная разность k+1­го порядка определяется как тогда   (7.7) Интерполяционным многочленом Ньютона называется многочлен Показано, что интерполяционный многочлен Лагранжа (7.3) совпадает с  интерполяционным многочленом Ньютона (7.8). Замечания  (7.8)  В формуле (7.8) не предполагалось, что узлы x0, x1,…, xn расположены в каком­то  определенном порядке. Поэтому роль точки x0 в формуле (7.8) может играть любая  из точек x0, x1,…, xn. Соответствующее множество интерполяционных формул  можно получить из (7.8), перенумеровав узлы. Например, тот же самый  многочлен Pn(x) можно представить в виде  (7.9) o Если   то (7.8) называется формулой  интерполирования вперед, а (7.9) ­ формулой интерполирования назад. o Интерполяционную формулу Ньютона удобнее применять в том случае, когда интерполируется одна и та же функция f(x), но число узлов интерполяции  постепенно увеличивается. Если узлы интерполяции фиксированы и  интерполируется не одна, а несколько функций, то удобнее пользоваться  формулой Лагранжа.   7.4. Сходимость интерполяционного процесса Обсудим следующий вопрос: будет ли стремиться к нулю погрешность  интерполирования f(x) – Ln(x), если число узлов n неограниченно увеличивать: 1. Свойства сходимости или расходимости интерполяционного процесса зависят как  от выбора последовательности сеток, так и от гладкости функции f(x). 2. Известны примеры несложных функций, для которых интерполяционный процесс  расходится. Так последовательность интерполяционных многочленов, построенных для непрерывной  функции   по равноотстоящим узлам на отрезке  [­1, 1], не сходится к функции  рис. 7.2 в качестве иллюстрации изображен график многочлена L9(x) при   ни в одной точке отрезка [­1, 1], кроме точек –1, 0, 1. На ,  построенного для функции   по равноотстоящим узлам на отрезке [­1,1]. Рис. 7.2. Сходимость интерполяционных многочленов 3. Чтобы избежать этих некорректностей, в практике вычислений обычно избегают  пользоваться интерполяционными многочленами высокой степени.   7.5. Задача обратного интерполирования Пусть функция y = f(x) задана таблично. Задача обратного интерполирования заключается в том, чтобы по заданному значению функции y определить соответствующее значение  аргумента x. Для случая неравноотстоящих значений аргумента x0, x1,…, xn задача может быть  непосредственно решена с помощью интерполяционного многочлена Лагранжа. В этом  случае достаточно принять переменную y за независимую и написать формулу (аналог  выражения (7.3)), выражающую х как функцию у: . (7.10) Можно также, считая у аргументом, использовать формулу Ньютона: Замечание. Обратное интерполирование корректно только для взаимно однозначных  функций.  (7.11) Пример 7.1. Исходная функция y = f(x) задана табл. 7.1: Таблица 7.1 x y 10 3 15 7 17 11 20 17 Необходимо найти значение функции y при x = 12; найти значение x, для которого y = 10. Решение. В качестве примера задачу прямого интерполирования в начале таблицы с  неравноотстоящими узлами решим по формулам Ньютона (7.8); для обратного  интерполирования применим формулу Лагранжа (7.10). y(12) = f(x0) + (x –x0)f(x0, x1) + (x –x0) (x –x1) f(x0, x1, x2) + + (x –x0) (x –x1) (x –x2) f(x0, x1, x2, x3) = 3 + 2∙0.8 + + 2∙(­3)∙0.02857 + 2∙(­3)∙(­5)∙(­0.002857) = 4.3429. x(10) = 10∙3∙(­1)∙(­7)/[(­4)(­8)(­14)] + 15∙7∙(­1)∙(­7)/[4∙(­4)(­10)] + + 17∙7∙3(­7)/[8∙4∙(­6)] + 20∙7∙4∙1/[14∙10∙6] = 16.641.   7.6. Отыскание параметров эмпирических формул методом наименьших квадратов При эмпирическом (экспериментальном) изучении функциональной зависимости одной  величины у от другой х производят ряд измерений величины у при различных значениях  величины х. Полученные результаты можно представить в виде таблицы, графика: X Y x1 y1 x2 y2 … … xn yn Задача заключается в аналитическом представлении искомой функциональной  зависимости, то есть в подборе функции, описывающей результаты эксперимента. Особенность задачи состоит в том, что наличие случайных ошибок измерений делает  неразумным подбор такой формулы, которая точно описывала бы все опытные значения,  то есть график искомой функции не должен проходить через все экспериментальные  точки. Эмпирическую формулу обычно выбирают из формул определенного типа:  (7.12) Таким образом, задача сводится к определению параметров a, b, c,… формулы, в то время как вид формулы известен заранее из каких­либо теоретических соображений или из  соображения простоты аналитического представления эмпирического материала. Пусть  выбранная эмпирическая зависимость имеет вид  (7.13) с явным указанием всех параметров, подлежащих определению. Эти параметры а0, а1, а2, …, аn нельзя определить точно по эмпирическим значениям функции y0, y1, y2,…, yk, так  как последние содержат случайные ошибки. Таким образом, речь может идти только о получении достаточно хороших оценок  искомых параметров. Метод наименьших квадратов (МНК) позволяет получить  несмещенные и состоятельные оценки всех параметров а0, а1, а2,…, аn.   7.7. Суть метода наименьших квадратов Дальнейшие рассуждения будем проводить в предположении, что все измерения значений  функции y0, y1, y2,…, yn произведены с одинаковой точностью. Тогда оценка  параметров а0, а1, а2,…, аn определяется из условия минимума суммы квадратов  отклонений измеренных значений yk от расчетных f(xk; а0, а1, а2,…, аn):  (7.14) Отыскание же значений параметров а0, а1, а2,…, аn, которые доставляют min значение  функции  (7.15) сводится к решению системы уравнений  (7.16) Наиболее распространен способ выбора функции f(xk; а0, а1, а2,…, аn) в виде линейной  комбинации:  (7.17) Здесь  аn – коэффициенты, определяемые методом наименьших квадратов. Запишем в явном  виде условие (7.16), учитывая выражение (7.17): базисные функции (известные); n << k; а0, а1, а2,…,  Из системы линейных уравнений (7.18) определяются все коэффициенты ak. Система  (7.18) называется системой нормальных уравнений, матрица которой имеет вид  (7.18) Здесь  (7.19)  (7.20) Матрица (7.19) называется матрицей Грама. Расширенную матрицу получим добавлением  справа столбца свободных членов: (7.21), где   (7.22) Основные свойства матрицы Грама 1. Матрица симметрична относительно главной диагонали, то есть  . 2. Матрица является положительно определенной. Следовательно, при решении  методом Гаусса можно воспользоваться схемой единственного деления. 3. Определитель матрицы будет отличен от нуля, если в качестве базиса  выбраны линейно независимые функции  ; в этом случае система (7.18) имеет единственное решение. В качестве базисных можно выбрать линейно независимые степенные функции Следует учесть, что n << k. Тогда для этих функций расширенная матрица Грама примет  вид  (7.23) Если выбрать n = k, то на основании единственности интерполяционного полинома  получим функцию  , совпадающую с каноническим интерполяционным полиномом   (7.24) степени k. При этом аппроксимирующая кривая пройдет через все экспериментальные  точки, и функция S будет равна нулю. Пример 7.2. Исходная функция y = f(x) задана в виде табл. 7.2: Таблица 7.2 x y 10 3 15 7 17 11 20 17 Аппроксимируем экспериментальные данные линейной либо квадратичной функцией.  Методом наименьших квадратов необходимо уточнить коэффициенты  аппроксимирующего полинома. Решение 1. При линейной аппроксимации исходную зависимость представим в  виде  определим a0 и a1. Расширенная матрица Грама в нашем случае имеет вид , где  . Методом наименьших квадратов  =>  ; а1 = 1.3774; а0= ­11.8491. Таким образом, аппроксимирующая функция равна Оценим погрешность формулы, и результаты этой оценки сведем в табл. 7.3: Таблица 7.3 x 10 15 17 y 3 7 11 f 1.9249 8.8119 y ­ f |y­f| / |y| 1.0751 0.3584 ­1.8119 0.2588 11.5667 ­0.5667 0.0515 20 17 15.6989 1.3011 0.07654 Для нашей линейной функции S1 = 6.4528. 2. Решим ту же задачу, аппроксимировав эмпирические данные полиномом второй  степени:  . Матрица Грама в этом случае имеет вид Все результаты сведены в табл. 7.4. Таблица 7.4 x 10 15 17 20 y 3 7 11 17 f 2.9511 7.3381 y ­ f |y­f| / |y| 0.0489 0.0163 ­0.3381 0.0483 10.6007 0.3993 0.0363 17.1101 ­0.1101 0.0065 S2 = 0.2883. Обсуждение результатов 1. Аппроксимировав эмпирические результаты более простой функцией  (линейной), мы получили погрешность в различных узловых точках, лежащую  в пределах от 5 до 35 %. 2. Более сложная формула квадратичной интерполяции обеспечивает  погрешность не более 5 %. 3. Косвенную оценку погрешности можно провести, сравнив значения S1 и S2. 4. Матрица Грама для полинома второй, третьей степени имеет простой вид и  может быть решена, например, методом Гаусса. Вопросы для самопроверки  В чем заключается задача интерполирования и аппроксимации?  Запишите интерполяционные формулы Лагранжа и Ньютона.  Какие требования предъявляются а) к интерполяционным полиномам;  б) к аппроксимационным полиномам?  Что такое разделенные разности?  В каких случаях применяются формулы Ньютона для интерполирования  а) вперед, б) назад?  Что можно сказать о сходимости интерполяционных полиномов?  Что такое обратное интерполирование, при каких условиях оно возможно  (корректно)?  В чем заключается идея метода наименьших квадратов?  Что такое матрица Грама, каковы ее свойства?  Что такое базисные функции? Можно ли в качестве базисных функций выбрать а)  линейно независимые функции; б) линейно зависимые функции? Среднеквадратичное приближение функций  Алгоритмы* На днях нужно было написать программу, вычисляющую среднеквадратичное  приближение функции, заданной таблично, по степенному базису — методом наименьших квадратов. Сразу оговорюсь, что тригонометрический базис я не рассматривал и в этой  статье его брать не буду. В конце статьи можно найти исходник программы на C#. Теория Пусть значения приближаемой функции f(x) заданы в N+1 узлах f(x0), ..., f(xN).  Аппроксимирующую функцию будем выбирать из некоторого параметрического  семейства F(x, c), где c = (c0, ..., cn)T — вектор параметров, N > n. Принципиальным отличием задачи среднеквадратичного приближения от задачи  интерполяции является то, что число узлов превышает число параметров. В данном  случае практически всегда не найдется такого вектора параметров, для которого значения аппроксимирующей функции совпадали бы со значениями аппроксимируемой функции во всех узлах. В этом случае задача аппроксимации ставится как задача поиска такого вектора  параметров c = (c0, ..., cn)T, при котором значения аппроксимирующей функции как можно  меньше отклонялись бы от значений аппроксимируемой функции  F(x, c) в совокупности  всех узлов. Графически задачу можно представить так Запишем критерий среднеквадратичного приближения для метода наименьших квадратов: J( c) = √ (Σi=0 N[f(xi) — F(x, c) ]2)  min→ Подкоренное выражение представляет собой квадратичную функцию относительно коэффициентов аппроксимирующего многочлена. Она непрерывна и дифференцируема  по c0, ..., cn. Очевидно, что ее минимум находится в точке, где все частные производные  равны нулю. Приравнивая к нулю частные производные, получим систему линейных  алгебраических уравнений относительно неизвестных (искомых) коэффициентов  многочлена наилучшего приближения. Метод наименьших квадратов может быть применен для различных параметрических  функций, но часто в инженерной практике в качестве аппроксимирующей функции  используются многочлены по какому­либо линейно независимому базису {φk(x), k=0,...,n}: F(x, c)  = Σk=0 n[ckφk(x)]. В этом случае система линейных алгебраических уравнений для определения  коэффициентов будет иметь вполне определенный вид: a00c0 + a01c1 +… + a0ncn = b0 a10c0 + a11c1 +… + a1ncn = b1 … an0c0 + an1c1 +… + anncn = bn akj =  Σi=0 N    [φk(xi)φj(xi) ], bj =  Σi=0 N[f(xi)φj(xi) ] Чтобы эта система имела единственное решение необходимо и достаточно, чтобы  определитель матрицы А (определитель Грама) был отличен от нуля. Для того, чтобы  система имела единственное решение необходимо и достаточно чтобы система базисных  функций φk(x), k=0,...,n была линейно независимой на множестве узлов аппроксимации. В этой статье рассматривается среднеквадратичное приближение многочленами по  степенному базису  {φk(x) = xk, k=0,...,n}. Пример А теперь перейдем к примеру. Требуется  вывести эмпирическую формулу для  приведенной табличной зависимости  f(х),используя метод наименьших квадратов. x y 0,75 2,50 1,50 1,20 2,25 1,12 3,00 2,25 3,75 4,28 Примем в качестве аппроксимирующей функцию y = F(x) = c0 + c1x + c2x2, то есть, n=2, N=4 Система уравнений для определения коэффициентов: a00c0 + a01c1 +… + a0ncn = b0 a10c0 + a11c1 +… + a1ncn = b1 … an0c0 + an1c1 +… + anncn = bn akj =  Σi=0 N[φk(xi)φj(xi) ], bj =  Σi=0 N[f(xi)φj(xi) ] Коэффициенты вычисляются по формулам: a00 = N + 1 = 5, a01 = Σi=0 a10 = Σi=0 a20 = Σi=0 b0 = Σi=0 Nxi = 11,25,  a11 = Σi=0 2 = 30,94,  a12 = Σi=0 Nxi Nxi 3 = 94,92,  a22 = Σi=0 2 = 30,94,  a21 = Σi=0 Nxi Nxi Nxi Nyi = 11,25,  b1 = Σi=0 Nxi = 11,25,  a02 = Σi=0 Nxiyi = 29,  b2 = Σi=0 Nxi 2yi = 90,21 Nxi 2 = 30,94 3 = 94,92 4 = 303,76 Решаем систему уравнений и получаем такие значения коэффициентов: c0 = 4,822, c1 = ­3,882, c2 = 0,999 Таким образом y = 4,8 — 3,9x + x2 График получившейся функции Релизация на C# А теперь перейдем к тому, как написать код, который бы строил такую матрицу. А тут,  оказывается, все совсем просто: private double[,] MakeSystem(double[,] xyTable, int basis) {   double[,] matrix = new double[basis, basis + 1];   for (int i = 0; i < basis; i++)   {     for (int j = 0; j < basis; j++)     {       matrix[i, j] = 0;     }   }   for (int i = 0; i < basis; i++)   { for (int j = 0; j < basis; j++)     {       double sumA = 0, sumB = 0;       for (int k = 0; k < xyTable.Length / 2; k++)       {         sumA += Math.Pow(xyTable[0, k], i) * Math.Pow(xyTable[0, k], j);         sumB += xyTable[1, k] * Math.Pow(xyTable[0, k], i);       }       matrix[i, j] = sumA;       matrix[i, basis] = sumB;     }   }   return matrix; } На входе функция получает таблицу значений функций — матрицу, в первом столбце  которой содержатся значения x, во втором, соответственно, y, а также значение  степенного базиса. Сначала выделяется память под матрицу, в которую будут записаны коэффициенты для  решения системы линейных уравнений. Затем, собственно, составляем матрицу — в sumA записываются значения коэффициентов aij, в sumB — bi, все по формуле, указанной выше в теоретической части. Для решения составленной системы линейных алгебраических уравнений в моей  программе используется метод Гаусса. Архив с проектом можно скачать по ссылке. Скриншот работы программы на примере, решенном выше: Используемые источники: Сулимова В.В. Методические указания по курсу «Вычислительный практикум» —  Тула, ТулГУ, 2009 — 65 с. Среднеквадратичное приближение 1. Наилучшее приближение. Интерполяция позволяет легко аппроксимировать функцию  . Однако точность такой  аппроксимации гарантирована лишь в небольшом интервале порядка нескольких шагов  сетки. Для другого интервала приходится заново вычислять коэффициенты  интерполяционной формулы. Нам же всегда желательно иметь единую приближенную  формулу  заданную и аппроксимирующую функции на большом отрезке.  пригодную длябольшого отрезка  . Поэтому далее будем сравнивать определены неточно  При интерполяции мы приравниваем значения у  — например, из эксперимента, ­ то точное приравнивание неразумно. Поэтому нередко   в узлах. Если  целесообразней приближать функцию не по точкам, а в среднем, т. е. в норме  Пусть заданы функция  принадлежащие линейному  и множество функций      нормированному пространствуфункций. Нас интересуют две      задачи. Первая — аппроксимация с заданной точностью: по заданному   найти такую    чтобы выполнялось неравенство   Второе — нахождение наилучшего приближения, т. е. функции   удовлетворяющей соотношению Существует ли наилучшее приближение и единственно ли оно (для данных функции и  множества)? Это имеет место не при любом выборе пространства и множества. Например, в пространстве   выберем функцию   и множество   тогда В самом деле, при  е. двум.  эта норма равна площади заштрихованной трапеции на рис. 10, а, т.   эта норма, согласно рис. 10, б, равна площади заштрихованной трапеции (которая При  опять равна двум) плюс площади заштрихованных треугольников. Значит, для любого с,  по модулю меньшего единицы,  приближение здесь существует, но оно не единственно. минимизирует норму отклонения, т. е. наилучшее Рис. 10. Выведем достаточное условие существования наилучшего приближения. Пусть  в линейном     пространствефункций выбрано множество, образованное функциями вида   где функции  есть линейное    можно считать линейно­независимыми. Это множество      подпространство нашего пространства. Изменим один из коэффициентов  суммы (37) на величину   из неравенства треугольника (1.3) следует т. е. норма   непрерывно зависит от  . Очевидно,   также есть непрерывная  функция коэффициентов  . Рассмотрим нормы как функции координат  . Сфера есть замкнутое ограниченное множество, поэтому   на этой сфере имеет точную  нижнюю грань   и в силу непрерывности достигает ее при некотором   Очевидно,    в противном случае  , что противоречит линейной независимости  Возьмем шар   где   — какое­то положительное число. В силу однородности нормы ­  функции вне этого, шара   и, следовательно,   Значит, вне этого шара норма  погрешности заведомо далека от нижней грани. Только внутри шара   достаточно  близки по норме. Но шар — ограниченное и замкнутое     множествозначений координат      поэтому непрерывная функция координат   достигает на нем точной нижней грани. Следовательно, в любом линейном аппроксимации (37) наилучшее приближение существует, хотя не во  всяком линейном     пространстве оно единственно.     нормированном пространстве при линейной На практике используются пространства   и С. В этом параграфе рассмотрим  приближения в пространстве  , т. е. среднеквадратичную аппроксимацию.