Лабораторная работа 5. Методы интерполирования функций

Узел интерпо-ляции	Задание
№	х	e x	e - x	sh x	ch x	sin x	cos x	ln x
	1,00	2,7183	0,3679	1,1752	1,5431	0,8415	0,5403	0,0000
	1,01	2,7456	0,3642	1,1907	1,5549	0,8468	0,5319	0,0100
	1,02	2,7732	0,3606	1,2063	1,5669	0,8521	0,5234	0,0198
	1,03	2,8011	0,3570	1,2220	1,5790	0,8573	0,5148	0,0296
	1,04	2,8292	0,3535	1,2379	1,5913	0,8624	0,5062	0,0392
Узел интерпо- ляции	Задание
№	X	e x	e - x	sh x	ch x	sin x	cos x	ln x
	1,05	2,8577	0,3499	1,2539	1,6038	0,8674	0,4976	0,0488
	1,06	2,8864	0,3465	1,2700	1,6164	0,8724	0,4889	0,0583
	1,07	2,9154	0,3430	1,2862	1,6992	0,8772	0,4801	0,0677
	1,08	2,9447	0,3396	1,3025	1,6421	0,8820	0,4713	0,0770
5	1,09	2,9743	0,3362	1,3190	1,6552	0,8866	0,4625	0,0862

Окончание табл. 4.1

Узел интерпо-ляции	Задание
№	X	e x	e - x	sh x	ch x	sin x	cos x	ln x
	1,10	3,0042	0,3329	1,3356	1,6685	0,8912	0,4536	0,0953
	1,11	3,0344	0,3296	1,3524	1,6820	0,8957	0,4447	0,1044
	1,12	3,0649	0,3263	1,3693	1,6956	0,9001	0,4357	0,1133
	1,13	3,0957	0,3230	1,3863	1,7093	0,9044	0,4267	0,1222
	1,14	3,1268	0,3198	1,4035	1,7233	0,9086	0,4176	0,1310

Вспомогательные материалы

1. Пример выполнения работы (п.1)

Пусть задана функция y = sin xна интервале [1,15; 1,19 ]:

При задании с шагом 0,2:

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

.

Вычислим по данной формуле значение в точке x = 1,16 . Имеем

x(0) = 1,15 ; x(1) = 1,17 ; x(2) = 1,19; y(0) =0,9128 ; y(1) = 0,9208; y(2) = 0,9284.

Подставляя данные в выражение для L(x), получаем:

.

Полученное значение точно совпадает с табличным для x = 1,16.

Оценим погрешность интерполирования по формуле (4.18). Имеем n = 2.

Максимальное значение косинуса не превосходит 1. Следовательно, М _n+1 равно 1.

.

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

Разработка программы на Mathcadе трудностей не представляет.

Пункт второго задания выполняется аналогично.

Пример Mathcad (программы для построения кубического сплайна):

Контрольные вопросы

1. Значение интерполяционного полинома в узловых точках?

2. Какова максимальная степень многочлена Лагранжа?

3. Какова максимальная степень многочлена Ньютона?

4. Как оценивается погрешность интерполяции?

5. В чем состоит различие сплайн-интерполирования от обычной интерполяции?

ЧИСЛЕННОЕ ИНТЕГРИРОВАНИЕ

Задача приближенного вычисления определенного интеграла (на отрезке или по многомерной области) фактически разбивается на две самостоятельные подзадачи: 1. Интегрирование таблично заданной функции (полученной, например, при проведении лабораторного эксперимента). В таком случае информация о гладкости подынтегральной функции отсутствует, весьма ограничены возможности в выборе узлов интегрирования. Для этой задачи наиболее эффективными будут квадратурные формулы интерполяционного типа и правило Рунге оценки погрешности. 2. Подсчет значения определенного интеграла от известной функции. При этом самая ресурсоемкая операция с точки зрения вычислений – подсчет значения функции. Желательно построить численный метод, позволяющий получать как можно более высокую точность при наименьшем количестве вычислений. При этом выбор узлов квадратурных формул целиком определяется исследователем. В этом случае наиболее эффективными окажутся квадратурные формулы типа Гаусса.

Квадратурные формулы интерполяционного типа (формулы Ньютона-Котеса). Геометрический смысл определенного интеграла функции f(x) заключается в площади фигуры, образованной этой функцией и осью OX, поэтому самый простой способ посчитать определенный интеграл от "хорошей", то есть гладкой, функции ­­– применить формулу прямоугольников или трапеций. C помощью этих формул площадь упомянутой искомой фигуры подсчитывается как сумма элементарных прямоугольников или трапеций, множеством которых заменяется подынтегральная функция f(x).

Для подсчета интеграла интервал интегрирования [a,b] разбивается на n равных частей (отрезков) точками (h – шаг разбиения ). При этом точка x_i может выбираться, например, как начало каждого элементарного отрезка, либо как его центр. Значение функции f(x) в точках разбиения х_i обозначим через у_i. Непрерывная подынтегральная функция заменяется сплайном – кусочно-полиномиальной функцией S(x), аппроксимирующей данную функцию. Очевидно, что при стремлении h к 0, множество прямоугольников (или трапеций) стремится к искомой фигуре, образованной подынтегральной функцией, а численный результат – к истинному значению интеграла.

Интегрируя функцию S(x) на отрезке [a,b], получим в зависимости от S(x) следующие квадратурные формулы.

1. Формула прямоугольников.

Если на каждой части [x_i-1,…,x_i] (i=1,2,…,n) деления отрезка [a,b] функцию f(x) заменить функцией, принимающей постоянное значение, равное, например, значению функции f(x) в серединной точке i-й части , то есть на каждой части отрезка функция имеет вид прямоугольника с шириной и высотой f(x_i), то функция S(x) будет иметь ступенчатый вид:

, х Є [x_i-1,…,x_i], (i=1,2,…,n).

,

Отсюда получаем квадратурную формулу прямоугольников:

. (5.1)

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

2. Формула трапеций (рис. 5.2).

Если предположить, что функция f(х) на отрезке интегрирования [a, b] достаточно близка к линейной, то данную функцию можно заменить на каждом отрезке [x_i-1,…,x_i] ее линейной интерполяцией по точкам (x_i-1 , у_i-1) и (x_i , у_i), то есть ломаной линией с вершинами в концах элементарных отрезков, на которые разбивается интервал интегрирования (аппроксимация искомого интеграла множеством элементарных трапеций с высотой (b – a) и основаниями f(a) и f(b)).

В результате получим кусочно-линейную функцию

, х Є [x_i-1,…,x_i], (i=1,2,…,n),

здесь у_i=f(x_i). .

Тогда формула трапеций

, (i=1,2,…,n)

Рис. 5.2

3. Формула Симпсона (рис. 5.3) – это формула парабол, которую можно получить при условии, что сплайн S(x), аппроксимирующий подынтегральную функцию f(x), представляет собой непрерывную функцию, составленную из примыкающих парабол. Пусть парабола на отрезке [x_i-1,…,x_i] проходит через точки (x_i-1, у_i-1), (x_i-1/2 , у_i-1/2) и (x_i , у_i). Используя построение интерполяционного многочлена Лагранжа второго порядка на данном отрезке, получим сплайн:

;

х Є [x_i-1,…,x_i], (i=1,2,…,n).

Для дальнейших преобразований введем переменную t Є [0;1] с помощью равенства х = x_i-1 +ht . Значениям t (равным 0, ½, 1) соответствуют значения х (равные x_i-1 ; x_i-1/2 ; x_i ) Выразим сплайн S(x) через новую переменную

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

Учитывая, что имеем:

.

В результате приходим к квадратурной формуле парабол:

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

.

Рис. 5.3.Оценка погрешности квадратурных формул

Рассмотрим интеграл по элементарному отрезку

,

где х₀ Є [x_i,x_i+1]— некая опорная точка, тогда для приближенного значения интеграла верно

.

Коэффициенты ξ, η, ... зависят от производных f'(х₀), f''(х₀), ....

С другой стороны, любую из рассмотренных квадратурных формул можно представить так:

Заменяя в этой формуле значения функций f в точках f_i, f_{i + 1/2}, f_{i + 1}, ее разложением по формуле Тейлора получим:

,

где х₀Є [x_i,x_i+1].

Сравнивая разложения для , легко заметить, что вместе с первым слагаемым совпадают и другие слагаемые до (m - 1)-го порядка, так что ξ = ξ₁, η = η₁, ...

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

,

где v – константа.

Если просуммировать локальные погрешности по всем интервалам [ x_i, x_{i + 1}], то получим оценку погрешности квадратурной формулы по всему отрезку [a, b]:

,

где h=(maxh_i)/i на неравномерной сетке, или h= (b-a)/n на равномерной. Число m называется порядком точности квадратуры.

Если подынтегральная функция имеет непрерывную вторую производную, то оценка погрешности для формулы:

1) прямоугольников

;

2) трапеций

.

Если подынтегральная функция имеет непрерывную производную четвертого порядка, то справедлива следующая погрешность формулы Симпсона:

.

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

Квадратурные формулы Гаусса

Поскольку формулы Ньютона-Котеса являются интерполяционными, очевидно, что они не могут успешно использоваться для получения формул высокой точности по причине неустойчивости интерполяционного процесса для многочленов высокого порядка. По этой причине обычно используются полиномы степени от нуля до трех (соответственно, формулы прямоугольников со средней точкой, трапеций, Симпсона). Вычисление с их помощью интегралов от функций, обладающих высокой степенью гладкости, например, близким к полиномам высокой степени, представляется нерациональным. В выражение для погрешности этих формул входят первая, вторая или четвертая производные. Погрешность определяется низким порядком производной при высокой степени гладкости интегрируемой функции. Этих недостатков лишены квадратуры Гаусса.

Формулировка задачи построения квадратурных формул, поставленная Гауссом, такова.

Для заданного количества точек найти такое расположение узлов и веса c_i, чтобы квадратурная формула

была точной для полиномов как можно более высокой степени.

Пусть отрезок интегрирования [a,b] непрерывной функции f(x) разбит на n равных частей точками (h – шаг разбиения) . Обозначим с помощью S(x) сплайн-функцию, аппроксимирующую подынтегральную функцию f(x).

Пусть на каждой части разбиения [x_i-1 ,x_i] (i=1,2,…,n) расположено m узлов [x_i1,…,x_im], в которых подынтегральная функция f(x) на каждой i-й части аппроксимируется многочленом степени p, х Є [x_i-1,…,x_i], (i=1,2,…,n). При этом на многочлен накладываются два ограничения:

а) значения многочлена и подынтегральной функции равны в узлах интерполяции: S_i(x_ij)=f(x_ij) (i=1,…,n; j=1,…,m);

б) определенный интеграл от функции на отрезке [x_i-1, x_i] выражается через значения подынтегральной функции f(x_ij) в узлах в виде их линейной комбинации

. (5.2)

Квадратурные формулы Гаусса для выбранной степени p сплайна будут определены, если из условий а) и б) удастся найти m неизвестных коэффициентов c_i и координаты m узлов x_ij ( j=1,…,m).

Задача решается одновременно для всех n частей разбиения отрезка [a,b], если выразить х Є [x_i-1,…,x_i], (i=1,2,…,n) через переменную t Є [-1;1]:

;

Положим ; ; .

Тогда (i=1,…,n; j=1,…,m) и соотношение (5.2) перепишем в виде

. (5.3)

Выведем квадратурную формулу Гаусса с тремя узлами (m=3). Для этого необходимо определить шесть величин: . Функция -многочлен степени p, общий вид которого

. (5.4)

Подставив соотношение (5.3) и (5.4), а также учитывая, что (j=1,2,3,…,p), получим тождество относительно коэффициентов a_k (k=0,1,2,…,p):

.

Шесть неизвестных будут определены однозначно из системы шести уравнений. В общем случае степень p аппроксимирующего многочлена всегда является нечетным числом и связана с числом узлов m соотношением p=2m-1. В частности, для трех узлов имеем многочлен пятой степени.

Множители при a_k в левой части тождества вычисляются так:

Приравняем подобные выражения в левых и правых частях тождества при одинаковых коэффициентах ₀, _{1, 5} получим:

; ;

; ; (5.5)

; .

Эту систему позволяет упростить следующее свойство ее решения: неизвестные t_i (i=1,2,…,m) системы 2m уравнений вида

,

где k=0,1,…,2m-1 являются нулями многочлена Лежандра ; нули принадлежат интервалу (-1;1) и расположены симметрично относительно середины интервала. В данном случае m=3 и . Находим нули многочлена из уравнения 5t³-3t=0. Подставим корни уравнения , t₂=0, в формулу(5.5), получим систему трех линейных уравнений относительно переменных :

, , .

Теперь подставим найденные значения в соотношение (5.2)

,

где , , .

Итак, квадратурная формула Гаусса с тремя узлами записывается в виде

.

Если подынтегральная функция имеет непрерывную производную шестого порядка, то для оценки погрешности формулы Гаусса с тремя узлами можно использовать неравенство:

.

При вычислении интеграла до достижения заданной точности ε методом двойного пересчета условие окончания вычислений имеет вид

, где k=2m (m – число узлов в квадратурной формуле Гаусса). Полагают, что с точностью ε. Формулы Гаусса обеспечивают высокую точность уже при небольшом количестве узлов.

Блок-схемы решения задачи рассмотренными методами приведены на рис. 5.4 и 5.5.