Интерполяционная формула ньютона mathcad

БлогNot. Полином Ньютона для неравноотстоящих узлов в MathCAD

Полином Ньютона для неравноотстоящих узлов в MathCAD

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

Смысл такого подхода к расчёту состоит в том, что при добавлении к данным нового узла формула многочлена удлиняется только на 1 слагаемое, так что не приходится все расчёты выполнять заново. Впрочем, если мы реализуем данный расчёт "в лоб", и не вручную, а с помощью MathCAD, то пусть себе считает заново 🙂

Итак, заданы векторы X и Y , содержащие исходные данные — координаты точек _i,Y_i> . Сначала нужно построить матрицу разделённых разностей, её размерность будет на 1 меньше, чем количество точек:

После этого можно написать подпрограмму для вычисления значения полинома Ньютона в точке x малое (MathCAD, как и C++, соблюдает регистр символов в именах переменных).

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

Осталось нарисовать картинку, в приложенном документе она есть.

Скачать пример в формате .XMCD — полином Ньютона для неравноотстоящих узлов (102 Кб)

04.07.2013, 17:40; рейтинг: 21757

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

В этом случае шаг таблицы h = х_i+ 1 — x_i (i = 0 , 1 , 2 , . n) = const является величиной постоянной. Для таких таблиц построение интерполяционных формул (как впрочем, и вычисление по этим формулам) заметно упрощается.

Прежде чем перейти к рассмотрению этого вопроса, познакомимся с понятием конечных разностей.

Пусть функция задана таблицей с постоянным шагом. Разности между значениями функции в соседних узлах интерполяции называют конечными разностями первого порядка:

Из конечных разностей первого порядка образуются конечные разности второго порядка:

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

D 2 y

D 3 y

D y ₀

D 2 y ₀

D 3 y ₀

D y ₁

D 2 y ₁

D 3 y ₁

D y ₂

D 2 y ₂

D y ₃

Конечные разности любого порядка могут быть представлены через значения функции. Действительно, для разностей первого порядка это следует из определения. Для разностей второго порядка имеем:

Для разностей 3-го порядка

Методом математической индукции можно доказать, что:

D k y_i = y_{i + k} — k y_{i +k —} 1 + y _{i + k — 2} — . + (— 1 ) k y_i.

Пример 4. Построить конечные разности для функции

D 2 P (x) = [3 (х + 1) 2 + 3 (х + 1) + 1] -(3 х 2 + 3 х + 1) = 6 х + 6,

Обратите внимание, что конечная разность третьего порядка функции P (x) постоянна.

Вообще, справедливо утверждение: если

P_n (x) =

Пример 5. Составить таблицу конечных разностей функции

у = 2 х 3 — 2 х 2 + 3 х — 1

от начального значения x₀ = 0, приняв шаг h = 1.

Полагая x ₀ = 0 , x ₁ = 1 , x ₂ = 2, находим соответствующие значения y ₀ = -1 , y ₁ = 2 , y ₂ = 13. Отсюда имеем:

Эти значения заносим в таблицу (Таблица 2).

D 2 y

D 3 y

Так как наша функция есть полином третьей степени, то третья разность ее постоянна и равна D 3 y_i = 2 ? 3! = 12.

Поэтому дальнейшее заполнение Таблицы 2 можно производить при помощи суммирования, используя формулы

Первая интерполяционная формула Ньютона

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

+ a ₃ (x — x ₀ ) (x — x ₁ ) (x — x ₂ ) + ... + a_n (x — x ₀ ) (x — x ₁ ) ? ... ? (x — x_n- 1 ).

Это полином степени n. Значение коэффициентов a ₀ , a ₁ , . a_n Найдем из условия совпадения значений исходной функции и полинома P_n(x) a узлах интерполяции. Полагая х = х ₀ , из (13) находим:

Далее, придавая х значение х ₁ и х ₂ последовательно, получаем:

a _{2 = = ! =!.}

Aналогично можно найти и другие коэффициенты. Общая формула имеет вид:

а_к =!, k = 0 , 1 , 2 . n

Подставляя эти выражения в формулу (13), получим следующий вид интерполяционного полинома:

… + !(x — x ₀ ) (x — x ₁ ). (x — x_n- 1 ).

Практически формула (14) применяется в несколько ином виде. Положим

t = ,

т. е. х = х _{0 + ht. Тогда}

Полученное выражение называется первой интерполяционной формулой Ньютона для интерполирования вперед.

Интерполяционную формулу (15) обычно используют для вычисления значений функций в точках левой половины рассматриваемого отрезка. Это объясняется следующим. Разности D k у_i вычисляются через значение функции у_i, у_{i +} 1 , . у_{i + k}, причем i + k ? n, поэтому при больших значениях i мы не можем вычислить разности высших порядков (k ? n — i). Например, при i = n — 3 в (15) можно учесть только D у, D 2 у, D 3 у.

Если в формуле (15) положить n = 1 , то получим формулу линейной интерполяции:

При n = 2 будем иметь формулу квадратичной интерполяции:

P₂ (x) = у ₀ + t D у ₀ + .

Если дана неограниченная таблица значений у, то число n в интерполяционной формуле (15) может быть любым. Практически в этом случае число n выбирают так, чтобы разность D n у_i была постоянной с заданной степенью точности. За начальное значение x₀ можно принять любое табличное значение аргумента х.

Если таблица значений функции конечна, то число n ограничено, а именно: n не может быть больше числа значений функции у, уменьшенного на единицу.

Пример 6. Приняв шаг h = 0,05, построить на отрезке [3,5; 3,6] интерполяционный полином Ньютона для функции у = е х , заданной таблицей

33,115

34,813

36,598

38,475

40,447

Составляем таблицу разностей (Таблица 3).

D 2 y

D 3 y

33,115

34,813

36,598

38,475

40,447

Заметим, что в столбцах разностей, следуя обычной практике, мы не указываем десятичные разряды (которые ясны из столбца значений функции). Так как разности третьего порядка практически постоянны, то в формуле (15) полагаем n = 3. Приняв x₀ = 3,50, y₀ = 33,115, будем иметь:

P₃ (x) = 33,115 + 1,698 t + 0,087+ 0,005

P₃ (x) = 33,115 + 1,698 t + 0,0435+ 0,00083,

.

На практике часто необходимо для функции, заданной таблично, подобрать аналитическую формулу, представляющую с некоторой точностью данные табличные значения функции. Такая формула называется эмпирической, причем задача построения ее неоднозначна.

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

Пример 7. Построить эмпирическую формулу для функции у, заданной таблично

D 2 y

Составляя таблицу разностей (Таблица 4), убеждаемся, что вторая разность постоянна.

Используя интерполяционную формулу Ньютона в форме (14) и учитывая, что h = 1, будем иметь:

у = 5,2 + 2,8 х —

Вторая интерполяционная формула Ньютона

Для правой половины рассматриваемого отрезка разности лучше вычислять справа налево. В этом случае:

t= ,

o. е. t 0 и интерполяционную формулу Ньютона можно получить в виде:

Формулу (16) называют второй интерполяционной формулой Ньютона для интерполирования назад.

Пример 8. Дана таблица значений y = lg x семизначных логарифмов

3,0000000

3,0043214

3,0086002

3,0128372

3,0170333

3,0211893

Составляем таблицу разностей (Таблица 5).

Таблица разностей функции у = lg x

D 2 y

D 3 y

3,0000000

43214

3,0043214

42788

3,0086002

42370

3,0128372

41961

3,0170333

41560

3,0211893

Используя подчеркнутые разности, в силу формулы (16) будем иметь:

lg 1044 = 3,0211893 + (-0,6) * 0,0041560 + * 0,0000401 + * 0,0000008 = 3,0187005.

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

Таким образом, первая интерполяционная формула Ньютона обычно используется для интерполирования вперед и экстраполирования назад, а вторая интерполяционная формула Ньютона, наоборот — для интерполирования назад и экстраполирования вперед.

Операция экстраполирования менее точна, чем операция интерполяции в узком смысле слова.

График интерполяционного полинома у = j (х) проходит через заданные точки, т. е., значения полинома j (х) и данной функции у = ¦ (х) совпадают в узлах х = х_i (i = 0 , 1 , . . n). Если функция ¦ (х) сама является полиномом степени n, то имеет место тождественное равенство ¦ (х) = j (х). В общем случае в точках, отличных от узлов интерполяции

Эта разность и есть погрешность интерполяции, и называется остаточным членом интерполяционной формулы. Оценим его значение.

Предположим, что заданные числа y_i являются значениями некоторой функции y = ¦ (х) в точках х = х_i . Пусть эта функция непрерывна и имеет непрерывные производные до n + 1 порядка включительно. Можно показать, что в этом случае остаточный член интерполяционного полинома Лагранжа имеет вид:

R_L (x) = ¦ (n+ 1 ) ( x )

где ¦ (n+ 1 ) ( x ) — производная (n+ 1 )-го порядка функции ¦ (х) в некоторой точке х = x , x _I [x ₀ , x_n]. Если максимальное значение этой производной равно:

M_{n + 1} = | ¦ (n+ 1 ) (x) | ,

то можно записать формулу для оценки остаточного члена интерполяционного полинома Лагранжа

| R_L (x) | R_P (x) = ¦ (n+ 1 ) ( x ) h n + 1 , . (18)

Остаточный член второй интерполяционной формулы Ньютона:

R_P (x) = ¦ (n+ 1 ) ( x ) h n + 1 , . (19)

Если предположить, что разности D n + 1 у_n почти постоянны для функции y = f(x) и h достаточно мало, и учитывая, что

¦ (n+ 1 ) (х) = ,

приближенно можно положить:

¦ (n+ 1 ) ( x ) » .

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

R_P (x) » D n+ 1 у _{0 .}

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

R_P (x) » D n+ 1 у_n .

Пример 9. Оценим ошибку интерполяции при отыскании lg 1044 в Примере 8.

Составленная таблица разностей (Таблица 5) показывает, что третьи разности функции можно считать практически постоянными, следовательно n = 3. По формуле (19) имеем:

R₃(x) = ¦ IV ( x ) h 4 ,

Так как f(x) = lg₁₀ х, то ¦ IV (х) = — , поэтому

| ¦ IV ( x ) | » 10 -10 .

Таким образом, интерполирование с учетом третьих разностей, т. е. с применением интерполяционного полинома третьей степени, не вносит погрешностей до десятого знака.

Следует подчеркнуть, что существует один и только один интерполяционный полином при заданном наборе узлов интерполяции. Формулы Лагранжа, Ньютона и др. порождают один и тот же полином (при условии, что вычисления проводятся точно).

Повышение точности интерполяции целесообразно проводить за счет уменьшения шага и специального расположения точек х_i. Выбор способа интерполяции определяется различными соображениями: точностью, временем вычисления, погрешностями округления и др.

Исправляем ошибки: Нашли опечатку? Выделите ее мышкой и нажмите Ctrl+Enter

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

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

Составим разности значений функции:

Эти значения называются первыми разностями (или разностями первого порядка) функции.

Можно составить вторые разности функции:

Аналогично составляются разности порядка

Конечные разности можно выразить непосредственно через значения функции. Например,

Аналогично для любого можно записать

Эту формулу можно записать и для значения разности в узле

Используя конечные разности, можно определить

Перейдем к построению интерполяционного многочлена Ньютона. Этот многочлен будем искать в следующем виде:

График многочлена должен проходить через заданные узлы, т.е. . Эти условия используем для нахождения коэффициентов многочлена:

Найдем отсюда коэффициенты

Аналогично можно найти и другие коэффициенты. Общая формула имеет вид

Подставляя эти выражения в формулу (9), получаем следующий вид интерполяционного многочлена Ньютона:

Конечные разности могут быть вычислены по формуле (7).

Формулу (10) часто записывают в другом виде. Для этого вводится переменная тогда

С учетом этих соотношений формулу (10) можно переписать в виде

Полученное выражение может аппроксимировать заданную функцию на всем отрезке изменения аргумента . Однако более целесообразно (с точки зрения повышения точности расчетов и уменьшения числа членов в (11) ограничиться случаем , т.е. использовать формулу (11) для Для других значений аргумента, например для вместо лучше взять значение . Таким образом, интерполяционный многочлен Ньютона можно записать в виде

Полученное значение называется первым интерполяционным многочленом Ньютона для интерполирования вперед.

Интерполяционную формулу (12) обычно используют для вычисления значений функции в точках левой половины рассматриваемого отрезка. Это объясняется следующим. Разности вычисляются через значения функции , , причем поэтому при больших значениях мы можем вычислить разности высших порядков . Например, при в (12) можно учесть только .

Для правой половины рассматриваемого отрезка разности лучше вычислять справа налево. В этом случае

т.е. , и интерполяционный многочлен Ньютона можно получить в виде

Полученная формула называется вторым интерполяционным многочленом Ньютона для интерполирования назад.

Рассмотрим пример использования многочленов Ньютона для интерполирования.

Пример 1. Функция y = f (x) задана таблицей значений с шагом h.

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

«>