Аппроксимация функций

Тема 3. АППРОКСИМАЦИЯ ФУНКЦИЙ 1. Интерполяция Формулы Ньютона. Дано: функция у(х), задана таблицей своих значений в n + 1 точках (узлах): х0 х1 …. xn у0 у1 …. yn Требуется определить (хотя бы приближенно) ее значение у(х) при х  хj, т.е. требуется построить аппроксимирующую функцию Y(x)  y(x). При дополнительных условиях Y(xj) = y(xj) (полное совпадение в узловых точках) и Y(x) = Pn (x) (многочлен n-й степени) имеем классическую задачу многочленной интерполяции. Для ее решения существует много методов, наиболее удобна интерполяционная формула Ньютона: P n (x) = y0 + (х – x0) 1[x0, x1] + (х – x0) (х – x1) 2[x0, x1, x2] + … + (х – x0)(х – x1) … (x– xn–1)) n[x0, x1,… xn], где величины 1[x0, x1], 2[x0, x1, x2], …, n[x0, x1,… xn] называются разделенными разностями, определяемыми по рекуррентным формулам: , j = 1, 2, … n; , j = 2, 3, … n. ….. , j = k+1, k+2, … n. Нижний индекс при  – порядок разности. Пример 3.1. Дана таблица значений некоторой функции при n = 4. Построим многочлен Ньютона. Строим таблицу разделенных разностей (табл. 3.1). Разности, стоящие в выделенных ячейках используются как коэффициенты формулы Ньютона: P n (x) = 0 + (х – 0)  0.693147 + (х – 0) (х – 1)  (-0.143841) + (х – 0)(х – 1) (x– 2)  0.028317 + (х – 0)(х – 1)(х – 2)(х – 3)  ( -0.004861). Вычислим Pn (1.5) = 0.918. Точное значение интерполируемой функции у = ln(1 + x) равно у(1.5) = 0.916, т.е. точность достаточно большая. Таблица 3.1 j 0[xj] = yj 1[x0, xj] 2[x0, x1, xj] 3[x0, x1, x2, xj] 4[x0, x1, x2, x3, x4] - - - - 1 0.693147 0.693147 - - - 2 1.098612 0.549306 -0.143841 - - 3 1.386294 0.462098 -0.115525 0.028317 - 4 1.609438 0.402359 -0.096929 0.023456 -0.004861 Иногда ставится обратная задача: требуется по известному у определить соответствующее значение х – задача обратной интерполяции. Для ее решения применяется обратная интерполяционная формула Ньютона: Qn (y) = x0 + (y – y0) 1[x0, x1] + (y – y0) (y – y1) 2[x0, x1, x2] + … + (y – y0)(y – y1) … (y– yn–1)) n[x0, x1,… xn], где величины 1[x0, x1], 2[x0, x1, x2], …, n[x0, x1,… xn] называются обратными разделенными разностями, определяемыми по рекуррентным формулам: , j = 1, 2, … n; , j = 2, 3, … n. ….. , j = k+1, k+2, … n. Свойства многочленов Ньютона 1) Если все узлы различны, то многочлен Pn (х) всегда существует и единственен. Для наличия аналогичных свойств у Qn (у) дополнительно требуется монотонность функции у(х). 2) Построенный многочлен не зависит от порядка чередования узлов в таблице – это вытекает из свойства 1) (и это очень удобно при подключении новых узлов для построения полинома более высокого порядка). Замечание. В практике вычислений полиномы степени выше пятой обычно не используют, т.е. число узлов интерполяции не должно превышать шести. Если при таком числе узлов не обеспечивается заданная точность, следует уменьшать расстояние между узлами. Программы. Приведем фрагменты программ интерполяции. 1. Разделенные разности. Введем массив D[0..n, 0..n], элементы которого являются разделенными разностями, первый индекс – номер узла, второй индекс – порядок разности (напомним, что 0[xi] = уi). for i := 0 to n do D[i,0] := y[i]; for j := 1 to n do for i := j to n do D[i,j] := ( D[i,j–1] – D[j–1,j–1] ) / ( x[i] – x[j–1] ); 2. Многочлен. Для вычисления значения Pn (z) удобно использовать обобщенную схему Горнера (обычная схема Горнера – см. “Классификация численных методов. Прямые методы”, пример 1.16). read(z) ; p := D[n, n] ; for i := n – 1 downto 0 do p := p*( z – x[i] ) + D[i,i] ; write (p) ; 2. Приложение интерполяции к вычислительным задачам 2.1. Численное решение уравнений. Пусть на отрезке [a, b] надо найти корень уравнения f(x) = 0. Если функция f(x) монотонна на [a, b], то одним из эффективных методов решения является обратная интерполяция. Предлагается следующий алгоритм. 1. Выбираем некоторый шаг h и просчитываем значения у = f(х) с этим шагом на [a, b]. 2. Задаемся некоторым значением n и выбираем на [a, b] n + 1 точку. Для повышения точности рекомендуется выбирать точки так, чтобы примерно половина из них была с одной стороны от корня (в них f(xi) < 0), а другая половина была с другой стороны от корня (в них f(xi) > 0). Тогда корень будет близок к середине отрезка интерполяции, т.е. там, где точность максимальна. 3. Строим Qn (y) и вычисляем х* = Qn (0). Это и будет приближенным корнем уравнения. 4. Для повышения точности можно применить итерационный вариант метода. Для этого задаемся некоторым  > 0 и проверяем, достигнута ли заданная точность. В качестве критерия можно использовать неравенства | xn – x0 | <  или | f(x*) | < . Первое (оценка по аргументу) означает, что интервал интерполяции, на котором находится корень, достаточно мал. Второе (оценка по функции) – что в найденной точке функция приблизительно равна 0. Если выбранный критерий точности удовлетворен, то стоп. 5. Если точность не достигнута, то к имеющимся узлам интерполяции добавляем новый – (х*, f(x*)). Можно применять одну из двух стратегий модификации набора узлов: а) стратегия добавления – не менять старые узлы; при этом общее число узлов и степень многочлена возрастает на 1. б) стратегия замены – исключить один из старых узлов, например, тот, в котором значение |уi| максимально (т.е. наиболее отлично от 0); общее число узлов и степень многочлена не изменяется. После этого вернуться к пункту 3. Пример 3.2. Найти корень уравнения f(x)  cos(x) – x = 0 на интервале [–1, 2.5] с точностью по функции 10–5. Положим h = 0.5 и составим таблицу значений у(х) (таблица 3.2): Таблица 3.2. x – 1.0 – 0.5 0.5 1.0 1.5 2.0 2.5 y 1.54 1.38 1.0 0.38 – 0.46 – 1.43 – 2.42 – 3.30 В качестве узлов интерполяции будем использовать точки при х = 0; 0.5; 1.0; 1.5. Строим матрицу D для обратных разделенных разностей: 0 0 0 0 0.5 -0.8033194 0 0 D = 1 -0.6850734 -0.1412264 0 1.5 -0.6174713 -0.1028578 -0.039573 Получаем многочлен Q3 (y) = 0 – 0.8033194 (у – 1) – 0.1412264 (у – 1) (у – 0.377583) – – 0.039573 (у – 1) (у – 0.377583) (у + 0.459698). Отсюда Q3 (0) = х* = 0.74312596; f(x*) = – 0.006769 – точность не достигнута. Заменяем новой точкой (0.74312596; – 0.006769) старый узел (1.5; –1.429263), в котором значение функции наиболее отличается от 0. Получаем новую матрицу D 0 0 0 0 0.5 -0.8033194 0 0 D = 1 -0.6850734 -0.1412264 0 0.743126 -0.7381297 -0.1696096 -0.0626661 и новый многочлен Q3 (y) = 0 – 0.8033194 (у – 1) – 0.1412264 (у – 1) (у – 0.377583) – – 0.0626661 (у – 1) (у – 0.377583) (у + 0.459698) Q3 (0) = х* = 0.7391176; f(x*) = – 0.000054 – точность не достигнута. Заменяем на новую точку (0.7391176; – 0.000054) старый узел (0; 1). Получаем новую матрицу D 0.7391176 0 0 0 0.5 -0.6331945 0 0 D = 1 -0.5675757 -0.0783713 0 0.743126 -0.5969738 -0.0942385 -0.0350323 Q3 (y) = 0.7391176 – 0.6331945 (у – 0.000054) – 0.0783713 (у – 0.000054) (у – 0.377583) – – 0.0350323 (у – 0.000054) (у – 0.377583) (у + 0.459698) Q3 (0) = х* = 0.73908513; f(x*) = 1.488  10 – 9 – точность достигнута. 2.2. Поиск минимума функции одной переменной. Пусть даны три точки на плоскости: A = (x0, y0); B = (x1, y1); С = (x2, y2); x0 < x1 < x2. Def. Три точки на плоскости (x0, y0), (x1, y1), (x2, y2) образуют выпуклую тройку, если при x0 < x1 < x2 выполнены также неравенства у1 < у0, у1 < у2. Ставится задача: найти приближенно точку минимума функции у(х) по трем точкам, образующим выпуклую тройку. Эта задача решается с помощью параболической интерполяции. Согласно этому методу строят многочлен Ньютона 2-й степени Р2 (х) (параболу) и находят точку его минимума, т.е. точку, в которой . Имеем: P2 (x) = y0 + (х – x0) 1[x0, x1] + (х – x0) (х – x1) 2[x0, x1, x2]. Найдем производную, приравняем ее к нулю и решим полученное линейное уравнение. Получим: . (3.1) Свойство. Если узлы интерполяции образуют выпуклую тройку, то хmin [x0, x2], т.е. точка минимума, получаемая методом параболической интерполяции, устойчива. Если минимизируемая функция у(х) задана аналитически, то решение можно получить с любой заданной точностью. Для этого используют итерационный вариант метода параболической интерполяции. Cтроят новую точку F = (хmin , у(хmin )) и проверяют, достигнута ли требуемая точность. Если да, то хmin – решение задачи, СТОП. В противном случае заменяют точкой F один из старых узлов, например тот, в котором значение уi максимально, и повторяют интерполяцию. В качестве критерия точности по функции можно использовать условие max { уi } – min { уi } < , т.е. значения функции в узлах совпадают с заданной точностью. В качестве критерия точности по аргументу можно использовать условие max { хi } – min { хi } < , т.е. значения аргумента лежат в достаточно малом интервале. Перед применением метода выбирают некоторый шаг h и просчитывают значения у(х) с этим шагом на [a, b] для обнаружения выпуклой тройки. Замечание 1. Так как на последних итерациях значения функции очень близки, то вычисления надо проводить с большим числом значащих цифр. Замечание 2. Если требуется найти максимум функции, то для того, чтобы не менять все выкладки и формулу (3.1), удобнее поменять знак функции и искать min (– y(x)). Пример 3.3. Найти минимум функции f(x) = 0.5 x2 – sin(x) на отрезке [–1, 2.5] с точностью по аргументу 10–3. Положим h = 0.5 и составим таблицу значений у = f(х) (таблица 3.3). График функции – на рис. 3.1. Таблица 3.3. x – 1.0 – 0.5 0.5 1.0 1.5 2.0 2.5 y 1.34 0.60 0.0 –0.35 – 0.34 0.13 1.09 2.53 В качестве узлов интерполяции будем использовать точки при х = 0; 0.5; 1.0, т.к. они образуют выпуклую тройку. Строим матрицу D для разделенных разностей: 0 0 0 D = -0.35443 -0.70885 0 -0.34147 -0.34147 0.73476 Отсюда F = (хmin , у(хmin )) = (0.732369; –0.4004509). После 1-й итерации точность не достигнута, следовательно, заменяем точкой F узел при х = 0, т.к. у(0) = 0 = max. Получаем новую матрицу D: -0.40045 0 0 D = -0.35443 -0.19807 0 -0.34147 0.22038 0.8369 (хmin , у(хmin )) = (0.7345207095; -0.4004712). Точность не достигнута, следовательно, заменяем точкой F узел при х = 1, т.к. у(1) = -0.341470985 = max, и т.д. После пятой итерации получаем следующую таблицу узлов: Таблица 3.4. х 0.7390872156 0.7393051349 0.7390850585 у -0.400488612 -0.400488572 -0.400488612 (хmin , у(хmin )) = (0.7390851332; -0.400488572); max { хi } – min { хi } = 2.20110–4 – точность достигнута. 3. Среднеквадратичное приближение. Метод наименьших квадратов Пусть, как и в п.1, функция у(х), задана таблицей своих значений в N точках (для удобства мы поменяли нумерацию). x1 х2 …. xN y1 у2 …. yN Требуется найти аппроксимирующую функцию вида Y(x) = b1 f1 (x) + b2 f2 (x) + … + bk fk (x), (3.2) где fj(x) – известные функции (базисные функции); bj – неизвестные коэффициенты; k < N. В качестве меры близости между у(х) и Y(x) задано условие , т.е. в отличие от задачи интерполяции, точного совпадения в узловых точках не требуется. Такая задача называется задачей среднеквадратичного приближения. Искомые коэффициенты определяются методом наименьших квадратов (МНК). Приведем матричную запись этого метода, как наиболее компактную. Обозначим: Ф = [Фij] = [fj(xi)] – матрица размером N  k (структурная или регрессионная матрица); b – вектор искомых коэффициентов; у – вектор значений функции у(х) в узлах хi. Тогда выражение b1f1(xi) + … + bkfk(xi) = Фb представляет собой N-мерный вектор. Для определения вектора b коэффициентов имеем систему линейных уравнений: (ФTФ) b = ФTy, (3.3) которую легко решить, например, методом Гаусса. В качестве меры точности полученного решения часто используют величину sСРЕД = – среднеквадратичная погрешность аппроксимации (b* – найденное решение системы (3.3)). Пример 3.4. Для данных примера 3.1 методом наименьших квадратов построить аппроксимирующий многочлен 2-й степени. Имеем f1 (x)  1; f2 (x) = x; f2 (x) = x2. Следовательно, матрица Ф имеет вид: Отсюда Решив систему (3.3), получим Иными словами, Y(х) = 0.0239 + 0.6937 х – 0.0756 х2. На рис. 3.3 изображены графики исходной и аппроксимирующей функций. 4. МНК при нелинейных моделях Предположим, что уравнение, связывающее вход и выход исследуемого объекта, нелинейно по параметрам, т.е. имеет вид: у = F(x, b) (3.4) Принцип наименьших квадратов легко обобщается и на нелинейный случай. Решением (вектором коэффициентов) задачи среднеквадратичной аппроксимации при нелинейной модели (3.4) будем понимать то значение вектора bЕk, для которого сумма квадратов отклонений (b) =  (F(xi, b) – yi)2 принимает минимальное значение. В отличии от линейного случая, в котором минимум находится сравнительно просто – как решение системы линейных уравнений, поиск минимума (b) в общем случае достаточно сложен. В нелинейной модели может существовать несколько локальных минимумов функции (b). По определению МНК-решение соответствует глобальному минимуму. Линеаризация. В некоторых случаях уравнение (3.4) можно преобразовать к линейному виду: g(у) = g(F(x, b)) = b1f1(x) + … + bkfkx). (3.5) Введя новую переменную z = g(y), получаем обычный МНК. Такая процедура называется линеаризацией. К ней прибегают, когда непосредственная минимизация (b) по каким-либо причинам затруднена. Решить линейную систему намного проще, чем минимизировать сложную функцию. Пример 3.5. 1) . Положим g(y) = 1/y, получим ax + b = 1/y. Отсюда получим линейный случай:  ( axi + b – 1/yi ) 2  min. 2) . Положим g(y) = ln(y), получим ax + b = ln(y). Отсюда получим линейный случай:  ( axi + b – ln(yi) ) 2  min. При наличии ошибок в исходных данных преобразование (3.5) обычно приводит к потере точности, т.к. минимизация проводится уже преобразованной суммы квадратов, чье оптимальное значение не соответствует минимуму (b). Пример 3.6. Истинная зависимость y от х имеет вид: . Имеем данные эксперимента при наличии случайных ошибок : xi = 0.8, 1,…, 3; у = (1.496, 1.039, 0.57, 0.5141, 0.5428, 0.3272, 0.2415, 0.2065, 0.374, 0.1348, 0.091, 0.2631) (см. рис. 3.4). В результате линеаризации получили модель . Результаты аппроксимации приведены на рис. 3.5 и 3.6. Из рис. 3.6 видно, что при начальных значениях х погрешность аппроксимации особенно велика. Взвешенный МНК. Расхождение между истинным значением функции и ее аппроксимацией, полученной в результате линеаризации, можно значительно уменьшить с помощью введения соответствующих весов i для экспериментальных точек, т.е. минимизировать сумму: 1(b) = i2 (g(F(xi, b)) – g(yi)) 2 . Опишем соответствующую процедуру, предложенную С. А. Айвазяном. Предположим, что экспериментальные значение у представляют собой сумму у = F(x, b) + . Разложим функцию g(y) по формуле Тейлора до линейных членов  в окрестности F(x, b), положив g(y) = g(F(x, b) + )  g(F(x, b)) + g’y  = b1f1(x) + … + bkfk(х) + g’y . Отсюда . Так как в МНК минимизируют сумму i 2, то получаем 1 (b) =  [b1f1(xi) + … + bkfk(хi) – g(yi) ]2 , т.е. . В матричной записи МНК это означает умножение каждой строки линеаризованной матрицы Ф и вектора g(y) на i, или умножение матрицы Ф и g(y) слева на диагональную матрицу R, в которой Ri i = i. Результат применения процедуры Айвазяна к примеру 3.6 изображен на рис. 3.7 (сравните с рис. 3.6). Получена аппроксимация .