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

1.

Численные методы анализа.
ч.5-6.
«Всё опыт, опыт! Опыт – это вздор.
Значенья духа опыт не покроет.
Всё, что узнали до сих пор,
искать не стоило. И знать не стоит.»
Монолог Бакалавра. «Фауст», Гёте

2.

5. Численное дифференцирование и
интегрирование функций
5.1. Постановка вопроса
Найти производные указанных порядков от функции f(x), заданной таблично,
либо имеющей сложное аналитическое выражение.
Данную функцию на интересующем отрезке [a,b] заменяют интерполирующей
функцией P(x) (чаще полиномом) и полагают
df ( x) dP ( x)
.
dx
dx
Если известна погрешность для интерполирующей функции
то погрешность производной выражается формулой
r ( x)
df ( x) dP ( x) dR ( x)
dx
dx
dx
То же самое относится и к производным высших порядков.
R ( x) f ( x) P ( x),

3.

5.2. Приближенное дифференцирование на основе
первой интерполяционной формулы Ньютона
Пусть функция y(x) задана в равноотстоящих точках xi (i=0,1,2,…n) отрезка [a,b]
с помощью значений yi=f(xi).
Заранее должно быть известно о существовании соответствующих производных. Для
нахождения производных функцию y(x) заменим интерполяционным полиномом
Ньютона, построенном для системы узлов x,j (j=0,1,2,…k, k ≤ n).
y ( x) y0 q y0
где
q
q (q 1) 2
q (q 1)(q 2) 3
q (q 1)(q 2)(q 3) 4
y0
y0
y0 ...
2!
3!
4!
x x0
,
h
h xi 1 xi
i 0,1,2,...
В качестве x0 следует брать ближайшее табличное значение аргумента.
3

4.

Перемножая биномы получим
q2 q 2
q 3 3q 2 2q 3
q 4 6q 3 11q 2 6q 4
y ( x) y0 q y0
y0
y0
y0 ...
2
6
24
dy dy dq 1 dy
.
dx dq dx h dq
Учтем
В результате:
dy ( x) 1
2q 1 2
3q 2 6q 2 3
4q 3 18q 2 22q 6 4
y0
y0
y0
y0 ... .
dx
h
2
6
24
Далее, поскольку
получим
d 2 y d dy d dy dq 1 d dy
,
2
dx
dx dx dq dx dx h dq dx
2
d 2 y ( x) 1 2
6
q
18q 11 4
3
2 y0 q 1 y0
y0 ... .
2
dx
h
12
Аналогично можно получить формулы и для производных более высокого порядка.
4

5.

5.3. Приближенное дифференцирование для равноотстоящих
точек (узлов), выраженных через значения функций в этих точках
на основе интерполяционной формулы Лагранжа
Для данной системы узлов построим интерполяционный полином Лагранжа.
x y
L ( x)
x x x
n
n
где
n 1
i
'
i 0
i
n 1
i
x x x x x ... x x ,
n 1
0
1
n
x x x x x ... x x x x ... x x .
'
n 1
i
i
0
i
1
Тогда в силу единственности решения
i
i 1
i
i 1
Ln xi yi
i
n
i 0,1,2,..., n
5

6.

x x0
q,
h
Полагая
x h
получим
n 1
n 1
q q 1 q 2 ... q n h n 1q n ,
x x x x x ... x x x x ... x x
'
n 1
i
i
0
i
1
i
i 1
i
i 1
i
n
h ni i 1 i 2 ...1 1 ... n i 1 h ni ! n i ! .
n i
n i
1 yi q n 1
Ln x
i 0 i ! n i ! q i
n
Тогда, для полинома Лагранжа имеем
Учитывая то, что
Получаем
dx
h
dq
n i
dy d
1 n 1 yi d q n 1
Ln x
dx dx
h i 0 i ! n i ! dx q i
6

7.

Погрешность вычисления
первой производной:
Для Rn получим
rn ( x)
d
d
dy d
Rn ( x) y ( x) Ln ( x)
Ln ( x),
dx
dx
dx dx
y n 1
x ,
Rn ( x)
n 1
n 1 !
где ξ = ξ(x) – промежуточное значение между точками x0,x1,…xn ,
y(n) – n-ая производная по x.
Тогда,
'
d
1
d n 1
Rn x
y n 1 n 1 x n 1 x
y
dx
dx
n 1 !
1 h n
n i
i ! n i ! n 1
y
n 1 !
Если число узлов нечетно и производная берется в средней точке, то выражение
для численного дифференцирования получается более просто и имеет
повышенную точность.
7

8.

5.4. Приближенное интегрирование функций. Общие замечания
Если функция f(x) непрерывна на отрезке [a,b] и известна ее первообразная F(x),
то определенный интеграл в пределах от a до b может быть вычислен по формуле
Ньютона-Лейбница
b
f ( x)dx F (b) F (a)
a
где
dF
f (x),
dx
Приближенные и, в первую очередь, численные методы вычисления определенных
интегралов применяются, когда:
1. Первообразная не может быть найдена аналитически или имеет очень сложный вид,
2. f(x) задана таблично (само понятие первообразной теряет смысл).
Задача численного интегрирования заключается в нахождении определенного
интеграла на основе ряда значений подынтегральной функции.
Численное вычисление однократного интеграла называется механической квадратурой,
двойного – механической кубатурой. Соответствующие формулы называются
квадратурными и кубатурными.
8

9.

5.6. Метод прямоугольников
Пусть функция y = f(x) задана на отрезке [a,b]. С помощью точек x0, x1…, xn
разобьем этот отрезок на n элементарных отрезков [xi-1, xi] (i = 1, 2,..., n) причем
x0= a, xn= b. На каждом из этих отрезков выберем точку ξi = xi-1 или ξi = xi и
найдем произведение si = f(ξi )∙(xi - xi-1) = f(ξ i )∙Δxi . Сумма этих произведений
является приближенным значением определенного интеграла
b
n
n
f ( x)dx s f ( ) x
i
i 1
a
i
i 1
Sn ,
i
Более точным является метод, называемый методом средних и использующий
значение функции в средних точках элементарных отрезков (в полуцелых узлах)
b
n
f ( x)dx f ( x
a
i 1
i
n
1
2
) xi f ( x
Если шаг задания узлов hi постоянный,
формула приобретает вид
i 1
i
b
1
2
)hi ,
(i=1, 2,..., n)
n
f ( x)dx h f ( x
a
i 1
i
1
2
).
9

10.

5.7. Метод трапеций
В этом методе используется линейная интерполяция функции y=f(x) в промежутках
между узлами.
b
n
f ( x)dx
a
i 1
f ( xi 1 ) f ( xi )
hi ,
2
(i=1, 2,..., n)
При постоянном шаге интерполяции
b
a
f ( x0 ) f ( xn ) n 1
f ( x)dx h
f ( xi ) .
2
i 1
10

11.

5.8. Уточненные значения интегралов
Погрешность численного метода в общем случае равна
b
Rn f ( x)dx S n
a
Главный член погрешности интеграла (I1), полученного методом прямоугольников
на отрезке [xi-1, xi]:
1
hi f 4 x 1 ,
24
i 2
а интеграла (I2), полученного методом трапеций, примерно в 2 раза больше и имеет
противоположный знак:
1 3 d 2 f ( xi )
hi
,
2
24
dx
На основании этого можно записать уточненную формулу для вычисления
определенного интеграла с использованием значений I1 и I2:
b
a
f x dx
2 I1 I 2
3
11

12.

5.9. Метод парабол (метод Симпсона)
Разобьем отрезок интегрирования [a,b] на четное число n равных частей с шагом h.
На каждом отрезке [x0, x2], [x2, x4],..., [xi-1, xi] подынтегральную функцию заменим
интерполяционным многочленом второй степени:
f ( x) i ( x) ai x 2 bi x ci .
(xi-1≤x≤xi)
В качестве φi(x) можно взять интерполяционный многочлен Лагранжа, проходящий
через точки Mi-1(xi-1,yi-1) , Mi(xi,yi) , Mi+1(xi+1,yi+1) :
i ( x)
x xi x xi 1 y x xi 1 x xi 1 y x xi 1 x xi y
xi 1 xi xi 1 xi 1 i 1 xi xi 1 xi xi 1 i xi 1 xi 1 xi 1 xi i 1
Элементарная площадь может быть вычислена аналитически и с учетом того,
что шаг интерполирования h постоянный, получим следующее выражение
xi 1
si i ( x)dx
xi 1
1
2
2h
xi 1
( x x )( x x
i
xi 1
i 1
) yi 1 2( x xi 1 )( x xi 1 ) yi ( x xi 1 )( x xi ) yi 1 dx
h
yi 1 4 yi yi 1
3
12

13.

Просуммировав все отрезки, получим:
n
h
S si y0 4 y1 2 y2 4 y3 2 y4 ... 2 yn 2 4 yn 1 yn .
3
i 1
Значение S принимается в качестве определенного интеграла. Окончательное
выражение для формулы Симпсона имеет вид :
b
a
n
n
1
2
2
h
f ( x)dx y0 4 y 2 j 1 2 y 2 j yn
3
j 1
j 1
Точность метода Симпсона составляет 6 знаков. Главный член погрешности этого
метода
h 4 d 4 f x
Rn
180 dx 4
имеет тот же порядок, что и комбинированный метод прямоугольников и трапеций, т.е.
на порядок лучше, чем для отдельно взятых методов прямоугольников и трапеций.
13

14.

5.10. Формула Ньютона-Кортеса
Пусть для данной функции y=f(x) необходимо вычислить определенный интеграл.
b
f x dx
a
Разобьем отрезок интегрирования [a,b] на n равных частей с шагом h. Будем считать,
что функция задана в узлах yi=f(xi), i=0, 1, 2,…, n. Заменим подынтегральную
функцию интерполяционным полиномом Лагранжа и получим приближенную
b
квадратурную формулу:
n
f ( x)dx A y ,
a
i 0
i
i
где Ai – некоторые постоянные коэффициенты. Введем обозначения
q
x x0
,
h
q n 1` q q 1 ... q n
n i
1
q n 1
Ln ( x)
yi
i 0 i ! n i ! q i
n
и представим полином Лагранжа в виде
14

15.

Ai b a H i ,
где постоянные коэффициента Hi называются коэффициентами Кортеса :
1 1
q n 1
Hi
dq.
n i ! n i ! 0 q i
n i
n
Коэффициенты Кортеса обладают следующими свойствами:
n
H
i 0
i
H i H n i .
1,
Окончательный вид квадратурной формулы Ньютона-Кортеса:
b
n
ydx b a H y ,
a
где
h
b a
,
n
i 0
i
i
yi f i a ih ,
i 0, 1, 2,..., n
Формулы методов прямоугольника, трапеций и Симпсона являются частыми
случаями формулы Ньютона-Кортеса.
15

16.

Формулы методов прямоугольника, трапеций и Симпсона являются частыми
случаями формулы Ньютона-Котеса.
Остаточный член формулы Ньютона-Котеса:
2 E n2 3
Rn O h ,
где E(n/2) – целая часть дроби n/2. Таким образом, нечетное число ординат
является более выигрышным.
16

17.

5.11. Квадратурная формула Гаусса
1
Полиномы Лежандра
P2(x)
y
n
1 dn
2
Pn ( x) n
x 1 ,
n
2 n ! dx
P0(x)
P3(x)
0
n 0, 1, 2, ...
P4(x)
P1(x)
Важные свойства полиномов:
Pn 1 1,
Pn 1 1 ,
n
-1
-1
0
1
x
1
P x Q dx 0,
n
k
где Qk –любой полином степени k<n.
1
Полином Лежандра Pn(x) имеет n действительных корней в интервале (-1, 1).
17

18.

Рассмотрим функцию f(t), заданную на отрезке [-1, 1].
Постановка задачи: подобрать точки t1, t2, …tn и коэффициенты А1, А2,…Аn, чтобы
квадратурная формула
1
n
f t dt A f t 1
i 1
1
i
i
была точной для всех полиномов f(t) наивысшей возможной степени N . Так как у нас
2n неизвестных, а полином степени 2n-1 определяется 2n коэффициентами, то
высшая степень полинома N = 2n-1. Для обеспечения приведенного сверху
равенства необходимо и достаточно, чтобы оно было верным при
f t 1, t , t 2 ,..., t 2 n 1.
1 1
k
1t dt k 1
1
Учитывая соотношение
k 1
2
k 1
0
при k четном,
при k нечетном
18

19.

Для решения поставленной задачи достаточно
определить t1, t2, …tn и А1, А2,…Аn, из
нелинейной системы 2n уравнений :
Ai 2
i 1
n
A
t
0
i i
i 1
................
n
2
2n 2
A
t
i i
2
n
1
i 1
n
2 n 1
A
t
0
i i
i 1
n
(2)
Далее применяется искусственный прием. Рассмотрим полиномы f(t),
сконструированные в том числе из полиномов Лежандра
f t t k Pn (t ),
k 0, 1,...,
n 1
Так как степени этих полиномов не превышают 2n-1, то на основании
системы (2) для них должна быть справедлива формула (1).
19

20.

Подстановка в эту формулу f(t) дает :
1
n
k
t
P
(
t
)
dt
A
t
i i Pn (ti ),
n
k
1
k 0, 1,...,
n 1
i 1
В силу ортогональности полиномов Лежандра
1
k
t
Pn (t )dt 0,
k n
1
или
n
k
A
t
i i Pn (ti ) 0
k 0, 1,...,
n 1
i 1
Это равенство будет заведомо справедливо , если положить
Pn (ti ) 0
i 1,..., n .
То есть для достижения наивысшей точности квадратурной формулы в качестве ti
взять нули соответствующих полиномов Лежандра, далее, подставив их в систему
(2), которая относительно Ai будет линейной, найти эти коэффициенты.
Подстановка найденных значений ti и Ai в выражение (1) даст квадратурную
формулу Гаусса.
20

21.

5.12. Дифференцирование и интегрирование в пакете MathCad
21

22.

6. Численное решение дифференциальных уравнений
6.1. Основные понятия
Дифференциальные уравнения делятся на:
1. обыкновенные (содержащие одну переменную),
2. уравнения в частных производных.
Обыкновенные дифференциальные уравнения содержат одну или несколько
производных искомой функции y=y(x) и могут быть записаны в виде
dy d 2 y d n y
F x, y, , 2 ,... n 0.
dx dx
dx
Наивысший порядок n входящей в уравнение производной называется порядком
дифференциального уравнения.
Уравнение, имеющее вид
dny
dy d 2 y d n 1 y
F x, y, , 2 ,... n 1 ,
dx n
dx dx
dx
называется уравнением, разрешенным относительно старшей производной.
22

23.

Линейными дифференциальными уравнениями называются уравнения, линейные
относительно искомой функции и её производных.
Решением дифференциального уравнения всякая функция y = φ(x), которая
после её подстановки в уравнение, превращает его в тождество. Графическое
представление решения – интегральная кривая.
Общее решение обыкновенного дифференциального уравнения порядка n
содержит n постоянных C1, C2,…Cn.
Частное решение дифференциального уравнения получается из общего, если
произвольным константам придать определенные значения.
Геометрическая интерпретация линейного дифференциального уравнения
первого порядка. Поскольку производная характеризует наклон касательной к
интегральной кривой в данной точке, то при dy/dx=k получаем уравнение линии
постоянного наклона, называемой изоклиной. Меняя k, получаем семейство изоклин.
Общее решение описывает бесконечное семейство интегральных кривых с
параметром С, а частному решению соответствует одна кривая этого семейства.
23

24.

Для выделения некоторого частного решения линейного дифференциального
уравнения первого порядка достаточно задать координаты некоторой точки (x0,y0)
на данной интегральной кривой.
Для выделения частного решения из общего решения дифференциального
уравнения порядка n следует задать столько дополнительных условий, сколько
произвольных постоянных C1, C2,…Cn в общем решении.
В зависимости от способа задания дополнительных условий существуют два типа
задач:
1.
Задача Коши: дополнительные условия задаются в одной точке (начальной
точке) и называются начальными условиями.
2.
Краевая задача: дополнительные условия задаются более, чем в одной точке
(как правило, на границах области существования решения), называются
граничными или краевыми условиями.
24

25.

6.2.1. Метод Эйлера
Решить дифференциальное уравнение
dy/dx=f(x,y)
численным методом , значит для заданной последовательности аргументов x0, x1,…,
xn и числа y0 , не определяя функцию y=F(x), найти такие значения y1,…, yn , что
yi=F(xi) (i=1,2,…n) и y0=F(x0).
Разобьем отрезок [a,b] на n равных частей и получим последовательность значений
аргумента xi = x0 +i∙h, где h - шаг интегрирования. Будем считать, что x0 и y0 заданы.
Функцию y=F(x) можно разложить в ряд Тейлора и, с точностью до членов O(h2),
записать
y1 y0 h f x0 , y0 ,
y2 y1 h f x1 , y1 ,
.................................
yn yn 1 h f xn 1 , yn 1
25

26.

6.2.1. Метод Рунге-Кутта
Этот метод является методом повышенной точности. Как и в методе Эйлера
yi= yi-1+Δyi-1,
(i=1,2,…n)
но функцию y=F(x) раскладывают в ряд Тейлора с точностью до членов h4,
включительно.
dy h 2 d 2 y h 3 d 3 y h 4 d 4 y
y y ( x h) y ( x) h
dx 2 dx 2 6 dx 3 24 dx 4
Производные dky/dxk определяются последовательным дифференцированием
уравнения dy/dx=f(x,y).
Вместо непосредственных вычислений производных в методе Рунге-Кутта
определяются 4 числа:
h
k
k3 h f ( x , y 2 ),
k1 h f ( x, y ),
2
h
k
k 2 h f ( x , y 1 ),
2
2
В результате :
yi
2
k
k 4 h f ( x h, y 3 ).
2
1
k1 xi , yi 2k2 xi , yi 2k3 xi , yi k4 xi , yi
6
26

27.

6.2.2. Метод Рунге-Кутта в пакете MathCad
27

28.

6.3. Приближенные методы решения краевой задачи
для линейного обыкновенного дифференциального
уравнения
2
d Y
dY
p
x
q x Y f x .
2
dx
dx
Рассмотрим уравнение вида
Краевая задача состоит в отыскании решения Y=Y(x) на отрезке [a,b],
удовлетворяющего граничным условиям
Y(a)=A, Y(b)=B
Для нахождения приближенного решения выбирается линейно независимая
(базисная) система дважды дифференцируемых функций φ0(x), φ1(x), φ2(x),…,
φn(x). При этом φ0(x) удовлетворяет данным граничным условиям, а φ1(x),
φ2(x),…, φn(x) – однородным. Искомое решение представляется в виде линейной
n
комбинации:
y ( x) 0 ( x) ai i ( x)
i 1
Невязка :
d2y
dy
x, a1 ,..., an 2 p q y f x
dx
dx
Коэффициента ai стараются подобрать так, чтобы невязка была минимальной.
28

29.

6.4.1. Метод коллокаций
В этом методе выбираются n точек xi, принадлежащих отрезку [a,b], называемых
точками коллокации , невязки ψ(x,a1,a2,…,an) в которых приравниваются нулю. В
результате получается система n алгебраических уравнений относительно
коэффициентов ai.
6.4.2. Метод наименьших квадратов
Основан на минимизации суммы квадратов невязок в заданной системе точек xi,
принадлежащих отрезку [a,b]. Из этого условия также получается система n
алгебраических уравнений относительно коэффициентов ai.
6.4.3. Метод Галеркина
Основан на требовании ортогональности базисных функций к невязке, которое
выражается в виде
b
x, a ,..., a x dx, i 1,2,..., n
1
n
i
a
29

30.

6.4.4. Метод стрельбы
Рассмотрим краевую задачу для уравнения второго порядка, разрешенного
относительно второй производной.
d 2Y
dY
f
x
,
Y
,
dx 2
dx
Решение будем искать на отрезке [0,1]. Граничные условия:
Y (0) y0 , Y (1) y1
Сущность метода стрельбы заключается в сведении краевой задачи к задаче Коши с
начальными условиями:
dY
Y ( 0) y 0 ,
dx
k tg
x 0
Считая решение задачи Коши Y=Y(x,α), зависящим от параметра α, ищется такая
интегральная кривая, которая выходит из точки (0,y0) и попадает в точку (1,y1). На
основании чего можно записать уравнение относительно α :
Y (1, ) y1 F 0
и решить его любым методом (например, делением отрезка пополам).
30

31.

The end
31