Глава 3. ЧИСЛЕННОЕ ДИФФЕРЕНЦИРОВАНИЕ И ИНТЕГРИРОВАНИЕ
Существует много машинных методов интегрирования и дифференцирования. Метод, более всего подходящий для данной задачи, в значительной степени зависит от информации о соответствующей функции. Особый интерес с точки зрения вычислительной математики при этом вызывают задачи для одной функции f(x) одного действительного переменного х Î [a, b], которые условно подразделяют на четыре категории.
-
Значения функции f(x) заданы только на фиксированном конечном множестве точек хi интервала [a, b].
-
Функция f(x) определена и может быть вычислена для любого действительного х из интервала [a, b], однако первообразной для нее не существует.
-
Определение функции может быть аналитически продолжено на комплексные значения х.
-
Имеется явная формула, пригодная для символического манипулирования.
Функции, входящие в первую категорию, чаще всего появляются в результате некоторого эксперимента для заданных хi, которые, как правило, распределены неравномерно или расположены в виде таблицы.
Поэтому ясно, что для функций первых двух категорий численное дифференцирование является более трудной вычислительной задачей, чем численное интегрирование. Очевидно причина в том, что математическая операция численного дифференцирования увеличивает любую ошибку, присутствующую в данных, в то время как численное интегрирование обычно сглаживает и уменьшает такие ошибки. Если значения функции известны или могут быть вычислены с большой точностью и если требуются производные относительно невысоких порядков, то часто методы, основанные на интерполяции сплайнами или полиномами, дают вполне удовлетворительные результаты. Однако если нужны производные высокого порядка или если значения функции зашумлены, то результаты могут быть, мягко говоря, неточными.
Примером функций, входящих в третью категорию, могли бы служить сложные выражения, составленные из тригонометрических или других элементарных функций. Для задач этого типа следует извлекать выгоду из комплексного расширения, если оно доступно. Тогда производные функции f(x) можно выразить через комплексные контурные интегралы и сглаживающий эффект численного интегрирования можно использовать для получения хороших приближений к производным высокого порядка [Линес, Моулер, 1967; Линес, Санде, 1971].
В основном настоящая глава посвящена методам работы с функциями 1, 2 и 3-й категории.
Задачи 4-й категории в данной работе не рассматриваются, но сведения о данных такого типа можно найти в работах Мозеса (1972).
Для того чтобы избежать путаницы с численным интегрированием обыкновенных дифференциальных уравнений, для численного приближения определенных интегралов используется понятие квадратура.
§ 1. ЧИСЛЕННОЕ ДИФФЕРЕНЦИРОВАНИЕ С ПОМОЩЬЮ ИНТЕРПОЛЯЦИОННЫХ ФОРМУЛ НЬЮТОНА И ЛАГРАНЖА
Пусть некоторая функция F(х) задана таблично в n+1 точке, расположенной на интервале [а, b]. Требуется вычислить первую и вторую производные данной функции в заранее определенных точках.
В качестве аппроксимирующей функции выберем интерполяционный многочлен. Если узлы интеpполяции расположены не равномерно, то таким многочленом могут быть полиномы Лагранжа или Лежандра, а если равномерно, то лучше использовать полином Ньютона, так как он дает меньшую вычислительную погрешность по сравнению с другими.
Пусть узлы хi, в которых известны значения функции F(хi), расположены равномерно на интервале [а, b], т.е. хi+1-хi = h, где i = 0, 1, ..., n. Тогда в качестве интерполирующей функции Рn(х) выберем полином Ньютона
,
где q = (х - х0) / h; h = хi+1 - хi.. Здесь использована запись первой интерполяционной формулы Ньютона (см. § 2, гл. 2). Раскроем скобки, приведем подобные, получим
(3.1)
Заметим, что .
Продифференцируем f(х) из выражения (3.1) и с учетом последнего замечания получим выражение для первой производной
(3.2)
Продифференцируем последнее выражение еще раз:
. (3.3)
Здесь учитывалось, что
.
Аналогично можно определить и следующие производные функции f(х). Но при вычислении производных более высоких порядков следует учитывать большее количество членов ряда Ньютона для обеспечения точности результата.
Если производные функции f(х) надо определить в узлах интерполяции, т.е. в х = хi , то тогда формулы (3.2), (3.3) для первой и второй производных заметно упрощаются, так как в этом случае q º 0:
f ’(x) = (D1 - D2/2 + D3/3 - D4 /4 + ... ) / h ; (3.4)
f ’’(x) = (D2 - D3 +11 D4 /12 -5 D5 / 6 + ... ) / h2 . (3.5)
Погрешность в определении производных можно приближенно оценить как [Плис, Сливина, 1983]
.
Для того чтобы получить значения производных функции f(х) в точках, расположенных в конце таблицы, следует, как и в § 2, гл. 2, воспользоваться второй интерполяционной формулой Ньютона. Повторив аккуратно и последовательно все рассуждения, получим для первой производной следующее выражение:
(3.6)
Если же для равноотстоящей системы узлов хi построить полином Лагранжа, то, записав его в виде
,
где q = (х - х0)/h - шаг интерполяции; h = хi+1 - хi - расстояние между узлами интерполяции; yi = f(хi)- значения функции f(х) в заданных узлах хi, и c учетом, что dх/dq = h, получим выражение для первой производной:
.
Для оценки погрешности последнего выражения можно воспользоваться формулой, приведенной в работе Плис и Сливиной (1983)
.
В заключение отметим, что рассматриваемые здесь формулы численного дифференцирования являются менее точными по сравнению с интерполяционными (см. § 1 и 2 в гл. 2). Однако их применяют в практических расчетах достаточно часто, так как они удобны, просты, не требуют специальных программ. И, что особенно важно, позволяют быстро сделать предварительные (прикидочные) расчеты.
Процедуры new1, new2, newton из § 1 и 2 (гл.2) вполне пригодны для вычисления производных. Но с учетом решаемой задачи они должны быть несколько модифицированы, хотя формальные параметры остаются у процедур без изменений.
{***Первая интерполяционная формула Ньютона**}
function new1 (key:integer;x,y,r1,r2,r3,r4: mas;
x1:real):real;
label 30;
var i,j,k : integer; q : real;
begin
i := 1;
j := 14;
if x1 < x[i] then
goto 30;
{ ****Метод двоичного поиска интервала****}
repeat
k := (i+j) div 2;
if x1< x[k] then
j := k;
if x1>=x[k] then
i := k;
until j <= i+1;
30:
q := (x1 - x[i]) / (x[2]-x[1]);
{ ****Выбор типа работы процедуры****}
case key of
0:new1:=y[i]+q*(r1[i]+(q-1)*(r2[i]/2.0+
(q-2)*r3[i]/6.0));
1: new1 := r1[i] + 0.5*(2*q-1)*r2[i] +
((3*q-6)*q+2)*r3[i]/6+