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

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

4.1.  Постановка задачи

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

Будем искать решение задачи для функции y = f (x), заданной на интервале [a, b] дискретными значениями с шагом h . Число разбиений N = (b – a) / h (Рис.4.1). Решение задачи будет приближенным. Оно зависит от числа разбиений интервала. Точность вычислений возрастает при увеличении числа разбиений, но до определенного предела, который связан с возрастанием ошибки округления.

Рис.4.1. Дискретное задание функции

Методическая ошибка вычисления может быть представлена в виде

F ( x n) – F n = c hm ,

где F ( x n) – точное решение, F n – приближенное решение,

c – постоянная, не зависящая от h. Методы решения могут иметь
m – порядок точности. При первом порядке точности ошибка решения пропорциональна шагу.

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

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

Представим исходную функцию в виде таблицы 4.1 соответствующих значений x i и y i и конечных разностей Dy i , D2 y i , D3 y i ...

Таблица 4.1.

x	y	D y	D2 y	D3 y	D4 y	D5 y
x0	y0	Dy0	D2 y0	D3 y0	D4 y0	D5 y0
x1	y1	Dy1	D2 y1	D3 y1	D4 y1
x2	y2	Dy2	D2 y1	D3 y2
x3	y3	Dy3	D2 y1
x4	y4	Dy4
x5	y5

В таблице конечные разности различных порядков рассчитываются по формулам:

Dy0 = y1 – y0 , Dy1 = y2 – y1 ,

D2 y0 = Dy1 – Dy0 = y2 – 2y1 + y0 , D2 y1 = Dy2 – Dy1 = y3 – 2y2 + y1,

D3 y0=D2 y1–D2 y0=y3–3y2+3y1-y0, D3 y1=D2 y2–D2 y1 = y4-3y3+3y2-y1,

– – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – –

Применим теорему Лагранжа к функции f (x)

,

где x – число, принадлежащее интервалу [ x i–1 , x i ]. Полагая приближенно x = x i , что соответствует кусочно–линейной интерполяции, получим

,

где h – шаг по оси x. Переходя к конечным разностям, будем иметь

(4.1)

Вычисление второй производной производится по схеме

Производная n порядка рассчитывается по формуле

(4.2)

Для h = 1 формулы упрощаются и производные оказываются равны своим конечным разностям. Этот метод дифференцирования можно считать грубым и получить приемлемую точность можно только путем уменьшения шага.

Пример 4.1.

Продифференцируем численно функцию y = x3 . Составим таблицу конечных разностей при шаге h = 1 .

Таблица 4.2.

x	y	Конечные разности
Dy	D2 y	D3 y	D4 y
– 3	– 27	19	– 12	6	0
– 2	– 8	7	– 6	6	0
– 1	– 1	1	0	6	0
0	0	1	6	6
1	1	7	12
2	8	19
3	27

Точные значения производных y¢ = 3 x2 и y¢¢ = 6 x.

Погрешности определения первой и второй производной в соответствии с таблицей 2 могут достигать 100 %. Для их снижения нужно уменьшать шаг по оси x. Большую точность можно получить, переходя от кусочно–линейной интерполяции к интерполяционным полиномам более высокого порядка.

Запишем интерполяционный полином в форме Ньютона

Найдем первую производную, продифференцировав полином по x.

(4.3)

Найдем вторую производную, продифференцировав (4.3) по x.

(4.4)

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

Производная может быть также вычислена путем дифференцирования интерполяционного полинома в форме Лагранжа. Запишем интерполяционный полином для трех точек x i–1 , x i, x i+1 , расстояние между которыми h .

(4.5)

Дифференцируем полином (4.5) по x:

(4.6)

При x = x i –1

(4.7)

Аналогично можно получить, подставляя в (4.6) x i и x i –1

Вторую производную получим, дифференцируя (4.6) по x

(4.8)

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

, (4.9)

. (4.10)

Пример 4.2.

Рассмотрим функцию

(4.11)

Найдем первую и вторую производные при x = 0.5 по формуле (4.9) и (4.10)

y¢ (0.5 , 0.01) = – 0,1760326626. Точное значение – 0,1760326634

y¢¢ (0.5 , 0.01) = – 0,2640489949. Точное значение – 0, 264048995

На рис. 4.2 показаны графики функции (4.11) и ее производных при
h = 0,01

Рис. 4.2. Графики функции (4.11) и ее производных:

1 – y(x) ; 2 – y¢(x) ; 3 – y¢¢(x) .

Влияние погрешности округления на общую точность дифференцирования рассмотрим на примере функции f(x) = e x

Приближенное значение

Рис.4.3. Точность определения производной от f(x) = e x

при x = 1 для различных значений шага h .

Погрешность вычислений D = f ¢(x) – f ¢ (x)* оценим при x = 1 постепенно уменьшая шаг h . На рис. 4.3 показана зависимость суммарной погрешности вычислений для персонального компьютера IBM Pentium 2 , позволяющего рассчитывать с точностью до 15 разрядов. Из рисунка видно, что минимальная абсолютная ошибка в данном случае составляет 10 –8 .

4.3. Численное интегрирование.

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

-  подинтегральная функция задана дискретным набором данных,

-  интеграл не имеет аналитического решения,

-  аналитическое решение интеграла слишком сложно.

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

В общем виде задача состоит в нахождении величины

(4.12)

Эта величина равна заштрихованной ( рис.4.4) площади между кривой f(x) и осью x. Если разделить промежуток [a, b] на конечное число участков с шагом h , то задача (4.12) сводится к большому количеству более простых задач, число которых

.

Рис.4.4. Численное интегрирование

Задачу (4.12) представим в виде:

Для её решения применим аппроксимацию функции f(x) в виде интерполяционного полинома по формуле Ньютона.

(4.13)

где Dy0 = y1 – y0, D2 y0 = D y1 – D y0 , D3 y0 = D2 y1 – D2 y0…
Dn y0 = Dn–1 y1 – Dn–1 y0 .

Конечные разности запишем в виде таблицы

x	y	D y	D2 y	D3 y	D4 y
x0	y0	Dy0	D2 y0	D3 y0	D4 y0
x1	y1	Dy1	D2 y1	D3 y1
x2	y2	Dy2	D2 y2
x3	y3	Dy3
x4	y4

Введем обозначение

x = x0 + t h ,

тогда

(x – x0) = t h; (x – x1) = h( t – 1) ;

и (4.13) можно представить в виде

(4.14)

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

Рассмотрим частные случаи: 1) t = 0 ; 2) t = 1 . На каждом участке значение функции будет постоянным, и составная сумма, которая определяет определенный интеграл, равна сумме прямоугольников.
В методе левых прямоугольников (t = 0) выбирается значение функции на левом конце участка.

I л = h ( f (a) + f (a + h) + ... + f (b – h)) = h ( y0 + y1 + ...+ y n–1). (4.15)

В методе правых прямоугольников (t = 1) выбирается значения функции на правом конце участка.

I п = h ( y1 + y2 + ...+ y n ) . (4.16)

Точное значение суммы находится между найденными значениями. Погрешность вычисления

D = ½ I л – I п ½ = h½ y0 – y n ½.

В этом случае ошибка пропорциональна первой степени шага h .

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

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

Для первого элементарного участка при x = x0 + t h , y =y0+Dy0t

. (4.17)

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

(4.18)

Оценим точность метода трапеций. Точное решение для элементарного участка имеет вид:

. (4.19)

Воспользовавшись формулой (4.17) для приближенного решения, можно записать

, (4.20)

где d - погрешность вычислений. Разложим в ряд Тейлора F(xi+1)

. (4.21)

Продифференцируем (4.21) по x и учитывая, что F'(xi+1) = f(xi+1), F'(xi) = f(x), F''(xi) = f '(xi) и т. д., получим

. (4.22)

Подставляя (4.21) и (4.22) в (4.20) и решая (4.20) относительно d, получим

Пренебрегая 0(h4), имеем

. (4.23)

Погрешность для всего промежутка интегрирования

, (4.24)

где – среднее значение второй производной. Из (4.22) следует, что методическая ошибка метода трапеций пропорциональна h 2 .

Разновидностью метода трапеций с переменным шагом является метод Ромберга, используемый системой MATHCAD. Он сочетает простоту метода трапеций с высокой точностью вычисления. Для метода трапеций (4.18) справедливо

,

где h = (b – a) / n и x i = a + i h. Уменьшать шаг, увеличивая n можно, используя прогрессию вида

n = 2 k-1 ; k = 1, 2, 3 ...

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

h k = (b – a) / 2 k-1.

Приближенное значение интеграла

,

при k = 1 ,

при k = 2 и т. д.

В то время как k увеличивается, погрешность вычисления R k уменьшается. Окончание расчета по условию .

Метод параболы (Симпсона)

Если в интерполяционном полиноме (4.13) учесть дополнительный член, то

(4.25)

При этом необходимо на каждом участке использовать два интервала и три узловые точки (Рис.4.5).

Рис.4.5. Метод парабол (Симпсона).

При n = 2 таблица конечных разностей имеет вид

x	y	D y	D2 y
x0	y0	Dy0	D2 y0
x1	y1	Dy1
x2	y2

Подстановка (4.25) в интеграл (4.12) дает

(4.26)

Формулу (4.25) называют формулой Кеплера. Она имеет ошибку

.

Записав также формулы для следующих двойных интервалов и сложив их, получим формулу Симпсона

. (4.27)
Общая ошибка имеет порядок

. (4.28)

Пример 4.3.

Вычислить с точностью 0,01 интеграл

Решение.

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

i	x	y	D y	D2 y	D3 y	D4 y
0	0	0	0,47	- 0,37	0,26	- 0,2
1	0,25	0,47	0,1	- 0,11	0,06
2	0,5	0,57	-0,01	- 0,05
3	0,75	0,56	- 0,06
4	1,0	0,5

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

Ошибка по формуле (4.27)

.

Значение интеграла определяем по формуле (4.26)

4.4. Операторы дифференцирования и интегрирования

Для вычисления производной в системе MATHCAD используется метод Риддера. Первая оценка производной делается по соответствующей разделенной разности. Другие оценки делаются при уменьшенном размере шага. При значении ошибки менее встроенной переменной TOL возвращается значение производной. В противном случае выдается сообщение "отсутствует сходимость".

Алгоритмы дифференцирования являются неустойчивыми по отношению к ошибкам округления или усечения. Особенно сильно это обстоятельство проявляется вблизи точек разрыва или других особенностей функции. Алгоритм Риддера дает 7 или 8 верных значащих цифр первой производной, при условии, что точка, в которой производится дифференцирование, не находится вблизи особенностей функции. Точность этого алгоритма убывает примерно на одну значащую цифру на единицу увеличения порядка производной.

В Маткаде имеются два оператора дифференцирования. Оператор производной первого порядка

и оператор производной высших порядков

.

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

Хотя операция дифференцирование возвращает только одно число, можно определить функцию как производную от другой функции.

f(x) : =

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

Имеется возможность определения производной в символьном виде. Для этого используется обычная запись оператора дифференцирования и команда [Ctrl,?]. Другой способ символьного определения производной - это выделение переменной и использование команды "Дифференцировать по переменной" из меню "Символика".

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

Оператор интегрирования предназначен для численного вычисления определенного интеграла методом Ромберга

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

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

Точность интегрирования методом Ромберга определяется сравнением последовательных значений интеграла при вычислениях с уменьшающимся шагом. Значение интеграла возвращается, когда разность двух последних вычисленных величин меньше встроенной переменной TOL. Если заданная точность не достигается, то выдаётся сообщение об ошибке "не сходится".

Для определения точного значения определенного интеграла или символьного выражения для неопределенного интеграла могут быть использованы возможности символьного процессора. Неопределенный интеграл может быть определен с использованием оператора интегрирования с выражением, выделенным рамкой, и команды
[Shift, F9], или после выделения переменной интегрирования по команде "Интегрировать по переменной" из меню "Символика".

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

Определённый интеграл определяется символьно по команде [Shift, F9] для окруженного рамкой оператора интегрирования с полями ввода пределов интегрирования и интегрируемого выражения. Сначала вычисляется неопределенный интеграл. Затем из значения первообразной функции от верхнего предела вычитается её значение от нижнего предела. Если символьное интегрирование вычислено успешно, а пределы интегрирования - целые числа, дроби или точные константы, подобно p, то процессор выдаёт точное значение интеграла. Если подынтегральное выражение или один из пределов интегрирования содержит десятичную точку, то символьный ответ будет числом, отображаемым с двадцатью значащими цифрами. Этот ответ обычно совпадает с результатом численного интегрирования.