Мур. Тема 2.4 Диф Интегр. Тема Численное вычисление производных и интегралов 4 Постановка задачи численных вычислений

• Название: Тема Численное вычисление производных и интегралов 4 Постановка задачи численных вычислений
• Дата: 16.03.2023
• Размер: 1.16 Mb.
• Имя файла: Тема 2.4 Диф Интегр.pdf
• Тип: Документы; #994263

Тема 2.4. Численное вычисление
производных и интегралов
2.4.1. Постановка задачи численных вычислений
производных
и конечных разностей
Производной функции y=f(x) называется одна из функций:
dy
df(x)
y ;
; f (x);
dx
dx








, которая, при каждом значении независимой переменной х, равна пределу отношения приращения функции к приращению независи- мой переменной х при стремлении Δх к нулю.
Процесс вычисления производной называется дифференциро-
ванием функции. Дифференцирование происходит в соответствии с математическими правилами дифференцирования.
Производная характеризуется порядком, то есть производная от первой производной называется производной второго порядка, от второй – третьего порядка и так далее.
Однако такое вычисление производных возможно только то- гда, когда функция задана аналитически и при этом достаточно про- ста. Именно поэтому в случае, если функция задана таблицей значе- ний, для нахождения производной функции в каком-либо узле задан- ной таблицы обычно используют интерполяционные многочлены невысоких степеней, которые хорошо приближают функцию к окрестности значений этого узла. Для этого обычно используют вы- ражения, основанные на таблице конечных разностей. Наиболее из- вестные из них – формулы Ньютона, применение которой возможно только, если узлы таблицы являются равноотстоящими.

Чтобы вычислить производные с использованием таблицы ко- нечных разностей, предположим, что решается следующая задача: при заданной таблице значений функции y
i
=f(x
i
), где i=0,1,…,n, где шаг изменения x
i
является постоянным (
h=x
i
-x
i-1
=const), тре- буется определить таблицу значений производных в точках
x
i
С использованием значений конечных разностей производные функции в точках x
i
можно определить, как:






i
i
i
i
y
y
y (x )
.
x
h
Где значения производных вычисляются с использованием пра- вых конечных разностей.
Очевидно, что приведенная выше формула позволяет опреде- лить значения производных во всех точках, кроме x
n
. Вычислить производную в точке x
n
можно по аналогичной формуле, в которой используются левые конечные разности:
i
i 1
i
y
y
y (x )
, i
1, 2,..., n.
h





Поскольку использование этих выражений требует вычисления конечных разностей, то точность вычисления производной напря- мую зависит от величины шага между узловыми точками таблицы.
Следует отметить, что по величине конечных разностей можно сделать вывод о степени интерполяционного полинома, описываю- щего таблично заданную функцию. Если для таблицы с равноотсто- ящими узлами конечные разности k
-го порядка постоянны или соиз- меримы с заданной погрешностью, то функцию можно представить многочленом k-й степени.
Рассмотрим, например, таблицу конечных разностей для мно- гочлена y=x
2
-3x+2.

Таблица 2.4.1 – Таблица конечных разностей
x
y

y

2
y

3
y
1.0 0.0
-0.16 0.08 0
1.2
-0.16
-0.08 0.08 0
1.4
-0.24 0
0.08 1.6
-0.24 0.08 1.8
-0.16
В данном примере все конечные разности второго порядка равны 0.08. Это говорит о том, что функцию, заданную таблично, можно представить многочленом второй степени.
Введя понятие конечных разностей, рассмотрим еще две формы записей интерполяционных полиномов.
Рассмотрим имеющиеся в системе Scilab средства дифференци- рования функций, заданных аналитически и таблично, а также сред- ства, позволяющие проводить вычисления конечных разностей.
2.4.2. Вычисление производных средствами Scilab
Вычисление производной от аналитической функции
Версия системы Scilab 6.1.1, к сожалению, не позволяет полу- чить аналитических выражений производных от функции f(x). Од- нако в ней имеется функция Scilab numderivative, предназначен- ная для вычисления значений производных 1-го или 2-го порядка в заданных точках от аналитических функций. Эта функция имеет сле- дующие форматы:
numderivative(f,x)
J=numderivative(f,x)
J=numderivative(f,x,h)
h=numderivative(f,x,h,order)
[J,H]=numderivative(f,x,h,order).
Входные параметры:

f – функция или список функций, от которых берутся произ- водные;
х – вектор значений независимых переменных, относительно которых вычисляются производные;
h – шаг, используемый в конечно-разностных приближениях.
Если h не указан, то шаг по умолчанию вычисляется в зависи- мости от х и порядка. Если h-матрица 1х1, то она расширяется до того же размера, что и х;
order – порядок конечных разностей (порядок по умолчанию равен 2), доступные значения order: 1,2.
Выходные параметры:
J – матрица (необязательный параметр), для функции от од- ной переменной размером 1х1, для функции от двух переменных –
2х2 (значения частных производных выводятся на ее диагонали);
H – матрица (необязательный параметр), в которую выводятся частные производные второго порядка.
--> // Вычисление производной от функции
--> // f(x) в точке х=3
--> // Описание исходной f(x)
--> deff('[y] = f(x)','y = 2^x – 4*x');
-->
--> numderivative(f, 3) // Вычисление производной
--> // функции f в точке ans =
1.5451774
-->
--> // Вычисление производной f(x),
--> // заданной аналитическим выражением
--> deff('[y1] = f1(x)','y1 = 2^x*log(2) - 4');
-->
--> f1(3) // Значение в точке х=3 ans =
1.5451774
Рисунок 2.4.2 – Вычисление производной в точке функцией
numderivative

На рисунке 2.4.2 приведены два варианта вычисления значения производной в точке х=3 от функций, описанных с использованием оператора deff. Причем, в первом примере значение производной вычислено с помощью функции numderivative, а во втором – для сравнения с использованием выражения производной, полученного заранее в явном виде. Как показали приведенные вычисления, чис- ловые значения производных в точке х=3, полученные двумя спосо- бами, полностью совпали.
На рис. 2.4.3 приведены примеры получения выражений произ- водной от полинома, созданного в аналитическом виде (использова- ние функции derivat только для полиномов) и значения производ- ных в точках x=[1 3 6] (использование функции numderiva-
tive).

--> // Вычисление производных полинома
--> p = poly([1,2,3], "x","c") //Определение полинома
p =
1 +2x +3x
2
--> p1 = derivat(p) // Определение производной от p
p1 =
2 +6x
--> ps = pol2str(p) // Преобразование полинома
--> // в символьную строку
ps =
1+2*x+3*x^2
-->
--> deff('[sd] = fd(x)','sd = 1 + 2*x + 3*x^2')
--> x = [1 3 6];
--> numderivative(fd, x) ans =
8. 0. 0.
0. 20. 0.
0. 0. 38.
-->
-->// Описание функции fp с ps
--> deff('sp = fp(x)', 'sp='+ ps)
--> numderivative(fp, x) ans =
8. 0. 0.
0. 20. 0.
0. 0. 38.
Рисунке 2.4.3 – Вычисление производных от полинома
Функция numderivative позволяет вычислить и значения частных производных. На рисунке 2.4.4 приведен пример вычисле- ния частных производных от функции трех переменных
f(x
1
,x
2
,x
3
)=x
1
2
x
2
+x
2
2
x
3
+x
3
2
x
1
при х
1
=1, x
2
=3 и х
3
=6.

--> // Вычисление частных производных функции
--> // трех переменных
-->
--> deff('s = f(x)', 's = x(1)^2*x(2) + x(2)^2*x(3)
+ x(3)^2*x(1)')
-->
--> x = [1 3 6]; // Вектор значений переменных
--> numderivative(f, x, 2) // Вычисление второй
--> //производной ans =
42. 37. 21.
Рис.2.4.4 – Вычисление значений частных производных от функции трех переменных с использованием функции numderivative
Построить касательную к функции f(x)=
(
3*x
2
-7
)
/
(
2*x+1
)
(рисунок 2.4.5)

Рисунок 2.3.4 – Построение уравнения и графика производной функции
Вычисление производной от табличной функции
Если функция задана таблично, то для вычисления значения производной в точке, необходимо предварительно получить значе- ния конечных разностей, поскольку в основе численного дифферен- цирования, как уже отмечалось выше, лежит аппроксимация функ- ции, от которой берется производная, заменяемая интерполяцион- ным многочленом Ньютона
2
3
0
0
0
y
( y
y / 2 !
y / 3!
...) / h
  
 
 

.???????????????
Именно поэтому численное дифференцирование еще называют
аппроксимированным дифференцированием.

Для получения конечных разностей, входящих в состав интер- поляционного многочлена, в Scilab используется функция diff, ко- торая может иметь один из следующих форматов:
dy=diff(y);
dy=diff(y,n);
dy=diff(y,n,dim),
где: y – вектор или матрица (вещественная, комплексная или поли- номиальная) табличных значений функции;
n– порядок конечных разностей (целое число), по умолчанию значение n равно 1;
dim– размерность матрицы, по которой происходит вычисле- ние конечных разностей. Размерность может принимать значения как 'r', 'c' или 1, 2 соответственно (по строкам или столбцам), так и'*' (по всем элементам матицы). По умолчанию значение dim
равно "*", то есть запись diff(y,n) равносильна записи
diff(y,n,'*').Если параметр dim имеет значение 'r', то это эк- вивалентно записи dim=1, а dim='c' – эквивалентно dim=2.
dy – скаляр или матрица значений конечных разностей.
Формат dy=diff(y) предназначен для вычисления разностей между соседними значениями y (y(2:$)-y(1:$-1)), поскольку по умолчанию n=1. При использовании формата dy=diff(y,n)
функция diff вычисляет n-конечные разности, а использование формата dy=diff(y,n,dim) позволяет вычислить n-ые разности при заданной размерности dim.
На рис.2.4.6 приведен пример вычисления значений производ- ной от функции y(x)=log10(x), заданной в виде таблицы. Зная аналитическое выражение для производной от функции y(x)
(y'=1/ln(10)/x), в примере проведена проверка результата, по- лученного при использовании функции diff :

--> // Вычисление значений производных для функции,
--> // заданной таблично
--> x = 10 : 5 : 25; h = 5;
--> y = log10(x); // Дифференцируемая функция
-->
--> // Конечные разности 1-го, 2-го и 3-го порядка
--> dy = diff(y); dy2 = diff(y,2); dy3 = diff(y,3);
-->
--> // Приближенные y'(x) в точках x = 10, 15, 20
--> y10 = (dy(1) - dy2(1)/2 + dy(1)/3)/h y10 =
0.0520729
--> Y15 = (dy(2) - dy2(2)/2) / h
Y15 =
0.0277906
--> Y20 = dy(3) / h
Y20 =
0.019382
--> // Значения производных функции log10(x)
--> (1/log(10)) ./ x(1:$-1) ans =
0.0434294 0.028953 0.0217147
Рисунок 2.4.6. Вычисление производных по таблице значений функции

На рисунке. 2.4.7 показана еще одна возможность приближен- ного вычисления производных функции y=sin(x), заданной таб- лично. Принимая во внимание, что для исходной функции известно аналитическое выражение производной (y'=cos(x)), с помощью функции norm (норма матрицы), можно получить значение макси- мального отклонения приближенных значений производных от их соответствующих точных значений, полученных по аналитической формуле.
--> // Аппроксимированное дифференцирование
--> x = 1 : 0.01 : 1.05; h = 0.01;
--> // Функция, от которой берется производная
--> y = sin(x);
--> // Аппроксимированное вычисление производной
--> dy = diff(sin(x)./h) dy =
0.536086 0.5276177 0.5190967 0.5105238
--> y1 = cos(x) // Производная от функции sin(x)
y1 =
0.5403023 0.5318607 0.523366 0.5148188 0.5062203
--> // Максимальный элемент вектора модулей разно-
стей
--> norm(dy - cos(x(1:$-1)), %inf) ans =
0.0043204
Рисунок 2.4.7 – Аппроксимированное дифференцирование
Поскольку результат выполнения функции norm(V,%inf) ра- вен значению максимального элемента вектора V по модулю, то в данном случае он равен значению 0.0043204.

2.4.3. Постановка задачи численного
интегрирования
Интегрирование – это распространенная математическая опе- рация для решения многих задач науки и техники. С его помощью проводят вычисление площадей, скорости, ускорения, перемещения и многое другое. Интегрирование простых функций можно легко провести аналитически, но чем сложнее подынтегральная функция, тем труднее, а иногда невозможно, получить точное значение инте- грала. Следует также отметить, что, если в курсах математического анализа подынтегральная функция обычно задана аналитически, то в прикладных задачах она может быть задана в виде таблицы, по- скольку здесь часто используются данные экспериментов.
Из курса математического анализа известно, что, если функция
f(x) непрерывна и дифференцируема на отрезке [a;b], то опреде- ленный интеграл от этой функции в пределах от a до b существует и может быть вычислен по формуле Ньютона-Лейбница:
b
a
s
f(x)dx,


s =
b a
f x dx
F b
F a где f x
F x
( )
( )
( ),
( )
'( ).




где: f(x) – подынтегральная функция;
а и b – пределы интегрирования.
Исходя из геометрической интерпретации, значение интеграла
s численно равна площади области, ограниченной графиком функции, осью абсцисс с ОX и двумя прямыми x=а и x=b.
Если первообразную функцию F(x) не удается выразить ана- литически через элементарные функции или если при проведении практических расчетов подынтегральная функция f(x) задается в виде таблицы, то это приводит к необходимости замены аналитиче- ского интегрирования численными методами.
Для функции f(x), заданной в прямоугольной системе коорди- нат на интервале [a;b], этот интеграл численно равен площади,

ограниченной кривой f(x), осью Ox и двумя ординатами ac и bd
(рисунок 2.4.8).
При численном интегрировании заданная область разбивается на несколько частей, вычисляется площадь каждой маленькой части, а затем они суммируются.
Рисунок 2.4.8 –
Задача численного интегрирования заключается в нахождении значения определенного интеграла через ряд значений подынте- гральной функции y
i
=f(x
i
), заданной в точках x
i
(i=0,1,…,n).
Причем, x
0
=a, x
n
=b. Чаще всего интервал разбивают на под интер- валы длиной h=x
i+1
-x
i
Применительно к однократному интегралу, формулы числен- ного интегрирования представляют собой квадратурные формулы вида: b
n i
i i 0
a f(x)dx
A f(x ),




где Ai – числовые коэффициенты, называемые весами квадратур-
ной формулы, а xi – точки отрезка – узлами квадратурной фор-
мулы, n – целое положительное число.
Известно, что искомый определенный интеграл можно пред- ставить в виде суммы интегралов:

n
i i 0
I
I



На каждом i-м отрезке функция интерполируется
(заменяется) некоторой другой легко интегрируемой функцией gi(x). В резуль- тате получаем следующую квадратурную формулу:






i i 1
x b
n i 1
a x
f(x)dx g(x)dx
Для решения поставленной задачи подынтегральную функцию
f(x) необходимо заменить приближенной функцией, которая мо- жет быть проинтегрирована в аналитическом виде. В качестве такой функции обычно используют полином Рn(х) с узлами интерполя- ции в точках
х0, х1…,хn. В этих точках значения функции и интер- поляционного полинома полностью совпадают f(xi)=Рn(xi).
Существует множество численных методов, которые отлича- ются способами разбиения области интегрирования на части и спо-
собами вычисления площади каждой из этих частей.
Суть всех формул численного интегрирования состоит в том, что на элементарных отрезках интегрирования подынтегральную функцию заменяют простейшим интерполяционным полиномом, ко- торый легко может быть проинтегрирован в аналитическом виде.
Так, например, для получения формул прямоугольников, трапеций
и Симпсона используют полиномы соответственно нулевой, первой и второй степени.
Для получения простых формул интегрирования используют полиномы нулевой, первой и второй степени и соответственно полу- чают формулы численного интегрирования: прямоугольников,
трапеций и Симпсона.
Очевидно, что замена функции f(x) интерполирующим полино- мом приводит к образованию погрешности вычисления значения ин- теграла b
1 1
a
I
g(x)dx R;
I
I R,


 


где I
1
– точное значение интеграла,
I – значение интеграла, вы- численного численным методом, а
 
1
R | I
I |
– погрешность метода.
Отметим, что увеличение числа под интервалами n (или умень- шение значения длины шага интегрирования h) ведет к уменьшению погрешности.
2.4.4 Методы и алгоритмы численного интегриро-
вания
Метод прямоугольников
Заменим подынтегральную функцию f(x) в пределах элемен- тарного отрезка [xi;xi+1] интерполяционным многочленом нуле- вой степени (рисунок 2.4.9), то есть постоянной величиной, равной либо f(xi), либо f(xi+1).
Рисунок 2.4-9 – Вычисление интегралов методом прямоугольников
Значение элементарного интеграла равно площади прямоуголь- ника, в первом случае I = h∙f(x i
), а во втором I=∙f(xi+1), где h=xi+1 - xi. Для определения значения интеграла на отрезке
[a;b] найдем суммы элементарных интегралов, взяв в первом слу- чае в качестве f(x) – значение подынтегральной функции в левом конце i-го отрезка, а во втором – в правом конце отрезка: x
d c
h x
i x
i+1
P
0
(x i
)
P
0
(x i+1
)
f(x)
y h/2

n 1
i i 0
I
h f(x ),





n i
i 1
I
h f(x ).




Первая Формула называется формулой левых прямоугольни-
ков, а вторая формула – формулой правых прямоугольников.
Для вычисления определенного интеграла может быть исполь- зована и формула средних прямоугольников, в которой на элемен- тарном отрезке интегрирования функция f(x) тоже заменяется ин- терполяционным многочленом нулевой степени, но равным значе- нию функции в середине отрезка: n 1
i 0
h
I
h f(a i h).
2




  

На рисунке 2.4.10 представлен алгоритм метода средних пря- моугольников, в котором оценка погрешности проводится с исполь- зованием метода двойного просчета.

Рисунок 2.4.10 – Алгоритм метода средних прямоугольников
Метод трапеций
Разобьем интервал интегрирования [a;b] на n равных отрезков
(рисунок 2.4.11) и восстановим из полученных точек a,х
1
,x
2
,…,b перпендикуляры до пересечения с графиком функции. Соединив по- следовательно точки пересечения, представим площадь полученной криволинейной трапеции как сумму прямолинейных трапеций, пло- щади которых легко подсчитать. Заменив подынтегральную функ- цию f(x) в пределах элементарного отрезка [x
i
;x
i+1
] интерполя- ционным многочленом первой степени, получим следующие фор- мулы для элементарных площадей:

Рисунок 2.4.11 – Иллюстрация метода трапеции

0 1
0 1
0 1
0 0
1 2
3 1
2 0
1 2
n 2
n 1
n 1
n n 2
n 1
y y
b a
S
(x x ), (x x )
x h
,
2
n y
y y
y b
a b
a y
y b
a
S
,
S
,
S
,
n
2
n
2
n
2
b a
y y
b a
y y
S
S
n
2
n
2












 


















































Тогда общая площадь равна:
0 1
2
n
0 1
n
(y
2y
2y
... y )
b a
I
S
S
... S
n
2


 



 


Отсюда получаем формулу трапеций: b
n 1 0
n i
i i
i 1
a b a
I
f(x)dx
(y y
2
y ),
где y f(x ).
2n










Алгоритм метода трапеций, представленный на рисунке
2.4.11, должен быть дополнен процедурой-функцией f(x), в кото- рой вычисляется значение подынтегральной функции.
Алгоритм метода основан на замене подынтегральной функции f(x) в пределах элементарного отрезка [x
i
;x
i+1
] ин- терполяционным многочленом первой степени. x
d c
a=x
0
b=x n
f(x)
X
i+1
P
1
(x)
x
1
x
2
X
n-1
x i
f(x)

Оценка погрешности проводится с использованием метода двойного просчета.
Рисунок 2.4.11 – Алгоритм метода трапеций
Алгоритм метода основан на замене подынтегральной функ- ции f(x) в пределах элементарного отрезка [x
i
;x
i+1
] интерполя- ционным многочленом первой степени.
Оценка погрешности проводится с использованием метода двойного просчета.

Метод Симпсона
Для получения формулы Симпсона применяется квадратичный интерполирующий полином, следовательно, за элементарный интер- вал интегрирования принимается отрезок [xi;xi+2]. Поэтому разобьем интервал интегрирования [a;b] наn отрезков, где n=2m – четное число (рисунок 2.4.12).
Рисунок 2.4.12 – Иллюстрация метода Симпсона
Для получения интерполирующей функции на интервале
[xi;xi+2] воспользуемся первой интерполяционной формулой
Ньютона, используя в качестве узлов интерполяции точки x
i
, х
i+1
и x
i+2
2
i i
i i
i i 1
i i 2 2
y y
f(x)
y
(x x )
(x x )(x x ),
x
[x ,x
].
1!h
2!h




 





В пределах отрезка [x i
;x i+2
], на котором подынтегральная функ- ция аппроксимирована многочленом, получим приближенную фор- мулу Симпсона: i 2
i 2
i i
x x
2
i i
i i
i i 1
i i 1
i 2 2
x x
y y
h f(x)dx
[y
(x x )
(x x )(x x )]dx
(y
4y y
).
1!h
2!h
3


















x y
f(x)
x i
X
i+1
X
i+2
P
2
(x)
y i+2
y i+1
y i

Для отрезка [x
0
; x
2
]
2 0
x
0 1
2
x h
f(x)dx
(y
4y y ).
3




Для отрезка [x
2
; x
4
]
4 2
x
2 3
4
x h
f(x)dx
(y
4y y ).
3




Тогда для всего интервала интегрирования [a;b] формула
Симпсона выглядит следующим образом:
0 1
2 2
3 4
2m 2 2m 1 2m
0 2m
1 3
2m 1 2
4 2m 2
h
I
[(y
4y y ) (y
4y y ) ... (y
4y y
)]
3
h
[y y
4(y y
y
) 2(y y
y
)],
3










 







 


 
или m
m
0 2m
2i 1 2i 2
i 1
i 2
h
I
(y y
4
y
2
y
),
3










при i
i y
f(x ).

Схема алгоритма метода Симпсона, представленная на ри- сунок 2.4.13, должна быть дополнена процедурой-функцией
f(x), в которой вычисляется значение подынтегральной функ- ции.

Рисунок 2.4.13 – Алгоритм метода Симпсона
Алгоритм метода Симпсона основан на замене подынте- гральной функции f(x) в пределах двух элементарных отрезков
[x
i
;x
i
+2] интерполяционным многочленом второй степени.
Количество отрезков должно быть четным (n=2m).

Оценка погрешности численного интегрирования
Замена подынтегральной функции интерполяционным полино- мом приводит к погрешности вычисления его значения R=|I
1
–I|, где: b
1
a
I
f(x)dx,
I
g(x)dx.




Очевидно, что вычислить эту погрешность можно только, если известно точное значение интеграла. Поэтому на практике принято проводить оценку погрешности численного интегрирования следу- ющим образом (подынтегральная функция задана таблично (Т) или аналитически (А)):

при использовании формул левых или правых прямо-
угольников т
А
л.пр п.пр b a b a
R
| y |,
R
h max | f '(x) |,
где x [a;b],
2 2








при использовании формулы средних прямоугольников т
2
А
2
ср ср b
a b
a
R
|
y |,
R
h max | f ''(x) |,
где x [a;b],
24 24








при использовании формулы трапеций т
2
А
2
тр ср b
a b
a
R
|
y |,
R
h max | f ''(x) |,
где x [a;b],
12 12








при численном интегрировании по формуле Симпсона: т
4
А
4
(4)
Сим
Сим b
a b
a
R
|
y |,
R
h max | f
(x) |,
где x [a;b].
180 180








В приведенных выше формулах: a, b – границы интервала ин- тегрирования; h=(b-a)/n – шаг интегрирования;

y
.
2
y

и
4
Δ у – среднее арифметическое, соответственно, первых, вторых и четвер- тых конечных разностей.
Поскольку в формуле погрешности для метода трапеций при- сутствует вторая производная, а в формуле Симпсона – четвертая, то формула трапеций точна только для линейных функций, а формула
Симпсона для линейных, квадратичных и кубических.
Из приведенных формул видно, что уменьшение шага интегри- рования (h) приводит к уменьшению погрешности. При этом, по- скольку квадратичная интерполяция представляет функцию с боль- шей точностью, чем линейная, то при использовании формулы
Симпсона требуемая точность достигается при меньших значениях n
(количестве разбиений), чем, например, при использовании фор- мулы трапеций и формулы прямоугольников.
Формулы для оценки погрешности могут быть также использо- ваны для выбора числа разбиений n или шага интегрирования h, не- обходимых для обеспечения заданной точности. Однако практиче- ское использование этих формул ограничено в связи с трудоемко- стью их вычислений, поэтому при реализации численных методов на
ПК используется прием, позволяющий получить оценку погрешно- сти в неявном виде.
Этот прием основан на двукратном вычислении значения инте- грала вначале с шагом h (где h=(b-a)/n), а затем с шагом h/2. Полу- ченные значения интегралов I
h и I
h/2
могут быть применены для оценки погрешности интегрирования по формуле: h
h / 2
k
I
I
R,
2 1



где: k=1–для формул левых и правых прямоугольников; k=2
–для формул трапеции и средних прямоугольников; k=4
–для формулы Симпсона.
Если полученная погрешность не удовлетворяет требуемой точности, то вычисляется значение интеграла при h=h/4 и снова оце- нивается погрешность, и т.д. до тех пор, пока не окажется, что по- грешность стала меньше заданной точности. Это правило называется
правилом Рунге (или правилом двойного просчета).

Пример: Вычислить значение интеграла


2
1
x
0
e
dx.
Предположим, что, подынтегральная функция задана таб- лично:
x
0.0
0.1
0.2
0.3
0.4
0.5
f(x)
1.0
0.9900
5
0.96078
9
0.913913
0.85214
4
0.778801
Воспользуемся формулой правых и левых прямоугольников, считая, что h=0.2, а n=(b-a)/h=5, имеем: л.пр
I
0.2(1 0.960789 0.852144 0.697676 0.527292)
0.807580,






n п.пр i
i 1
I
h f(x ) 0.2 (0.960789 0.852144 0.697676 0.527292 0.367879)
0.681156.










Воспользуемся формулой трапеций:

 





ТР
0.2
I
(1 2(0.960789 0.852144 0.697676 0.527292) 0.367879)
0.744368.
2
Воспользуемся формулы формулой средних прямоугольни-
ков: ср
I
0.2 (0.99005 0.913913 0.77880 0.612626 0.444858)
0.7480496







Воспользуемся формулой Симпсона при m=2∙n=10(2∙5) и шаге h=0.1: сим
0.1
I
(1 4(0.99005 0.913913 0.778801 0.612626 0.444858)
3 2 (0.960789 0.852144 0.697676 0.527292) 0.367879)
0,746822.

 





 





0.6
0.7
0.8
0.9
1.0
0.69767
6
0.61262
6
0.527
292
0.4485
8
0.367879

Оценим погрешности каждого из полученных значений, ис- пользуя известное аналитическое выражение подынтегральной функции f(x):
1 2
(4)
4
M
max f (x)
0.86,
M
max | f (x) | 2,
при x [a,b],
M
max f
(x)
12.









Следовательно,
2 2
л.пр
1
ср
2 2
2 4
5
тр
2
сим
2
(b a)
(b a)
R
h M
0.045,
R
h M
0.33 10 ),
2 24
(b a)
(b a)
R
h M
0.66 10 , R
h M
0.66 10 .
12 180






















Анализируя значения погрешностей, можно с уверенностью сказать, что самый точный результат получен с использованием фор- мулы
Симпсона.

2.4.5. Численное вычисление определенных
интегралов средствами Scilab
Вычисление определенных интегралов – inttrap
В системе Scilab для вычисления определенного интеграла
b
a
y(x) x


по формуле трапеций используется функция
inttrap, име- ющая формат:
s=inttrap(y)
s=inttrap(x,y),
где y – подынтегральная функция y(x), заданная таблично;
x – вектор значений независимой переменной – необязатель-
ный параметр, если он отсутствует, то элементы вектора х прини- мают значения номеров вектора y.
Функция inttrap вычисляет площадь области, ограниченной функцией y(x), которая описана набором точек (x,y). Простей- ший пример использования функции inttrap приведен на рисунке 2.4.14.
Здесь подынтегральная функция задается 11-ю узловыми точ- ками, где х меняется от 1 до 10 с шагом 0.5, а функция y принимает в этих точках те же значения.
--> // Пример использования функции inttrap
-->
--> x = 1 : 0.5 : 10; // Вектор аргументов
--> y = 1 : 0.5 : 10; // Вектор значений
--> // подынтегральной функции
--> s = inttrap(x, y) s =
49.5
Рисунок 2.4.14 – Пример использования функции inttrap
На рисунке 2.4.15 приведены несколько примеров вычисления значений определенных интегралов методом трапеций. В последнем примере для функции интегрирования использован формат

inttrap(y). Поэтому здесь в соответствии с форматом функции
inttrap произошла замена элементов вектора х номерами элемен- тов вектора y, в результате чего получено другое значение инте- грала. Система Scilab не выдала ошибки, поскольку в данном случае была вычислена площадь совершенно другой фигуры.
--> // Вычисления определенных интегралов функцией
--> // Пример1
--> a = 5; b = 13; // Пределы интегрирования
--> x = a : b; y = sqrt(2*x - 1);
--> // Подынтегральная функция
--> inttrap(x, y) // Вычисление интеграла с шагом 1 ans =
32.655571
-->
--> // Пример2
--> h = 0.1; x = a : h : b; y = sqrt(2*x -1);
--> inttrap(x,y) //Вычисление интеграла с шагом 0.1 ans =
32.666556
-->
--> // Пример3
--> inttrap(y)
--> // Элементы вектора x равны номерам элементов y
ans =
326.66556
Рисунок 2.4.15 – Вычисление определенных интегралов функцией inttrap
Вычисление определенных интегралов – integrate
Для вычисления интеграла по формуле Симпсона в Scilab при- меняется функция integrate, имеющая следующий формат:
integrate('f','x',a,b)
integrate('f','x',a,b,,er1)
integrate('f','x',a,b,,er1,er2),

где: f – функция, задающая подынтегральное выражение;
x – независимая переменная;
a,b – пределы интегрирования – действительные числа;
er1 и er2 – необязательные параметры – абсолютная и от- носительная погрешности вычисления интеграла – действительные числа.
Примеры использования функции integrateпоказаны на ри- сунке 2.4.16. Особенность использование integrate состоит в том, что при обращении к этой функции шаг интегрирования не задается, а используются параметр er1 и (или) er2 (Пример1). Если погреш- ность вычисления интеграла отсутствует, то вычисление проводится с погрешностью, установленной по умолчанию (Пример2). По умол- чанию значение er1 равно 1.D-8
, а значение er2 равно 1.D-14.
--> //Вычисление интеграла с заданной погрешностью
--> // Пример1 Вычисление интеграла
--> // с абсолютной погрешность .0001
--> q1 = integrate('(2*x-)^0.5', 'x',5,13, 0.0001) q1 =
32.666667
-->
--> // Пример2 Вычисление интеграла с погрешность
--> // по умолчанию
--> q2 = integrate('(2*x - 1)^0.5', 'x', 5, 13) q2 =
32.666667
-->
--> // Пример3 Вычисление интеграла с погрешностью
--> // по умолчанию, а подынтегральное выражение
--> // задано внутренней функцией
--> deff('y = fint((2 * x - 1)^0.5)', 'x', 5, 13)
--> q3 = integrate('fint', 'x', 5, 13) q3 =
32.666667
Рисунок 2.4.16 –. Вычисление определенных интегралов integrate

Вычисление определенных интегралов функцией intg
Для интегрирования в Scilab имеется универсальная функция
intg:
[s,er]=intg(a,b,f,er1,er2),
где: a, b – пределы интегрирования;
f – имя подынтегральной функции, которая может быть задана с помощью внешней функции или в виде набора дискретных точек и должна быть непрерывной;
er1 и er2 – необязательные параметры – абсолютная и отно- сительная погрешности интегрирования. Если погрешность вычис- ления интеграла отсутствует, то вычисление проводится с погреш- ностью, установленной по умолчанию (значение er1 равно
1.D-8, а значение er2 равно
1.D-14
).
Функция
intg возвращает значение интеграла (s) и оценку абсолютной ошибки вычислений (er).
Внешнюю функцию можно задать функцией deff или
function:

если f – функция, то её определение должно иметь видy =
f(t);

если f – список, то этот список должен быть описан как:
list(f,x1,x2,…).
Рассмотрим несколько примеров (рисунок 2.4.17).

--> // Загрузка сценария РИС23404 и функции intg
--> exec('РИС23404.sce ');
--> disp([" I ", " er"], ...
> [I1, er1], [I2, er2], [I3, er3], ...
> [I4, er4], [I5, er5]);
! I er !
-3.58957861311877350 0.00000000001621547
-3.58957861320974690 0.00001678794315296
-3.58957862082392730 0.0015909785599586
-3.58957862082391800 0.00159097855998436
-3.58957861311831690 0.00000000781318699
Рисунок 2.4.17 – Вычисление определенных интегралов функцией intg
Точность полученного результата зависит от заданной абсо- лютной погрешности. Поэтому очевидно, что в Примере1 и При-
мере5 обеспечена максимальная точность.

Опишем подынтегральную функцию y=cos(x)
2
/(1-x)в функции fi(x) и вычислим значение определенного интеграла (с пределами a=2 и b=6) сначала с помощью функции intg, а затем – функции inttrap, разбив при этом отрезок интегрирования на 5, 10 и 20 частей (рисунок 2.4.18).
--> // Сравнение вычисления интеграла intg и
inttrap по точности
-->
--> // Описание подынтегральной функции
--> deff('[y] = fi(x)', 'y = cos(x).^2 ./ (1-x)')
-->
--> a = 2; b = 6;
--> // Вычисление интеграла с помощью функции intg
--> [s,ir] = intg(a, b, fi) ir =
5.844D-10 s =
-0.8781216
-->
--> // Вычисление интеграла с помощью inttrap
--> //s1–5 точекh1=0.8);s2–10(h2=0.4);s3–0(h3=0.2);
--> h1 = 0.8; x = a:h1:b; y = fi(x); s1 = inttrap(x, y)
--> h2 = 0.4; x = a:h2:b; y = fi(x); s1 = inttrap(x, y)
--> h3 = 0.2; x = a:h3:b; y = fi(x); s1 = inttrap(x, y)
--> z = [h1,s1; h2,s2; h3,s3] z =
0.8 -0.8490697 0.4 -0.8711676 0.2 -0.8764039
Рисунок 2.4.18 – Влияние величины шага на точность интегрирования

Из полученных результатов (рис.2.4.18) следует, что, если зна- чение интеграла, полученное с использованием функции intg
можно принять практически за точное значение (относительная по- грешность 5.844D-10), то с помощью функции inttrap, исполь- зующей таблицу значений подынтегральной функции, мы только приблизились к этому значению. При этом из полученных результа- тов явно следует, что чем меньше значение шага интегрирования, тем точнее результат. А поскольку величина шага влияет на количе- ство используемых узловых точек, то становится очевидным, что увеличение числа, используемых в расчете узловых точек приводит к получению более точного значения интеграла.