Использование рациональной интерполяции часто целесообразнее интерполяции полиномами в случае функций с резкими изменениями характера поведения или особенностями производных в точках.
5.6.9. Метод наименьших квадратов (МНК)
При обработке экспериментальных данных, полученных с некоторой погрешностью, интерполяция становиться неразумной. В этом случае целесообразно строить приближающую функцию таким образом, чтобы сгладить влияние погрешности измерения и числа точек эксперимента. Такое сглаживание реализуется при построении приближающей функции по методу наименьших квадратов.
Рассмотрим совокупность значений таблично заданной функции fi в узлах хi при i=0,1,…,n. Предположим, что приближающаяся функция F(x) в точках х1, х2, …, хn имеет значения . Будем рассматривать совокупность значений функции f(x) и функции F(x) как координаты двух точек n-мерного пространства. С учетом этого задача приближения функции может быть определена другим образом: найти такую функцию F(x) заданного вида, чтобы расстояние между точками M(f1, f2, …, fn) и было наименьшим. Воспользовавшись метрикой евклидова пространства, приходим к требованию, чтобы величина
была наименьшей, что соответствует следующему:
(5.12)
то есть сумма квадратов должна быть наименьшей. Задачу приближения функции f(х) теперь можно формулировать иначе. Для функция f(х), заданной таблично, необходимо найти функцию F(x) определенного вида так, чтобы сумма квадратов (5.12) была наименьшей.
Выбор класса приближающихся функций определяется характером поведения точечного графика функции f. Это могут быть линейная зависимость, любые элементарные функции и т.д.
Практически вид приближающей функции F можно определить, построив точечный график функции f(х), а затем построить плавную кривую, по возможности наилучшим образом отражающую характер расположения точек. По полученной кривой выбирают вид приближающей функции.
Когда вид приближающей функции выбран, то последующая задача сводится к отысканию значений параметров функции. Рассмотрим метод нахождения параметров приближающей функции в общем виде на примере приближающей функции с тремя параметрами f=F(x, с, b, с). Тогда имеем
(5.13)
Сумма квадратов разностей соответствующих значений функций f и F будет иметь вид:
(5.14)
Сумма (5.14) является функцией φ(a, b, с) трех переменных a, b, с. Задача сводится к отысканию ее минимума. Для этого используем необходимое условие экстремума:
или
(5.15)
Решив эту систему (5.15) трех уравнений с тремя неизвестными относительно параметров a, b, с, получим конкретный вид искомой функции F(x, a, b, с). Изменение количества параметров не приведет к изменению сущности метода, а отразится только на количестве уравнений в системе (5.15).
Как следует из начальных условий, найденные значения функции F(x, а, b, с) в точках x1, х2, …, хn будут отличаться от табличных значений y1, у2, …, уn. Значение разностей
fi-F(xi, a, b, c,) = εi, i=1, 2, …, n
будет определять отклонение измеренных значений f от вычисленных по формуле (5.14). Для найденной эмпирической формулы (5.14) в соответствии с исходными табличными данными можно найти сумму квадратов отклонений
(5.16)
Она, в соответствии с принципом наименьших квадратов для заданного вида приближающей функции и ее найденных параметров (параметры a, b, с), должна быть наименьшей. Из двух разных приближений одной и той же табличной функции, следуя принципу наименьших квадратов, лучшим нужно считать тот, для которого сумма (5.16) имеет меньшее значение.
5.6.10. Тригонометрическая интерполяция рядами Фурье
При тригонометрической интерполяции используются тригонометрические полиномы — линейные комбинации тригонометрических функций sin(nx) и cos(nx). Этот вид интерполирования применяется для процессов, которые отражают циклические процессы, связанные с периодическими функциями [52–54]. Известно, что такие функции удобно представлять в виде тригонометрического ряда или его частичной суммы с достаточной степенью точности.
Функциональный ряд вида
(5.17)
называется тригонометрическим. Его коэффициенты аn и bn — действительные числа, не зависящие от х. Если этот ряд сходится для любого х из промежутка [-π, π], тогда он определяет периодическую функцию f(x) с периодом Т=2π. Ряд вида (5.17) называется рядом Фурье для интегрируемой на отрезке [-π, π] функции f(х), если коэффициенты его вычисляются по следующим правилам:
(5.18)
(5.19)
(5.20)
В практических расчетах, как правило, ограничиваются конечным числом первых членов ряда Фурье. В результате получается приближенное аналитическое выражение для функции f(х) в виде тригонометрического полинома N-го порядка
Но соотношения для вычисления коэффициентов Фурье (5.18)–(5.20) пригодны для случая аналитического задания исходной функции. Если функция задана в виде таблицы, то возникает задача приближенного отыскания коэффициентов Фурье по конечному числу имеющихся значений функции.
Таким образом, формулируется следующая задача практического, гармонического анализа: аппроксимировать на интервале (0, T) тригонометрический полином N-го порядка функцию у=f(х), для которой известны m ее значений уk=f(хk) при хk=kТ/m, где k=0, 1, 2, …, m-1.
Тригонометрический полином для функции, определенной на интервале (0, Т), имеет вид:
(5.21)
Коэффициенты аn и bn определяются следующими соотношениями:
(5.22)
(5.23)
Применяя в соотношениях (5.22)–(5.23) формулу прямоугольников для вычисления интегралов по значениям подынтегральных выражений в точках хk=kT/m, где k=0, 1, 2, …, m-1, имеем
(5.24)
(5.25)
Таким образом, тригонометрический полином (5.21), коэффициенты аn и bn находятся по формулам (5.24)–(5.25), служит решением поставленной задачи. При этом, коэффициенты (2.44)–(2.45) минимизируют сумму квадратов отклонений
В случае, когда m=2N коэффициенты аn и bn для n=0, 1, 2, …, N определяется соотношениями (5.24)–(5.25), а коэффициент aN определяется соотношением:
Сам же полином QN(x) становится интерполяционным полиномом, так как в этом случае при любом bN выполняется соотношения QN(xk)=yk для всех хk=kТ/m, где k = 0, 1, 2, …, m-1.
5.7. Аппроксимация зависимостей в Maple
5.7.1. Аппроксимация аналитически заданных функций
В Maple 9.5 если функция задана аналитически, то наиболее простым способом нахождения ее аппроксимирующей зависимости является применение функции convert, которая позволяет представить функцию в виде иного выражения, чем исходное. Например, при опции polynom осуществляется полиномиальная аппроксимация. Это поясняют следующие примеры (файл aprox):
> convert(taylor(exp(х),х,5),polynom);
> f:=х->(х^3+х)/(х^2-1);
> convert(f(x),parfrac, x);
На рис. 5.9 представлен пример полиномиальной аппроксимации хорошо известной статистической функции erfc(x). Для полинома задана максимальная степень 12, но ввиду отсутствия в разложении четных степеней максимальная степень результата оказывается равна 11.
Рис. 5.9. Пример полиномиальной аппроксимации функции erfc(x) — сверху построены графики исходной функции и полинома, снизу график абсолютной погрешности
Как видно из приведенного рисунка, в интервале изменения x от -1,4 до 1,4 аппроксимирующее выражение почти повторяет исходную зависимость. Однако затем график аппроксимирующей функции быстро отходит от графика исходной зависимости и погрешность аппроксимации резко возрастает. При этом он ведет себя иначе даже качественно, никоим образом не показывая асимптотическое поведение, характерное для исходной зависимости. Это говорит о том, что полиномиальная аппроксимация плохо подходит для экстраполяции (предсказания) зависимостей.
