Лекция 5.

Кусочная интерполяция. Численное дифференцирование (ЧД).

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

Пример.

Необходимо вычислить f (x) для x [xi ; xi+1 ].

Кусочно-линейная интерполяция

Для вычисления f (x) используется линейное приближение

f (x) ≈ fi + fi +1h− fi (x − xi ), где h = xi+1 − xi .

Очевидно, что при подстановке x = xi или x = xi+1 получается тождество.

Кусочно-квадратичная интерполяция.

Привлекается еще одна табличная точка ( xi−1 или xi+2 ), и строится поли-

ном второй степени:

тождество.

Кусочно-кубическая интерполяция.

38

Нетрудно убедиться, что при fi−1 , fi , fi+1 , fi+2 в последних двух форму-

лах стоят биномиальные коэффициенты.

Более высокие степени ИП для КИ, как правило, не используют.

Узлы интерполяции при кусочной интерполяции берутся вблизи интересующего нас узла xi . (Рис. 5.1)

Рис. 5.1 Узлы интерполяции при кусочной интерполяции

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

Возможные обобщения приближения функций с помощью интерполяции.

По аналогии с ИП можно искать приближение таб-

лично заданной функции в виде обобщенного ИП:

Pn (x)= a0φ0 (x)+ a1φ1(x)+... + anφn (x)

по системе линейно независимых функций {φk (x):k = 0,1,...,n}.

Исходя из условий интерполяции (совпадения значения полинома с табличными значениями функции), для неопределенных коэффициентов {ai } по-

лучим систему линейных уравнений:

a0φ0 (xi )+ a1φ1(xi )+... + anφn (xi )= fi , i = 0,1,..., n .

Для существования и единственности решения необходимо, чтобы детерминант удовлетворял условию

Например, периодическую функцию может оказаться удобно приближать

39

метрическая интерполяция. Часто используется, если функция разложима в ряд Фурье.

Если в табличных точках заданы не только значения функции, но и ее производные, то для приближения функции используются полиномы Эрмита. Для построения используется следующие условия:

•В табличных точках полином принимает заданные значения;

•Производные от полинома также принимают заданные значения.

Втаких случаях полиномы Эрмита обеспечивают лучшее приближение сравнительно с обычной интерполяцией.

Интерполяция кубическими сплайнами.

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

Строятся они следующим образом. На каждом отрезке [xk ; xk +1 ] ищется полином третьей степени P3k (x)= ak + bk x + ck x2 + dk x3 со своими коэффициен-

тами. Потребуем, чтобы в табличных точках полином принимал заданные значения:

.

Кроме того, необходимо следующее условие:

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

Отметим еще задачу обратной интерполяции.

40

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

жа и в форме Ньютона выписывается стандартным образом.

Численное дифференцирование.

Итак, пусть в точках xi , i = 0,1,..., n известны значения функции

f (x) :{fi ,i = 0,1,..., n}.

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

Дифференцируя полином нужное число раз, получаем требуемую форму-

лу.

Например, f (x) = Pn (x)+ Rn (x),

Погрешность формулы характеризуется k -той производной от ошибки

интерполяции: ∂k Rn (x).

∂xk

Замечание.

Анализ погрешности формул численного дифференцирования, опирающийся на оценку остаточного члена затруднен тем, что его выражение, как следует из его вывода (Лекция 3)

Rn (x)= f(n++1(x)′) ωn+1(x),

n 1 !

n

где ωn+1(x)= (x − x0 )(x − x1 )...(x − xn )= ∏(x − xi )

i=0

содержит неявную зависимость x′ = x′(x) [x0 ; xn ] .

41

Поэтому будем иначе анализировать точность формул численного дифференцирования.

Примеры.

Приближение производных функции f (x) в окрестности табличной точки

x = xi . Пусть xi = x0 +ih , h = const .

Приближение полиномом первой степени

f (x) ≈ P1(x)= fi + fi+1h− fi (x − xi )

совпадает с функцией в точках x = xi и x = xi+1 .

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

f (x) ≈ P1(x)= fi + fi −hfi−1 (x − xi )

совпадает с функцией в точках x = xi−1 и x = xi .

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

Приблизим f (x) в рассматриваемой окрестности точки полиномом второй степени:

(совпадает с табличными значениями функции в т. x = xi−1 , x = xi , x = xi+1 ),

отсюда:

При x = xi :

42

Получено так называемое центральное разностное отношение.

Дифференцируя 2 раза, получим

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

Причем порядок стремления к нулю есть O(h), т.е. пропорционально пер-

вой степени величины h . Это первый порядок аппроксимации производной разностным отношением.

Вычитая из формулы (1) формулу (2), получим

43

Это аппроксимация второго порядка (малости по h ). Складывая формулы (1) и (2) , получим для второй производной:

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

производных более высокого порядка.

Неустойчивость формул численного дифференцирования. Рассмотрим влияние погрешности входных данных на результат вычисле-

ния производных по формулам ЧД.

Погрешность формулы

методическая неустранимая ошибка погрешность

44

Рис.5.2 График функцииΦ(h)

Ясно, что при h → 0 , например, h δ, неустранимая погрешность Ο(1).

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

Возьмем значение функции в 10 точках с погрешностью входных данных

σ = 0.01 , т.е. ~pi = pi +σ rand(1,1) .

Построим ИП, например, используя функцию MATLAB y=polyval( ~p, x ), где ~p - коэффициенты полинома в точках x .

Используя функцию MATLAB k = polyder( ~p) , найдем коэффициенты поли-

нома, получающиеся при дифференцировании данного.

45

Ясно, что количество действий может быть уменьшено. Можно последо-

вательно вычислять и запоминать x2 , x3 ,..., xn , тогда потребуется (n −1) опера-

ция умножения. Далее вычисляются a j x j - n операций умножения и n опера-

ций сложения. Количество действий Φ2 = (2n −1) + n . Уже при n > 2 Φ2 < Φ1 .

Можно еще эффективнее:

Pn (x) = (...(an x + an−1)x + an−2 )x +...)x + a0 .

Для вычисления каждой скобки 1 операция умножения и 1 операция сложения. Количество действий Φ3 = n + n .

Φ3 < Φ2 < Φ1 при n > 2 .

Количество арифметических действий - важнейшая характеристика численного метода.

Часто оценивают не количество, а просто порядок количества арифметических операций.

46

Φ2 = 3n +O(n)

Φ3 = 2n

Заметим, что Pn (x) определяется величинами a0 ,a1,..., an и x .

Поэтому в общем случае для вычисления Pn (x) требуется не менее n дей-

ствий – оценка снизу.

Т.о. Φ2 ,Φ3 методы являются оптимальными по порядку, т.к. Φ2 ,Φ3 = O(n) и

для каждого метода Φ ≥ n .

Методы с числом действий O(n) называются экономичными.

47

Лекция 6. Численное интегрирование.

b

Определенный интеграл ∫ f (x)dx возможно точно найти по формуле Нью-

a

b

тона - Лейбница ∫ f (x)dx = F (a) − F (b) только в тех редких случаях, когда инте-

a

грал является табличным.

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

Общая схема простейших квадратурных формул.

Геометрический смысл определенного интеграла ∫b f (x)dx - площадь.

a

Для приближенного вычисления этой площади разобьем отрезок [a,b] на некоторое число N отрезков, так чтобы a = x0 < x1 <... < xn = b .

Определение.

xn - назовем узлами интегрирования, hn = xn+1 − xn - шагами интегрирования, fn = f (xn ) .

В частном случае hn = const;h = b N− a .

Затем заменим функцию y = f (x) кусочно-интерполяционным многочле-

ном:

Pk (x) = Pk (x, f , x0 , x1,..., xn ) .

Функция Pk (x) на каждом отрезке [xn , xn+1] есть многочлен степени не вы-

ше k . Он выбирается по точкам, ближайшим к данному отрезку x [xn , xn+1] .

Напомним формулы кусочной интерполяции (Из лекции 5). Кусочная интерполяция нулевого порядка

48

f (x) ≈ fn или f (x) ≈ fn+1 , или f (x) ≈ fn+12 ≡ f (xn + h2n ) (если есть значение функции в точке n + 12 ).

Кусочно-квадратичная интерполяция

(используется еще одна точка x = xn−1, xn , xn+1 ).

Таким образом, у нас есть интерполяционный полином степени не выше

k .

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

Определение.

Квадратурные формулы ( ) с равноотстоящими узлами называются фор-

49

ухудшается. Геометрическая интерпретация на Рис. 6.1

Рис. 6.1 Геометрическая интерпретация формулы прямоугольников

Формула трапеции. На элементарном отрезке x [xn , xn+1]

P1(x) = fn + fn+1 − fn (x − xn ) xn+1 − xn

(площадь трапеции Рис. 6.2).

Рис. 6.2 Геометрическая интерпретация формулы трапеции

50