Приближение функций

Рассмотрим задачи, в которых требуется найти способ приближенного вычисления значения функции в заданной точке. Эти задачи возникают в тех случаях, когда прямое вычисление данной функции у = f(x)для произвольного аргумента невозможно или слишком трудоемко.

Например, функция f(x)определяется как решение сложной задачи. Здесь могут быть известны некоторые свойства функции у = f(x),например, непрерывность и дифференцируемость.

Даже если функция легко вычисляется, может воз­никнуть необходимость ее замены, например, при вычис­лении некоторых определенных интегралов или специ­альных функций математической физики (ниже будут приведены примеры). В этом случае вместо подынтег­ральной функции нужно подобрать другую функцию, от которой интеграл легко вычисляется. Разумеется, новая функция должна быть приближенно равна в некотором смысле подынтегральной функции.

Если значения функции определяются в результате дорогостоящих экспериментов, могут быть найдены ее значения только в некоторых точках, а для вычисления значения в произвольной точке требуется приближенный метод. При этом может быть известен вид функции, но неизвестны параметры, входящие в определение функ­ции. В этом случае задача сводится к определению пара­метров известной функции.

Аппроксимацией(приближением) функции f(x)назы­вается нахождение такой функции g(x)(аппроксимирую­щей функции),которая была бы близка заданной. Крите­рии близости функций могут быть различны.

В том случае, когда приближение строится на дискрет­ном наборе точек (x_i, y_i), i = 1, 2, ... , n, аппроксимацию называют точечнойили дискретной.

При точечной квадратичной аппроксимациипара­метры а₁, а₂, ... , а_т аппроксимирующей функции g = g(x, а₁, а₂, ... , а_т), т £ п,определяются из условия

.

Если аппроксимация проводится на непрерывном мно­жестве точек (отрезке), она называется непрерывнойили интегральной.

При интегральной квадратичной аппроксимациифункции у = f(х) на отрезке [а, b] параметры аппрокси­мирующей функции g(x, а₁, а₂,... , а_m) определяются из условия

.

Примером непрерывной аппроксимации может слу­жить использование конечного числа слагаемых разло­жения функции в ряд Тейлора, т. е. замена функции мно­гочленом.

Наиболее часто встречающимся видом точечной аппрок­симации на дискретном наборе из (n+1)-й точки (x_i, y_i), i = 0, 1, ... , п является интерполяциямногочленом n-го порядка Р_п(х),коэффициенты которого определяются из условий

y_i = P_n(x_i), i = 0, 1, ... , п.

Применяя интерполяционный многочлен, можно вы­числить значения функции f(x)между узлами (провести интерполяцию в узком смысле),а также определить зна­чение функции за пределами заданного интервала (про­вести экстраполяцию).Следует иметь в виду, что по­грешность экстраполяции может быть велика.

В том случае, когда интерполяционный многочленедин для всей области интерполяции, говорят, что интер­поляция глобальная.Если между различными узлами интерполяционные многочлены различны, говорят о ку­сочнойили локальной интерполяции.Простейшим слу­чаем локальной интерполяции является кусочно-линей­ная интерполяция, когда в качестве интерполяционной функции выбирается полином первой степени, т.е. узловые точки соединяются отрезками прямой.

Интерполяция

Под задачей интерполирования,или интерполяциипринято понимать следующую:

Пусть известны значения y_i функции f(x)в точках х_i i = 0, 1,..., п,называемых узлами интерполяции.Требу­ется найти такую функцию F(x)(интерполирующая фун­кция),значения которой в узлах интерполяции совпада­ют со значениями f(x):

F(x_i) = y_i, i = 0,l, ... , n. (4.1)

Предполагается, что интерполирующая функция F(x)зависит от конечного числа параметров, которые находят из условий (4.1). Если зависимость интерполирующей функции F(x)от параметров линейна, то говорят о линей­ной интерполяции,в противном случае – о нелинейной интерполяции.В последнем случае нахождение парамет­ров из системы (4.1) может быть трудной задачей.

Формула у = F(x)называется интерполяционнойи используется для вычисления значения функции f(x)в точке х,не совпадающей с узлами интерполяции. Эта операция называется интерполированием(интерполяци­ей) функции.При этом, если точка х не принадлежит от­резку [a, b], a = min(x₀, х₁,... , х_п), b = mах(x₀, х₁,... , х_п),то говорят об экстраполировании функции.

Часто в качестве интерполирующей функции исполь­зуется многочлен такого порядка Р_п(х)(интерполяционный многочлен),удовлетворяющий условию

y_i = P_n(x_i), i = 0, l,..., п. (4.2)

Так как многочлен n-го порядка определяется своими коэффициентами, то число параметров равно n +1 и ус­ловия (4.2) представляют систему линейных уравнений.

(4.3)

Определителем системы (4.3) является определитель Вандермонда

(4.4)

который отличен от нуля, если среди узлов интерполяции нет совпадающих. Отсюда следует, что задача построения интерполяционного многочлена разрешима и существует единственный интерполяционный многочлен, удовлетво­ряющий условиям (4.2). Как увидим ниже, формы запи­си интерполяционного многочлена могут быть разными. Отметим, что задача отыскания произвольного много­члена,удовлетворяющего условиям (4.2), некорректна. Если степень многочлена меньше п,то такого многочле­на не существует, но если больше п,то таких многочле­нов бесконечное множество.

4.1.Интерполяционные формулы Ньютона

Первой конечной разностьюфункции у = f(x)называ­ется выражение

Dу = f(x + Dх) – f(x). (4.5)

Конечная разность второго порядкаопределяется формулой

D²у = D(Dy)= D[f(х + Dх) - f(х)] = f(x + 2Dх) - 2 f(х + Dх) + f(x). (4.6)

Используя формулу бинома Ньютона, можно вывести формулу конечной разности п-го порядка:

Dⁿу = D(D^n-1y)= f(х + nDх) - f[х + (n - 1)Dx] + f[x + (n – 2)Dх] - ...+ (-1)ⁿ f(x), где . (4.7)

Справедливо утвержде­ние: конечная разность(n+1)-го порядка многочлена п-го порядка равна нулю.

Пусть функция f(x)задана своими значениями y_i = f(x_i)в равноотстоящих точках x_i = х₀ + ih, i = 0, 1,..., п.

Таблица конечных разностейфункции f(х) записыва­ется в одной из двух форм: горизонтальной или диаго­нальной таблицы разностей. Так как конечная разность использует 2 значения, столбец Dy_i содержит п значений, D²y_i – на одно значение меньше и т. д. Если задано п +1 значений функции, то таблица конечных разностей со­держит п столбцов, причем последний столбец содержит только одно значение. Общий вид горизонтальной таб­лицы, конечных разностейприведен в табл. 1.

Таблица 1

Обобщенной степенью п-го порядкачисла х называется произведение

x^[n] = x(x - h)(x -2h)(x - 3h)...[x - (n - 1)h], (4.8)

где h = const, а для n = 0 полагают x⁽⁰⁾ = 1.
Найдем конечные разности для обобщенной степени:

Dх^[n] = (х + h)⁽ⁿ⁾ - x⁽ⁿ⁾ = (x + h)x(x - h)(x -2h)...[x - (n - 2)h] - x(x - h)(x -2h)...[x - (n - 1)h] = x(x - h)(x -2h)(x - 3h)...[x - (n - 2)h]{( x + h) - [x - (n - 1)h]} = nh x^{[n - 1]}

Мы получили следующую формулу:

Dх^[n] = nh x^{[n - 1]} (4.9)

Пользуясь методом математической индукции, можно вывести и общую формулу

D^kх^[n] = n(n -1)...[n -(k - 1)] h ^kx^{[n - k]}. (4.10)

Выведем первую интерполяционную формулу Ньютона. Найдем многочлен Р_п(х), удовлетворяющий условиям Р_п(х_i)= y_i для i = 0, 1,..., n. Будем искать многочлен Р_п(х)в следующем виде:

P_п(х)= а₀ + a₁(х – х₀) + a₂(x – x₀)(x – x₁) + ... + a_n(x – x₀)(x – x₁)...(x – x_n) =

= а₀ + a₁(х – х₀)^[1] + a₂(x – x₀)^[2] + ... + a_n(x – x₀)^[n].

Из условия Р_п(х₀)= y₀ следует a₀ = y₀.

Найдем первую конечную разность многочлена Р_п(х):

DР_n(х) = a₁h + 2hа₂(x – x₀)^[1] + ... + na_n(x – x₀)^[n-1].

Отсюда при х = х₀получим DР_n(х₀) = a₁h = Dy₀, т.е. a₁ = Dy₀/h. Аналогично можно получить общую формулу коэффициентов многочлена Ньютона:

a_k = . (4.11)

Теперь можно записать первую интерполяционную формулу Ньютона:

P_п(х)= y₀ + (х – х₀)^[1] + (x – x₀)^[2] + ... + (x – x₀)^[n] (4.12)

С помощью замены переменной q = (х - x₀)/h первую интерполяционную формулу Ньютона можно предста­вить в виде

P_п(х)= y₀ + q + + ... + . (4.13)

Аналогичными рассуждениями выводится вторая интерполяционная формула Ньютона:

P_п(х)= y_n + (х - х_n) + (x - x_n)(x - x_n-1) +...+ (x - x_n)(x - x_n-1)...(x - x₁), (4.14)

которая с помощью замены q = (х – x_n)/h приводится к виду

P_п(х)= y_n + q + + ... + . (4.15)

Первая интерполяционная формула Ньютона приме­няется для интерполирования вблизи начала таблицы (около х₀), а вторая – для интерполирования вблизи кон­ца таблицы (около х_n).

4.2.Интерполяционная формула Лагранжа

Интерполяционные формулы Ньютона (4.12), (4.15) пригодны лишь для равноотстоящих узлов интерполирования. Рассмотрим интерполяционную формулу Лагран­жа, которая применима при любом расположении узлов. Пусть дана система точек х₀, х_1, ..., х_n, принадлежащих некоторому отрезку [а, b], и известны соответствующие значения функции y_i = f(x_i).Найдем многочлен L_n(x),удовлетворяющий условиям L_n(х_i) = y_i.

При построении многочлена Лагранжа используются многочлены n-й степени р_i(х), принимающие значение 1 в точке x_i и нулевые значения в остальных узлах интер­поляции x_j, j ¹ i.Так как x_j,при j ¹ i являются корнями многочлена р_i(х), то справедливо разложение р_i(х) на мно­жители

р_i(х) = C_i(х - х₀)(х - х₁)...(х - х_i-1)(х - х_i+1) ...(х - x_n), i = 0, 1, ... , n.

Из условия p_i(x_i) = 1 находим значение константы C_i:

C_i= 1/(х_i - х₀)(х_i - х₁)...(х_i - х_i-1)(х_i - х_i+1) ...(х_i - x_n) и получаем выражение для многочленов p_i(x_i):

р_i(х) = . (4.16)

Интерполяционный многочлен Лагранжаимеет вид

L_n(x) = = × (4.17)

Формулу Лагранжа можно преобразовать так, чтобы упростить вычисления. Вынесем за знак суммы произве­дение

П_n+1(x) = (х - х₀)(х – х₁) ... (х - х_n),

которое является многочленом степени п +1. Получим

L_n(x) = П_n+1(x) .

Теперь формулу Лагранжа можно записать в виде

L_n(x) = П_n+1(x) , (4.18)

где D_i = (х_i - х₀)(х_i - х₁)...(х_i - х_i-1)(x - x_i)(х_i - х_i+1) ...(х_i - x_n).

Обратите внимание на то, что в произведении D_i из n + 1 сомножителя вместо скобки (х_i - x_i)присутствует множитель (х – x_i).

4.3.Остаточный член интерполяционной формулы Лагранжа

Определение 4.1. Многочлен Р_п(х) равномерно прибли­жает на отрезке [а, b] функцию f(x) с точностью до ε, если выполняется неравенство

| f(x) – P_n(x)| £ ε. (4.21)

Приведем без доказательства теорему Вейерштрасса.

Теорема 4.1.Если функция f(x)непрерывна на отрез­ке [а, b],то для любого ε > 0 найдется многочлен Р_n(х) достаточно высокой степени n, который равномерно при­ближает на отрезке [а, b]функцию f(x)с точностью до ε, т. е. выполняется (4.21).

Определение 4.2. Многочлен Р_п(х)называется много­членом наилучшего приближениядля функции f(x)на отрезке [а, b],если для любого многочлена Q_n(x) степени n выполняется неравенство

| f(x) – P_n(x)| £ | f(x) – Q_n(x) |. (4.22)

Многочлены Чебышеваопределяются рекуррентными формулами

Т₀(х) = 1, Т₁(х) = х, T_n+1(x) = 2хТ_п(х)- Т_п-1(х)при n ³ 1. (4.23)

Выпишем многочлены Чебышева Т_п(х) для п = 2, 3, 4, 5:

Т₂(х)= 2хТ₁(х) – T₀ = 2х² – 1,

Т₃(х) = 2хТ₂(х)– T₁ = 4х³ – 3х,

Т₄(х) = 8х⁴ – 8х² + 1,

Т₅(х) = 16х⁵ – 20х³ + 5х.

Старший коэффициент многочлена Т_n(х), т. е. коэффи­циент при х^п,равен 2^n-1.

Справедливо представление многочленов Чебышева через тригонометрические функции

T_n(x)= cos(n×arccos х), при п ³ 0. (4.24)

Поэтому |Т_n(х)| £ 1 для | х | £ 1. (4.25)

Корни многочлена Чебышева принадлежат отрезку [-1; 1]:

x_k = , k = 0, 1, ... , n - 1. (4.26)

Многочлены Чебышева Т_п(х) с четными индексами являются четными функциями, а с нечетными – нечет­ными функциями.

Многочлены Чебышева Т_n(х) обладают следующим замечательным свойством. Если умножить Т_п(х)на 2^1-n, получится многочлен, наименее уклоняющийся от нуляна отрезке [-1; 1].

Теорема 4.2. Если Р_п(х)– произвольный многочлен степени п со старшим коэффициентом 1, справедливо не­равенство

| P_n(x) | ³ | 2^1-n×T_n(x) | = = 2^1-n , (4.27)

где = 2^1-n×T_n(x).

В некоторых учебниках по численным методам многочленом Чебышева называется многочлен

Т_n(х) = 2^1-n×T_n(x), наименее уклоняющийся от нуляна отрезке [-1; 1].

Поясним смысл теоремы 4.2 на примерах.

Если п = 2, то теорема 4.2 утверждает, что наибольшее значение любой квадратной функции вида х² + рх + q на отрезке [-1; 1] не меньше 2^1-2 = 0,5.

Наибольшее значение любого многочлена вида х³ + рх² + qx + r на отрезке [-1; 1] не меньше 2^1-3 = 0,25.

Среди всех квадратных функций вида х² + рх + q наи­менее уклоняющимся от нуля на отрезке [-1; 1] много­членом является функция

2^1-2×T₂(x) = х² - 0,5.

Среди всех многочленов вида х³ + рх² + qx + r наиме­нее уклоняющимся от нуля на отрезке [-1; 1] является многочлен

2^1-3Т₃(х)= (4х³ - 3х)/4 = х³ - 0,75х.

Для произвольного отрезка [а, b]многочлен со стар­шим коэффициентом 1, наименее уклоняющийся от нуля, получается из заменой

.

Этот многочлен имеет вид

= (b - a)ⁿ×2^1-2n×T_n

Корнями многочлена являются точки

x_k = , k = 0, 1, ... , n - 1. (4.29)

Многочлены Чебышева используются для минимиза­ции погрешности интерполяционной формулы за счет оптимального выбора узлов интерполяции. В формуле остаточного члена (4.19) у многочлена П_n+1(х)старший коэффициент равен 1. Для минимизации погрешности интерполяционной формулы для функции f(x) на отрезке [а, b]нужно взять в качестве узлов интерполяции точки

x_k = , k = 0, 1, ... , n. (4.30)

являющиеся корнями многочлена . Тогда по­грешность интерполяции оценивается неравенством

| f(x) – L_n(x)|£ . (4.31)

4.5. Интерполяция сплайнами

Сплайн– гладкая кривая, проходящая через задан­ные точки (х_i, y_i), i = 0, 1, ... , п. Интерполяция сплай­нами заключается в том, что на каждом отрезке [х_i-1, x_i]используется многочлен определенной степени. Наиболее часто применяется многочлен третьей степени, реже – второй или четвертой. При этом для определения коэффициентов многочленов используются условия непре­рывности производных в узлах интерполяции.

Интерполяция кубическими сплайнамипредставляет собой локальную интерполяцию, когда на каждом отрез­ке [х_i-1, x_i], i = 1, 2, ... , п применяется кубическая кри­вая, удовлетворяющая некоторым условиям гладкости, а именно, непрерывности самой функции и ее первой и вто­рой производных в узловых точках. Использование куби­ческой функции вызвано следующими соображениями. Если предположить, что интерполяционная кривая соот­ветствует упругой линейке, закрепленной в точках (х_i, y_i),то из курса сопротивления материалов известно, что эта кривая определяется как решение дифференциального уравнения f^(IV)(x) = 0 на отрезке [х_i-1, x_i](для простоты из­ложения мы не рассматриваем вопросы, связанные с физи­ческими размерностями). Общим решением такого уравне­ния является многочлен 3-й степени с произвольными коэффициентами, который удобно записать в виде
S_i(x) = а_i + b_i(х - x_i-1) + с_i(x - x_i-1)²+ d_i(x - x_i-1)³,
х_i-1 £ х £ х_i, i = 1, 2, ... , п.(4.32)

Коэффициенты функций S_i(x)определяются из усло­вий непрерывности функции и ее первой и второй произ­водных во внутренних узлах x_i, i = 1, 2,..., п - 1.

Из формул (4.32) при х = х_i-1 получим

S_i(x_i-1) = y_i-1 = a_i, i = 1, 2,..., п,(4.33)

а при х = х_i

S_i(x_i) = а_i + b_ih_i + с_ih_i²+ d_ih_i³,(4.34)

i = 1, 2,..., n.

Условия непрерывности интерполяционной функции записываются в виде S_i(x_i) = S_i-1(x_i), i = 1, 2, ... , n - 1 и из условий (4.33) и (4.34) следует, что они выполнимы.

Найдем производные функции S_i(x):

S'_i(x) = b_i + 2с_i(х - x_i-1) + 3di(х – x_i-1)²,

S"_i(x) = 2c_i + 6d_i(x - x_i-1).

При x = x_i-1, имеем S'_i(x_i-1) = b_i, S" (x_i-1) = 2с_i, а при х = х_i получим

S'_i(x_i) = b_i + 2с_ih_i + 3dih_i², S" (x_i) = 2с_i + 6d_ih_i.

Условия непрерывности производных приводят к уравнениям

S'_i(x_i) = S'_i+1(x_i) Þ b_i + 2с_ih_i + 3dih_i² = b_i+1,

i = l, 2,... , п - 1. (4.35)

S"_i (x_i) = S"_i+1(x_i) Þ 2с_i + 6d_ih_i = 2c_i+1,

i = l, 2,..., n - 1. (4.36)

Всего имеем 4n – 2 уравнений для определения 4n не­известных. Чтобы получить еще два уравнения, исполь­зуют дополнительные краевые условия, например, требо­вание нулевой кривизны интерполяционной кривой в концевых точках, т. е. равенства нулю второй производ­ной на концах отрезка [а, b] а = х₀, b = х_n:

S"₁(x₀) = 2c₁ = 0 Þ с₁ = 0,

S"_n(x_n) = 2с_n + 6d_nh_n = 0 Þ с_n + 3d_nh_n = 0. (4.37)

Систему уравнений (4.33)–(4.37) можно упростить и получить рекуррентные формулы для вычисления коэф­фициентов сплайна.

Из условия (4.33) имеем явные формулы для вычисле­ния коэффициентов a_i:

a_i = y_i-1, i= 1,..., n. (4.38)

Выразим d_i через c_i с помощью (4.36), (4.37):

; i = 1, 2,...,n; .

Положим с_n+1= 0, тогда для d_i получим одну формулу:

, i = 1, 2,...,n. (4.39)

Подставим выражения для а_i и d_i в равенство (4.34):

, i = 1, 2,..., n.

и выразим b_i, через с_i:

, i = 1, 2,..., n. (4.40)

Исключим из уравнений (4.35) коэффициенты b_i и d_i с помощью (4.39) и (4.40):

,

i = 1, 2,..., n -1.

Отсюда получим систему уравнений для определения с_i:

(4.41)

Система уравнений (4.41) может быть переписана в виде

(4.42)

Здесь введено обозначение

, i =1, 2,..., n - 1.

Решим систему уравнений (4.42) методом прогонки. Из первого уравнения выразим с₂через с₃:

c₂ = a₂c₃ + b₂, , . (4.43)

Подставим (4.43) во второе уравнение (4.42):

h₂(a₂c₃ + b₂) + 2(h₂ + h₃)c₃ + h₃c₄ = g₂,

и выразим с₃ через с₄:

с₃ = a₃с₄ + b₃, (4.44)

Предполагая, что с_i-1= a_i-1c_i + b_i-1 из i-го уравне­ния (4.42) получим

c_i = a_iс_i+1 + b_i

, i = 3,..., n – 1, a_n = 0, (4.45)

, i = 3,..., n.

Сформулируем алгоритм интерполяции с помощью кубического сплайна.

Исходные данные: значения функции у₀, у₁,..., у_п в узлах х₀, х₁,..., х_п (х₀ < х₁< ... < х_n).

0. Вычислим значения h_i и g_i:

h_i = x_i – x_i-1 i = 1, 2,..., n,

, i =1, 2,..., n - 1. (4.46)

1. Прямой ход прогонки для вычисления коэффици­ентов c_i. Вычислим коэффициенты прогонки a_i и b_i:

, ,

, i = 3,..., n – 1, a_n = 0, (4.47)

, i = 3,..., n.

2. Обратный ход прогонки для вычисления коэффици­ентов с_i. Вычислим коэффициенты c_i:

c_n+1 = 0,

c_i = a_iс_i+1 + b_i, i = n, n -1,..., 2, (4.48)

c₁ = 0.

3. Вычисление коэффициентов а_i, b_i, d_i:

a_i = y_i-1,

, (4.49)

i = 1, 2,..., n.

4. Вычисление значения функции с помощью сплай­на. Для этого найти такое значение i,что данное значе­ние переменной х принадлежит отрезку [x_i-1, x_i] и вы­числить

S_i(x) = а_i + b_i(х - x_i-1) + с_i(x - x_i-1)²+ d_i(x - x_i-1)³. (4.50)

4.6. Аппроксимация. Метод наименьших квадратов