Численное интегрирование и дифференцирование в среде MATLAB. Лабораторная работа 7 по дисциплине Математическое моделирование в среде matlab Численное интегрирование и дифференцирование в среде matlab

• Название: Лабораторная работа 7 по дисциплине Математическое моделирование в среде matlab Численное интегрирование и дифференцирование в среде matlab
• Анкор: Численное интегрирование и дифференцирование в среде MATLAB
• Дата: 26.11.2021
• Размер: 1.01 Mb.
• Имя файла: MTB_7.pdf
• Тип: Лабораторная работа; #282737

Составитель: Афонин В.В.
Лабораторная работа № 7
по дисциплине «Математическое моделирование в среде MATLAB»
Численное интегрирование и дифференцирование
в среде MATLAB
Цель работы: изучить методы численного дифференцирования и численного нахо- ждения определенных интегралов с рассмотрением практических примеров.
ТЕОРЕТИЧЕСКАЯ ЧАСТЬ
Если определенный интеграл содержит подынтегральную функцию (инте- грант) типа табличного вида, то аналитический расчет может быть выполнен по формуле Ньютона – Лейбницы. Задача численного интегрирования возникает в слу- чаях, когда подынтегральная функция не является табличной или крайне сложно подобрать соответствующую подстановку. Кроме того, часто возникает задача ин- тегрирования экспериментальных данных, аналитическое описание которых неиз- вестно.
В среде MATLAB к самой распространенной функции вычисления опреде- ленного интеграла относится trapz, в которой заложен метод трапеций. Геометриче- ский смысл определенного интеграла (рис. 7.1) – площадь криволинейной трапеции, ограниченной осью OX, кривой f(x), и прямыми x = a и x = b.
Рис. 7.1. Пример расчета определенного интеграла с помощью криволинейных трапеций

Составитель: Афонин В.В.
- 2 -
Численные методы интегрирования основаны на различных способах оценки площади, составленной из суммы криволинейных трапеций, поэтому полученные формулы численного интегрирования называются квадратурными (формулами вы- числения площади). В соответствии с рис. 7.1, расчет площади криволинейной тра- пеции может быть в простейшем случае может быть в виде следующих этапов. От- резок [a, b] делят на n необязательно равных частей – элементарных отрезков. При- нято такое деление отрезка называть сеткой, а точки деления по оси абсцисс x
0
,
x
1
,…, x
i
, ..., x
n
– узлами сетки. Если сетка равномерная, то
n
a
b
h


– шаг сетки, при интегрировании – шаг интегрирования, а координата i-го узла вы- числяется по формуле
...,
,
2
,
1
,
0
,
n
i
ih
a
x
i



Полная площадь криволинейной трапеции, соответственно значение опреде- ленного интеграла I состоит из n элементарных криволинейных трапеций – элемен- тарных площадей S
i
:
1 0




n
i
i
S
I
Рассмотрим пример использования библиотечной функции trapz по вычисле- нию полуволны синусоиды, которая задана дискретными значениями для эмуляции экспериментального одномерного массива. Программный код решения примера в
MATLAB: clear,clc x = linspace(0, pi, 500); % 500 точек разбиения интервала интегрирования y = sin(x);
S = trapz(x, y); fprintf('\n\tI = S = %0.6f\n', S)
Результат:
I = S = 1.999993
Проверка с использованием формулы Ньютона – Лейбницы:
2 1
1 1
)
1
(
))
0
cos(
(
))
(cos(
)
cos(
)
sin(
0 0



















x
dx
x
S
I
Как видно, полученные результаты близки друг другу. С увеличением точек разбиения интервала интегрирования результат, получаемый с помощью trapz, будет неограниченно приближаться к аналитическому результату, т. е. к 2.

Составитель: Афонин В.В.
- 3 -
Кроме функции trapz в MATLAB имеются еще функции интегрирования оп- ределенных интегралов, включая двойных и тройных: integral, integral2, integral3.
Эти функции в качестве первого аргумента принимают указатель/дескриптор на функцию, которая описывает подынтегральную функцию, или на имя анонимной функции. Рассмотрим тот же пример с полуволной синусоиды: clear, clc f = @(x) sin(x);
S = integral(f, 0, pi); fprintf('\n\tI = S = %0.15f\n', S)
Результат:
I = S = 2.000000000000000
Рассмотрим еще пример с неберущимся определенным интегралом, в котором интегрант представляет собой функцию плотности стандартного нормального зако- на распределения непрерывной случайной величины на всей действительной оси: clear,clc f = @(x) exp(-x.^2/2);
S1 = integral(f, -inf, inf);
S2 = integral(f, -inf, inf, 'RelTol', 1e-12, 'AbsTol', 1e-12); fprintf('\n\tS1 = %0.15f\n\tS2 = %0.15f\n', S1/sqrt(2*pi), S2/sqrt(2*pi))
Результат:
S1 = 1.000000000000038
S2 = 1.000000000000000
В функции integral заложен численный метод Симпсона, при котором на каж- дом отрезке подынтегральную функцию заменяют параболой, проходящей через три точки - через начало отрезка, его середину и последнюю точку отрезка, при этом разбиение всего отрезка интегрирования предполагается на равные части. Метод
Симпсона более точен, чем метод трапеций, реализованный в trapz.
Результат S2 получен для integral с дополнительными аргументами, опреде- ляющие относительную (RelTol) и абсолютную (AbsTol) точности вычисления за- данного интеграла, который имеет следующий вид:






2 1
2
/
2
dx
e
I
x

Результаты интегрирования можно получить с применением символьных пе- ременных, например

Составитель: Афонин В.В.
- 4 - clear,clc syms x f1 = sin(x);
I1 = int(f1, 0, pi); fprintf('\n\tI1 = %0.15f\n', double(I1))
%% Добавление переменной «y» для вычисления двойного интеграла syms y f2 = 4 - x^2 - y^2;
I2 = int(int(f2, y, 0, 3/2), x, 0, 1); fprintf('\n\tI2 = ') disp(I2)
Результаты:
I1 = 2.000000000000000
I2 = 35/8
Вид двойного интеграла следующий:
 



2
/
3 0
1 0
2 2
)
4
(
2
dx
y
x
dy
I
Примечание. При вычислении кратных интегралов следует правильно определить порядок интегрирования, т. е. по каким переменным со своими пределами.
Формулы для численного вычисления интегралов называются квадратурными формулами.
Широко известны следующие формулы: формула прямоугольников (правых, левых, средних прямоугольников), формула трапеций (реализована в trapz, cumtrapz системы MATLAB), формулы Симпсона, Ньютона – Котеса (реализованы в integral, integral2, integral3 системы MATLAB).
Если отрезок интегрирования [a, b] разбить на n частей с помощью точек (уз- лы квадратурной формулы) x
i
, i = 0, 1, 2, ..., n (см. рис. 7.1), то
)
(
)
(
1 1
 





n
i
x
x
b
a
i
i
dx
x
f
dx
x
f
I
Приведем выражения для программной реализации перечисленных формул.
Формула прямоугольников. В ней частичный интеграл приближается сле- дующей площадью прямоугольника:

Составитель: Афонин В.В.
- 5 -


,
,
2
,
1
,
2
/
2
,
,
)
(
1 1
2
/
1 1
2
/
1 1
n
i
h
x
x
x
x
x
x
h
h
x
f
dx
x
f
i
i
i
i
i
i
i
i
x
x
i
i















(7.1)
С учетом (7.1) получается составная формула прямоугольников:


)
(
2
/
)
1
(
2
/
3 2
/
1
h
f
f
f
dx
x
f
n
b
a






(7.2)
Точность формулы (7.2) не превышает О(h
2
). Ясно, чем меньше шаг h, тем точнее будет вычислен заданный интеграл.
Формула трапеций. Для нее частичный интеграл приближается следующей формулой трапецией:

  
,
2
)
(
2
/
1 1
h
x
f
x
f
dx
x
f
i
i
x
x
i
i





(7.3) где f
(x
i–1/2
), f
(x
i
) – стороны прямоугольной трапеции, h – высота трапеции.
Суммируя частичные интегралы вида (7.3), получим составную формулу тра- пеций:
)
(
2
)
(
)
(
)
(
1 1
0
h
x
f
x
f
x
f
dx
x
f
n
i
i
n
b
a











(7.4)
Формула Симпсона. Ее получают, когда на частичном интервале интегриро- вания подынтегральную функцию интерполируют многочленом Лагранжа второй степени через узловые точки (x
i–1
, f
i–1
), (x
i–1/2
, f
i–1/2
), (x
i
, f
i
):


).
(
,
4 6
)
(
2
/
1 1
1
i
i
i
i
i
x
x
x
f
f
f
f
f
h
dx
x
f
i
i








(7.5)
Применяя (7.5) ко всем частичным интервалам, получаем составную формулу
Симпсона в виде

Составитель: Афонин В.В.
- 6 -


)].
(
4 2
[
6 4
6
)
(
2
/
)
1
(
2
/
3 2
/
1 1
1 0
1 2
/
1 1




















n
n
i
i
n
n
i
i
i
i
b
a
f
f
f
f
f
f
h
f
f
f
h
dx
x
f
(7.6)
Чтобы избавиться от дробных индексов в формуле (7.6) следует положить
x
i
= a + 0,5hi, i = 0, 1, ... , 2n, hn = b – a. Тогда составная формула Симпсона приобре- тает следующий вид:
)].
(
4
)
(
2
[
6
)
(
1 2
3 1
2 2
4 2
2 0














n
n
n
b
a
f
f
f
f
f
f
f
f
n
a
b
dx
x
f
(7.7)
В руководствах по численным методам показано, что скорость уменьшения погрешности при уменьшении шага h у формулы Симпсона (7.7) или (7.7) будет больше, чем у формулы прямоугольников или формулы трапеций.
Численное дифференцирование, например, функции одной переменной осно- вывается на определении производной, т. е.
)
(
)
(
lim lim
)
(
0 0
x
x
f
x
x
f
x
y
dx
x
df
dx
dy
x
x













(7.8)
При конечных значениях приращений аппроксимация производной выражает- ся с помощью конечных разностей, когда функция f (x) задана в табличном виде
y
0
, y
1
, ..., y
n
в узлах x
0
, x
1
, ..., x
n
. Предполагается, чтошаг h между узлами постоянный.
В зависимости от способа вычисления конечных разностей используют следующие формулы вычисления производной в одной и той же точке:
формула левых разностей
;
,
,
0 1
1 0
1 1
h
y
y
dx
dy
h
x
y
y
y







(7.9)
формула правых разностей
;
,
,
1 2
1 1
2 1
h
y
y
dx
dy
h
x
y
y
y







(7.10)
формула центральных разностей
2
,
,
0 2
1 0
2 1
h
y
y
dx
dy
h
x
y
y
y







(7.11)

Составитель: Афонин В.В.
- 7 -
Формула центральной разности (7.11) имеет погрешность (погрешность усе- чения) порядка О(h
2
). Существует формула численного дифференцирования цен- тральной разности порядка О(h
4
), для которой должны выполняться условия при- надлежности узлов: x – 2h, x – h, x, x + h, x + 2h

[a, b]. Эта формула имеет следую- щий вид:
12
)
2
(
)
(
8
)
(
8
)
2
(
)
(
h
h
x
f
h
x
f
h
x
f
h
x
f
dx
x
df









(7.12)
Используя формулы (7.9), (7.10), (7.11), можно определить старшие производ- ные. Например, для 2-производной в точке получим
2
/
)
(
/
)
(
/
/
2 0
1 2
0 1
1 2
1 2
1 2
1 2
h
y
y
y
h
h
y
y
h
y
y
h
dx
dy
dx
dy
dx
dy
dx
d
dx
y
d

















(7.12)
Следует заметить, что в целом формулы численного дифференцирования с определенной точностью нельзя применять к некоторым функциям, например, к вы- соко осциллирующимся и ряда других. Поэтому формулы численного дифференци- рования применяются с осторожностью с предварительным анализом, как аналити- ческих функций, так и тех, которые заданы таблично.
ПРАКТИЧЕСКАЯ ЧАСТЬ
Пример 1. Численно продифференцировать синусоидальную функцию и построить графики.
Программный код решения примера clear, clc, close x = linspace(0, 2*pi, 50); %% 50 узловых точек h = x(2)-x(1); %% постоянный шаг между узлами y = sin(x); dy = diff(y)/h;
%% whos y dy
F = figure;
F.Name = 'Пример численного дифференцирования';
F.Color = 'w'; grid on line(x, y, 'LineWidth', 1.25) line(x(1:end-1), dy, 'color', 'r','LineWidth', 2) xlabel('\it\fontsize{12}\fontname{Times}x') legend('\fontsize{11}\fontname{Times}sin(x)',...
'\fontsize{11}\fontname{Times}d(sin(x))/dx', 'Location', 'Best')

Составитель: Афонин В.В.
- 8 -
Пример выполнения программы показан на рис. 7.2.
Рис. 7.2. Синусоида и ее производная
Задание 1 1. Определите значение производной функции sin(x) в точке Z/(20

), где Z –
номер Вашей фамилии в алфавите. Сравните значение функции косинуса в этой же точке, т. е. cos(Z/(20

)).
2. В программе вместо табличной функции diff примените формулы левых, пра- вых, центральных разностей. Определите шаг h, при котором среднеквадрати- ческая погрешность между значениями производной, полученной по форму- лам конечных разностей по сравнению с diff не превышала бы 10е–3.
Пример 2. Численно и аналитически вычислить объем тела, полученного вра- щением фигуры Ф вокруг указанной оси координат:
y = 1/(1 + x2), y = x/2, x = 0; Oy.
Построить фигуру вращения.
Примечание. Для вычисления объема тела
y
V
, полученного вращением фигуры Ф, следует воспользоваться формулой
,
)
(
2 1
2



b
a
y
dx
y
y
x
V

где
1
,
0
,
2
/
),
1
/(
1 1
2 2





b
a
x
y
x
y
Пример построения фигуры вращения Ф показан на рис. 7.3.

Составитель: Афонин В.В.
- 9 -
Рис. 7.3. Пример построения фигуры вращения фигуры Ф: y = 1/(1 + x2), y = x/2, x = 0; Oy.
Пример 3. Численно и аналитически определить объем тела, ограниченного плоскостями y = 0, z = 0, z = 3 и поверхностью
2 2
x
x
y


Схематичное изображение тела, объем которого следует определить, показано на рис. 7.4.
Рис. 7.4. Схема тела для определения его объема

Составитель: Афонин В.В.
- 10 -
Индивидуальные задания. Численно заданным методом/формулой вычис- лить тестовые определенные интегралы и построить их подынтегральные функции в соответствии с приводимым ниже списком. Если во втором столбце списка указан не номер формулы, то это означает имя библиотечной функции MATLAB.
Примечание. Приветствуются аналитические решения – для сравнения.
Фамилия Имя Отчество
Метод интегрирования
Тестовый интеграл
Абрамов Дмитрий Александрович
Формула (7.7).


2 1
4 2
1
dx
x
x
Абылказымов Эмир Тагаевич
Применить trapz.


3
/
1 0
3 2
1
x
dx
Алкхасонех Мохаммад Моиен Мохаммад
Применить integral.


1 0
4 1 x
dx
Афанасьев Артѐм Александрович
Применить integral.
dx
x
x



1 0
6 4
1
)
1
(
Баранов Олег Алексеевич
Применить trapz.
dx
e
x


)
2
ln(
0 1
Борьков Данил Игоревич
Формула (7.4).
dx
x
x



1 0
2 1
)
1
ln(
Горбункова Екатерина Дмитриевна
Применить integral.


3 2
)
ln(
1
x
dx
Гуляев Дмитрий Игоревич
Применить trapz.
dx
x

4
/
0
)]
ln[cos(

Гурин Олег Игоревич
Применить integral.
dx
x
x


100 10
)
1
ln(
Денисова Валерия Андреевна
Применить trapz.
dx
x


2
/
0
)
sin(
6
,
0 1


Составитель: Афонин В.В.
- 11 -
Дитяткин Алексей Андреевич
Применить int.
dx
x
x
x



0 2
)
(
cos
1
)
sin(
Дотолев Кирилл Витальевич
Применить integral.
dx
x
x
x


1
,
1 1
,
0 2
)
sin(
3
Забиров Иван Евгеньевич
Формула (7.2).


6
,
1 8
,
0 2
1 2x
dx
Иноземцев Семен Александрович
Применить trapz.


7
,
2 2
,
1 2
2
,
3
x
dx
Казакова Дарья Андреевна
Применить integral.


2 1
2 3
,
1 2x
dx
Карасев Илья Сергеевич
Применить trapz.


2
,
1 2
,
0 2
1
x
dx
Качалкин Степан Сергеевич
Применить integral.


4
,
1 8
,
0 2
3 2
dx
x
dx
Киселев Антон Алексеевич
Применить trapz.


2
,
1 4
,
0 2
5
,
0 2
x
dx
Колядина Елизавета Валерьевна
Применить integral.


1
,
2 4
,
1 2
1 3x
dx
Масляев Дмитрий Сергеевич
Формула (7.4).


2 0
2 2
2
)
1
(
dx
x
x
Мишин Павел Александрович
Формула (7.12).
dx
x
x

2
/
0
)
sin(

Немчинова Полина Андреевна
Применить integral
dx
x
x
x


1 0
2 3
1
)
arcsin(

Составитель: Афонин В.В.
- 12 -
Нестеренко Макар Денисович
Применить integral.


4
,
2 2
,
1 2
5
,
0
x
dx
Огрин Егор Викторович
Применить trapz.
dx
x


3 0
2
)
(
cos
1
Пешев Данила Сергеевич
Применить integral.
dx
x

5
,
0 0
)
(
tg
Правосудов Андрей Романович
Формула (7.12)

2
,
2 2
,
1 2
)
ln(
6
,
2
dx
x
x
Романов Евгений Алексеевич
Применить trapz.
dx
x
x

7 1
)
exp(
Рубцов Игорь Александрович
Применить integral.

10 4
)
lg( x
dx
Самарина Анастасия Сергеевна
Применить integral.
dx
x
x


4 2
2 4
100
Сергеев Илья Сергеевич
Формула (7.4)
dx
x




0
)
25
,
2
)
cos(
3 1
ln(
Степанова Инна Сергеевна
Формула (7.12)
dx
x


8 1
Томилина Ангелина Александровна
Применить integral.


3 0
)
(
arctg
dx
x
x
Ухрѐнков Юрий Борисович
Формула (7.2)



75
,
0 0
2 1
)
1
(
x
x
dx