Интерполяционный полином Ньютона

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

P n(x)= A 0+ A 1(x - x 0)+ A 2(x - x 0)(x - x 1)+...+ A n(x - x 0)(x - x 1)...(x - x n-1)

(5.6)

Коэффициенты этого полинома A ₀, A ₁, A ₂,..., A _n определяются из условий Лагранжа (5.3).

Полагаем x = x ₀. Тогда в (5.6) все слагаемые, кроме A ₀, обращаются в нуль, следова­тель­но,

A ₀ = f ₀.

Затем полагаем x = x ₁, тогда из (5.3) имеем:

f ₀ + A ₁(x ₁- x ₀)= f ₁,

откуда находим коэффициент A ₁:

A ₁ = или A ₁ = f ₀₁.

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

При x = x ₂ полином (5.6) принимает вид:

P _n(x)= f ₀+ f ₀₁(x - x ₀)+ A ₂(x - x ₀)(x - x ₁),

откуда с учетом (5.3) получаем:

f ₂ = f ₀+ f ₀₁(x ₂- x ₀)+ A ₂(x ₂- x ₀)(x ₂- x ₁) или f ₂ - f ₀- f ₀₁(x ₂- x ₀) = A ₂(x ₂- x ₀)(x ₂- x ₁),

следовательно, коэффициент A ₂:

A ₂= = = = ,

где . Обозначая = f ₀₁₂ (разделенная разность второго порядка),

окончательно получаем выражение для A ₂:

A ₂ = f ₀₁₂ .

Аналогично, при x = x ₃, находим коэффициент A ₃:

,

где ; .

Методом математической индукции можно получить для любого A _k (k=0,...,n) следующее выражение:

.

Полученные результаты сведены в представленной ниже таблице.

Следует отметить, что добавление новых узлов в исходных данных не изменяет уже вычисленные коэффициенты; таблица будет лишь дополняться новыми строками и столбцами.

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

 

		Разделенные разности
x	f
x 0	f 0	= A 0
x 1	f 1	f 01=	= A 1
x 2	f 2	f 02=	f 012=	= A 2
x 3	f 3	f 03=	f 013=	f 0123=	= A 3
x 4	f 4	f 04=	f 014=	f 0124=	f 01234=	= A 4

 

Рис.5.4. Алгоритм интерполяции полиномом Ньютона

не нужно организовывать двумерный массив. Более того, все разделенные разности, в том числе и диагональные элементы (коэффициенты полинома) можно помещать по мере вычисления на место исходных значений f k. Следовательно, в программе можно ограничиться только двумя массивами: Х -для узлов и F -для значений аппроксимируемой функции. Массив F по мере выполнения программы будет заполняться разделенными разностями, и в конце концов в нем будут получены коэф­фициенты полинома A 0, A 1, A 2,..., A n. Таким образом, данный метод аппроксимации, так же как и в случае канонического полинома, дает в качестве результата коэффициенты интерполяционного полинома. После определения коэффициентов A 0, A 1, A 2,..., A n вычисление значений полинома Ньютона при конкретных аргументах x рекомендуется выполнять также по схеме Горнера, для чего полином (5.8) надо преобразовать к виду: P n(x)= A 0+(x-x 0)(A 1+(x-x 1)(A 2+...+(x-x n-1) A n)...)). На рис.5.4 представлена блок-схема вычисления коэффициентов полинома Ньютона и его значений в точках интерполяции по схеме Горнера.

 

 

Метод наименьших квадратов

Последовательность действий при аппроксимации экспоненциальной зависимостью (5.21) выглядит так:

1) вычисление логарифмов значений аппроксимируемой функции ;

2) вычисление коэффициентов и по формулам (5.23);

3) вычисление коэффициентов c и d, по формулам (5.24);

4) вычисление значений по формуле (5.21).

В практике обработки экспериментальных данных могут быть ситуации, когда применение лагранжевой аппроксимации (полиномиальной или сплайновой) не оправдано или в принципе невозможно. Первым примером такой ситуации могут служить случаи, когда набор экспериментальных данных был получен со значительной погрешностью, либо на измеряемую (зависимую) величину влияли некоторые дополнительные, не учитываемые факторы. Для демонстрации этой ситуации на рис.5.5 представлены экспериментальные точки, истинная неизвестная кривая f(x) и аппроксимирующая кривая (x), полученная одним из методов лагранжевой аппроксимации. Второй пример, представленный на рис.5.6, демонстрирует ситуацию, когда экспериментальные замеры в каждом узле проводились неоднократно и, вследствие погрешности измерительных приборов либо каких-либо других факторов, дали разные результаты. В этом случае применение лагранжевой аппроксимации в принципе невозможно, так как каждому узлу xi соответствует несколько разных значений fi.

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

Рис.5.5.

Рис.5.6.

Как и в описанных выше методах аппроксимации считаем известными значения экспериментальных данных в узлах f (x _i) = f _i и через (x) обозначим непрерывную аппроксимирующую функцию. В узлах значения функций f (x) и (x) будут отличаться на величину _i = f (x _i) - (x _i). Отклонения _i могут принимать как положительные, так и отрицательные значения. Чтобы не учитывать знаки, возведем каждое отклонение в квадрат, а для оценки близости функций f (x) и (x) возьмем сумму этих квадратов

Q = = .

(5.11)

Метод построения аппроксимирующей функции (x) из условия минимума величины Q называется методом наименьших квадратов (далее - МНК).

Наиболее распространен способ выбора функции (x) в виде линейной комбинации

(x) = с 0 0(x) + с 1 1(x) + … + сm m(x),

(5.12)

где ₀(x), ₁(x), …, _m(x) - базисные функции; ;

с _0, с _1, …, с_m - коэффициенты, определяемые при минимизации величины Q.

При обработке экспериментальных данных, полученных с погрешностью в каждой узловой точке, обычно начинают с аппроксимации функцией , представленной одной-двумя базисными функциями. После определения коэффициентов с _k вычисляется ве­ли­чина Q по формуле (5.11). Если окажется, что , то необходимо расширить базис добавлением новых базисных функций . Расширение базиса необходимо про­­должать до тех пор, пока не выполнится условие .

Выбор конкретных базисных функций зависит от свойств аппроксимируемой функции f (x), таких, как периодичность, экспоненциальный или логарифмический характер, симметричность, наличие асимптот и т.д. Различные варианты базисов рассматриваются достаточно подробно в [1]. Здесь рассмотрим лишь частный случай, когда аппрксимирующая функция представлена двумя базисными функциями, т.е.

 

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

, ,

где символ "надчеркивание" обозначает среднее значение:

где		- среднее значение узловых точек аппроксимации;
		- средний квадрат значений узловых точек аппроксимации;
		- среднее значение аппроксимируемой функции в узловых точках;
		- среднее произведений значений аппроксимируемой функции в узловых точках на значения соответствующих узловых точек.

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

.

(5.21)

Прологарифмируем значения аппроксимируемой функции f (x) в узловых точках:

, i=0,1,2,…,n

и для реализации по формулам (5.20) найдем линейную аппроксимирующую функцию

(5.22)

Формулы (5.20) выглядят при этом так:

.

(5.23)

Чтобы теперь осуществить переход от функции к функции , надо пропотенцировать обе части равенства (5.22):

Величину обозначим , а коэффициенты с и d, входящие в (5.21), вычисляются по формулам:

.

(5.24)

 

Методы прямоугольников

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

В любом случае значение частичного интеграла определяется как произведение длины интервала разбиения на выбранную константу, т.е. как площадь прямоугольника. В зависи-мости от способа выбора аппроксимирующей константы различают методы левых, средних или правых прямоугольников (рис.6.4).

Левые	Средние	Правые

Введем следующие обозначения: точку a на оси OX обозначим через x ₀, точку b - через x _n, а точки разбиения промежутка [ a,b ] - через x ₁, x ₂,..., x _n-1. Предполагается, что длина интервала разбиения постоянна на всем [ a,b ]. Обозначим ее через h:

; x _i = x _i-1 + h, i =1,2,..., N.

Тогда в методе левых прямоугольников площадь каждого i -го прямоугольника

S i = h f (x i), i = 0,1,2,..., n -1,

(6.2)

а для всего промежутка [ a,b ]:

Аналогично, в методе правых прямоугольников

S i = h f (x i), i = 1,2,..., n; .

(6.3)

и в методе средних прямоугольников

S i = h ), i = 0,1,2,..., n -1; ,

(6.4)

где , i = 0,1,2,..., n -1.

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

На рис.6.5. приведена блок-схема вычисления определенного интеграла методом средних прямоугольников.

Рис.6.5. Алгоритм метода средних прямоугольников

Алгоритмы для методов левых и правых прямоугольников отличаются от изображенного на рис.6.5 лишь одним блоком (он выделен жирной линией). Для метода левых прямоугольников здесь должно стоять X=A, для метода правых прямоугольников должно быть X=A+h.

Оценим точность этих методов. В методе средних прямоугольников для каждого интервала разбиения получаем c учетом выражения для S _i в (6.4):

.

(6.5)

Для оценки R _i разложим функцию f (x) в ряд Тейлора около средней точки

(6.6)

В малой окрестности точки этот ряд с высокой точностью представляет функцию f (x) при небольшом количестве членов разложения. Поэтому, подставляя под знак интеграла вместо f (x) ее тейлоровское разложение (6.6) и интегрируя его почленно, можно вычислить интеграл с любой наперед заданной точностью. T.е. точное значение интеграла на интервале [ x _i, x _i+1] рав­но:

Подставим пределы интегрирования:

или, так как :

Все члены полученного при интегрировании ряда, имеющие (x - x _i) в четной степени, обращаются в нуль. Поэтому получаем:

(6.7)

Сравнивая (6.5) и (6.7), можно записать выражение для погрешности R _i:

При малой величине шага интегрирования h основной вклад в значение R _i дает первое слагаемое, которое называется главным членом погрешности вычисления интеграла на интервале [ x _i, x _i+1] и обозначается R _0i:

.

(6.8)

Главный член полной погрешности для интеграла на всем промежутке [ a,b ] определится как сумма:

.

(6.9)

Здесь использован тот же метод средних прямоугольников, но для функции .

Степень шага h, которой пропорциональна величина R ₀, называется порядком метода интегрирования. Как видно из (6.9 ), метод средних прямоугольников имеет второй порядок.

Аналогично проведем оценку метода левых прямоугольников. Разложим подынтегральную функцию в ряд Тейлора в окрестности точки x = x _i:

Интегрируя это разложение почленно на интервале [ x _i, x _i+1] получаем

Здесь первое слагаемое есть приближенное значение интеграла, вычисленное по методу левых прямоугольников (см. формулу (6.2)), а второе слагаемое является главным членом погрешности:

.

(6.10)

Тогда на всем промежутке интегрирования [ a,b ] главный член погрешности R ₀ получается суммированием частичных погрешностей R _0i:

,

(6.11)

т.е. метод левых прямоугольников имеет первый порядок. Метод правых прямоугольников также имеет первый порядок.

Сравнение (6.9) и (6.11) показывает, что метод средних прямоугольников имеет меньшую погрешность по сравнению с методом левых или правых прямоугольников и за счет коэффициента в знаменателе (24 > 2), и за счет интеграла от производной, т.к. для большинства функций выполняется неравенство

.

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

Метод трапеций

В этом методе подынтегральная функция f (x) на интервале [ x _i, x _i+1] заменяется полиномом первой степени, т.е. наклонной прямой линией. Обычно эта прямая проводится через значения f (x) на границах интервала (рис.6.6). В этом случае приближенное значение частичного интеграла определяется площадью трапеции:

Рис.6.6. Геометрическая интерпретация метода трапеций

, т.е. , а численное значение интеграла на всем [ a,b ] . Это вычислительная формула метода трапеций.

(6.12)     (6.13)

Блок-схему алгоритма метода трапеций предлагается студентам разработать самим.

Оценим погрешность R _i. Для этого разложим функцию f (x) в ряд Тейлора около точки x _i:

(6.14)

Тогда

(6.15)

С помощью разложения (6.14) вычислим подынтегральную функцию в точке x _i+ h:

откуда

(6.16)

Подставляя произведение (6.16) в выражение (6.15), получим

(6.17)

Сравнивая (6.12) и (6.17), получаем выражение для главного члена погрешности частичного интеграла

.

Тогда главный член полной погрешности метода трапеций имеет вид

,

(6.18)

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