Аппроксимацией(приближением) функции называется нахождение такой функции
(аппроксимирующей функции), которая была бы близка заданной.
Критерии близости функций и могут быть различные.
Основная задача аппроксимации — построение приближенной
(аппроксимирующей) функции, в целом наиболее близко проходящей около данных точек или около данной непрерывной функции. Такая задача возникает при наличии погрешности в исходных данных (в этом случае нецелесообразно проводить функцию точно через все точки, как в интерполяций) или при желании получить упрощенное математическое описание сложной или неизвестной зависимости.
Концепция аппроксимации
48
Близость исходной и аппроксимирующей функций определяется числовой мерой
— критерием аппроксимации (близости). Наибольшее распространение получил квадратичный критерий, равный сумме квадратов отклонений расчетных значений от "экспериментальных" (т.е. заданных), — критерий близости в заданных точках:
Здесь у
i
— заданные табличные значения функции; у
i
расч
— расчетные значения по аппроксимирующей функции; b
i
— весовые коэффициенты, учитывающие относительную важность i-и точки (увеличение b,. приводит при стремлении уменьшить R куменьшению, прежде всего отклонения в i-й точке, так как это отклонение искусственно увеличено за счет относительно большого значения весового коэффициента).
Квадратичный критерий обладает рядом "хороших" свойств, таких, как дифференцируемость, обеспечение единственного решения задачи аппроксимации при полиномиальных аппроксимирующих функциях.
Другим распространенным критерием близости является следующий:
Этот критерий менее распространен в связи с аналитическими и вычислительными трудностями, связанными с отсутствием гладкости функции и ее дифференцируемости.
Выделяют две основные задачи:
1) получение аппроксимирующей функции, описывающей имеющиеся данные, с погрешностью не хуже заданной;
2) получение аппроксимирующей функции заданной структуры с наилучшей возможной погрешностью.
Чаще всего первая задача сводится ко второй перебором различных аппроксимирующих функций и последующим выбором наилучшей.
24.
Метод аппроксимации алгебраическим полиномом. Алгоритм метода.
При анализе эмпирических данных возникает необходимость найти в явном виде функциональную зависимость между величинами x и y , которые получены в
49 результате измерений Как правило общий вид этой функциональной зависимости известен а некоторые числовые параметры закона неизвестны
Пусть например функция задана в виде
)
,
(
i
i
y
x
...,
,
2
,
1
n
i
Задача состоит в аппроксимации неизвестной функциональной зависимости между x и y многочленом заданной степени m :
m
j
j
j
m
m
x
p
x
p
x
p
p
x
y
0 1
0
)
(
Для решения этой задачи применяют критерий наименьших квадратов Согласно этому критерию коэффициенты многочлена нужно выбрать такими чтобы сумма квадратов отклонений найденного многочлена от заданных значений функции была минимальной Другими словами коэффициенты
m
p
p
p
,...,
,
1 0
должны минимизировать функцию
n
i
i
m
i
m
i
m
y
x
p
x
p
p
p
p
p
F
1 2
1 0
1 0
)
,....,
,
(
В точке минимума функции F её производные
i
p
F
/
обращаются в нуль
Дифференцируя F и приравнивая нулю производные получим систему линейных алгебраических уравнений:
,
,
1 1
2 2
1 2
1 1
0 1
1 1
1 2
1 3
1 1
2 0
1 1
1 2
1 2
1 1
0
n
i
i
m
i
m
n
i
m
i
n
i
m
i
i
n
i
m
i
n
i
m
i
n
i
i
i
m
n
i
m
i
n
i
i
n
i
i
n
i
i
n
i
i
m
n
i
m
i
n
i
i
n
i
i
y
x
p
x
p
x
p
x
p
x
y
x
p
x
p
x
p
x
p
x
y
p
x
p
x
p
x
np
Отметим что общее количество вычисляемых сумм равно
1 3
m
Получена СЛАУ относительно неизвестных
m
p
p
p
,...,
,
1 0
. Определитель этой системы отличен от нуля В алгоритме, представленном на рис.6.1, для решения методом Гаусса СЛАУ приведена к виду:
50
..,
,
,
1
,
1 1
0 0
1
,
1 1
1 11 0
10 1
,
0 0
1 01 0
00
m
m
m
mm
m
m
m
m
m
m
m
m
a
p
a
p
a
p
a
a
p
a
p
a
p
a
a
p
a
p
a
p
a
51
25.
Формула Ньютона-Котеса для вычисления определенных интегралов.
Частные случаи.
Рис.6.1 – алгоритм аппроксимации степенным полиномом по критерию наименьших квадратов
начало
)
1
)(
1
(
0
,
,
,
n
i
y
x
n
i
i
m
m
r
)
1
(
0
0 1
,
m
r
a
)
1
)(
1
(
0
n
i
i
r
i
m
r
y
x
a
1
,
r
i,
)
2
)(
1
(
0
m
j
j
i
j
x
sum
0
j
sum
)
1
)(
1
(
0
n
i
j
i,
1 1
m
r
)
1
(
0
m
c
)
1
(
0
c
r
rc
sum
a
r
c,
значений
табличных
ввод
полинома
о
симирующег
аппрок
m
порядка
ввод
)
0
(
уравнений
ических
алгебра
линейных
m
из
системы
частей
правых
вычисление
1
sum
массиве
в
в
результато
сохранение
и
системы
частях
левых
в
сумм
вычисление
полинома
симирущего
аппрок
параметров
m
значений
получение
и
уравнений
системы
решение
1
данных
табличных
исходных
график
на
нанесение
полинома
графика
построение
,
конец
системы
частей
правых
и
тов
коэффициен
матрицу
в
sum
массива
из
сумм
х
вычисленны
перенос
52
Формулы Ньютона-Котеса относятся к случаю равноотстоящих узлов интегрирования Они получены как результат интерполяции подынтегральной функции полиномом Лагранжа и в общем случае имеют вид
b
a
n
i
i
i
x
f
H
a
b
dx
x
f
0
),
(
)
(
где
;
ih
a
x
i
;
..,
,
1
,
0
n
i
;
/
)
(
n
a
b
h
n
i
n
i
dq
i
q
n
q
q
q
q
i
n
i
n
H
0
)
).....(
2
)(
1
(
)!
(
!
)
1
(
1
- коэффициенты Котеса
Очевидно, коэффициенты Котеса не зависят от подынтегральной функции.
Для n от 1 до 8 они имеют значения: n=1: H
0
= H
1
= 0.5; n=2: H
0
= H
2
= 0.1666667, H
1
= 0.6666667; n=3: H
0
= H
3
= 0.125, H
1
= H
2
= 0.375; n=4: H
0
= H
4
= 0.0777778, H
1
= H
3
= 0.3555556, H
2
= 0.1333333; n=5: H
0
= H
5
= 0.0659722, H
1
= H
4
= 0.2604166, H
2
= H
3
= 0.1736111; n=6: H
0
= H
6
= 0.0488095, H
1
= H
5
= 0.2571428, H
2
= H
4
= 0.0321428,
H
3
= 0.3238095; n=7: H
0
= H
7
= 0.0434606, H
1
= H
6
= 0.2070023, H
2
= H
5
= 0.0765626,
H
3
= H
4
= 0.1729745; n=8: H
0
= H
8
= 0.0348853, H
1
= H
7
= 0.2076895, H
2
= H
6
= -0.0327336,
H
3
= H
5
= 0.3702292, H
4
= -0.1601410.
Правильность вычисления коэффициентов Котеса можно проверить по их свойствам:
1)
n
i
i
H
0
,
1 2)
i
n
i
H
H
Поскольку интервал интегрирования может быть большим то его делят на k равных частей и к каждой из них применяют соответствующую формулу Формула для
53 всего интервала интегрирования становится составной Обычно используют
3
,
2
,
1
n
n
n
Частные случаи
Формула трапеций (n=1)
b
a
k
k
f
f
f
f
f
h
dx
x
f
,
2 2
2 2
)
(
1 2
1 0
где
k
a
b
h
- шаг интегрирования
);
(
i
i
x
f
f
;
0
ih
x
x
i
;
0
a
x
...,
,
1
,
0
k
i
2. Формула Симпсона (n=2)
,
2 4
3
)
(
2 2
4 2
1 2
3 1
2 0
b
a
k
k
k
f
f
f
f
f
f
f
f
h
dx
x
f
где
;
2k
a
b
h
h
2
- шаг интегрирования
);
(
i
i
x
f
f
;
0
ih
x
x
i
;
0
a
x
2
...,
,
1
,
0
k
i
3. Формула Ньютона (“трех восьмых”, n=3)
,
3 2
8 3
)
(
1 3
2 3
2 1
3 3
6 3
3 0
b
a
k
k
k
k
f
f
f
f
f
f
f
f
f
h
dx
x
f
где
;
3k
a
b
h
h
3 - шаг интегрирования
);
(
i
i
x
f
f
;
0
ih
x
x
i
;
0
a
x
3
...,
,
1
,
0
k
i
В этих формулах шаг интегрирования постоянный и поэтому они эффективны в методах с автоматическим выбором шага интегрирования
Процедура автоматического выбора шага интегрирования основана на двукратном вычислении интеграла на каждом шаге: вначале с h (значение интеграла
1
S
), а затем - с
2
/
h
(значение интеграла
2
S
). Если
2 1
S
S
меньше заданной наперед ошибки то считается что интегрирование на данном шаге выполнено правильно В
54 противном случае h делится пополам и повторяется описанная выше процедура интегрирования
26.
Формула Гаусса для вычисления определенных интегралов. Формулы
прямоугольников – срединных, левых, правых.
Для увеличения точности здесь взяты неравноотстоящие узлы
i
x
, что оказалось возможным при интерполировании подынтегральной функции полиномом Лежандра:
,
)
(
)
(
1
b
a
n
i
i
i
x
f
A
a
b
dx
x
f
где
;
i
i
X
a
b
a
x
...,
,
2
,
1
n
i
Ниже приведены абсциссы
i
X и коэффициенты
i
A
формулы Гаусса для значений пределов интегрирования
,
0
a
,
1
b
при которых формула Гаусса имеет вид
)
(
)
(
1 0
1
n
i
i
i
x
f
A
dx
x
f
Абсциссы и коэффициенты формулы Гаусса
n=1: X
1
= 0.5; A
1
= 1.0; n=2: X
1
= 0.2113249, X
2
= 0.7886751; A
1
= A
2
= 0.5; n=3: X
1
= 0.1127017, X
2
= 0.5, X
3
= 0.8872983;
A
1
= A
3
= 0.2777778, A
2
= 0.4444444; n=4: X
1
= 0.0694318, X
2
= 0.3300095, X
3
= 0.6699905, X
4
= 0.9305682;
A
1
= A
4
= 0.1739274, A
2
= A
3
= 0.3260726; n=5: X
1
= 0.0469101, X
2
= 0.2307653, X
3
= 0.5, X
4
= 0.7692347,
X
5
= 0.9530899;
A
1
= A
5
= 0.1184634, A
2
= A
4
= 0.2393143, A
3
= 0.2844444; n=6: X
1
= 0.0337652, X
2
= 0.1693953, X
3
= 0.3806904, X
4
= 0.6193096,
55
X
5
= 0.8306047, X
6
= 0.9662348;
A
1
= A
6
= 0.0856622, A
2
= A
5
= 0.1803808, A
3
= A
4
= 0.2339570; n=7: X
1
= 0.0254460, X
2
= 0.1292344, X
3
= 0.2970774, X
4
= 0.5,
X
5
= 0.7029226, X
6
= 0.8707656, X
7
= 0.9745540;
A
1
= A
7
= 0.0647425, A
2
= A
6
= 0.1398527,
A
3
= A
5
= 0.1909150, A
4
= 0.2089796; n=8: X
1
= 0.0198551, X
2
= 0.1016668, X
3
= 0.2372338, X
4
= 0.4082827,
X
5
= 0.5917173, X
6
= 0.7627662, X
7
= 0.8983332, X
8
= 0.9801449;
A
1
= A
8
= 0.0506143, A
2
= A
7
= 0.1111905,
A
3
= A
6
= 0.1568533, A
4
= A
5
= 0.1813419.
Если отрезок интегрирования
a
b разбить на k равных частей и к каждой из них применить формулу Гаусса для фиксированного n , то
,
)
(
)
(
)
(
1 1
)
(
1 1
k
j
n
i
j
i
i
b
a
k
j
x
x
x
f
A
h
dx
x
f
dx
x
f
j
j
где
k
a
b
h
/
- шаг интегрирования;
)
( j
i
x
- i -я абсцисса j -го шага интегрирования;
,
0
a
x
,
1
h
x
x
j
j
;
...,
,
2
,
1
k
j
,
1
)
(
h
X
x
x
i
j
j
i
;
...,
,
2
,
1
n
i
i
i
A
X ,
- абсциссы и коэффициенты формулы Гаусса
56
При одинаковом числе ординат формула Гаусса дает большую точность чем другие формулы но она менее эффективна в методах с автоматическим выбором шага интегрирования
Самой простой квадратурной формулой является формула прямоугольников (
1
n
5 0
1
X
0 1
1
A
)
b
a
b
a
f
a
b
dx
x
f
2
)
(
Если отрезок интегрирования большой и разбит на k равных частей то составная формула прямоугольников имеет вид
,
)
(
)
(
1
b
a
k
j
i
x
f
h
dx
x
f
где
k
a
b
h
/
;
2
/
1
h
a
x
h
x
x
j
j
1
k
j
...,
,
3
,
2
На практике используются также модификации общей формулы: формула
“левых” прямоугольников (
a
x
1
) и формула “правых” прямоугольников (
h
a
x
1
).
27.
Алгоритм получения зависимости от шага интегрирования
фактической ошибки вычисления интеграла по формуле прямоугольников
57
Рис.7.1 – алгоритм получения зависимости фактической ошибки интегрирования функции по формуле прямоугольников от шага интегрирования
начало
b
a,
)
(
)
(
a
F
b
F
s
20
)
1
(
1
k
0
,
y
k
a
b
h
h
x
b
x
h
a
x
,
,
2
/
)
(x
f
y
x
h
y
*
h
y
s
,
k
конец
58
28.
Алгоритм получения зависимости от шага интегрирования
фактической ошибки вычисления интеграла по формуле Гаусса.
Рис.7.2 – алгоритм получения зависимости фактической ошибки интегрирования функции методом Гаусса от шага интегрирования
начало
)
1
)(
1
(
0
,
,
,
,
,
n
i
X
A
n
b
a
i
i
)
(
)
(
a
F
b
F
s
20
)
1
(
1
k
0
,
y
k
a
b
h
h
x
h
b
x
a
x
;
2
/
;
)
(
i
i
hX
x
f
A
y
x
i,
h
y
*
h
y
s
,
k
конец
)
1
)(
1
(
0
n
i
59
Алгоритмы для получения зависимости фактической ошибки вычисления интеграла от шага интегрирования
Эта группа алгоритмов представлена на рис.7.1…7.6. Алгоритмы на рис.7.1, 7.3,
7.4, 7.5 соответствуют частным случаям формул Гаусса и Ньютона-Котеса. Алгоритмы на рис.7.2, 7.6 построены по формулам Гаусса и Ньютона-Котеса в общем случае.
Исходными данными являются: подынтегральная
)
(x
f
и первообразная
)
(x
F
функции, пределы интегрирования
b
a,
(вводятся). В алгоритме использованы переменные: s – точное значение интеграла, вычисленное по формуле Ньютона-
Лейбница, k - число шагов интегрирования, h - шаг интегрирования, y - приближенное значение интеграла. Как результат, выводятся фактическая ошибка вычисления интеграла
y
s и шаг интегрирования h (или
nh
).
29.
Алгоритм получения зависимости от шага интегрирования
фактической ошибки вычисления интеграла в частных случаях формулы
Ньютона-Котеса.
Ааааа, обоснование алгоритмов смотри в вопросе 25
60
Рис.7.3 – алгоритм получения зависимости фактической ошибки интегрирования функции по формуле трапеций от шага интегрирования
начало
b
a,
)
(
)
(
a
F
b
F
s
20
)
1
(
1
k
0
,
y
k
a
b
h
h
x
h
b
x
a
x
;
2
/
;
)
(
)
(
h
x
f
x
f
y
x
2
/
*
h
y
h
y
s
,
k
конец
61
Рис.7.4 – алгоритм получения зависимости фактической ошибки интегрирования функции по формуле Симпсона от шага интегрирования
начало
b
a,
)
(
)
(
a
F
b
F
s
20
)
1
(
1
k
0
,
2
y
k
a
b
h
h
x
h
b
x
a
x
2
;
;
)
2
(
)
(
4
)
(
h
x
f
h
x
f
x
f
y
x
3
/
*
h
y
h
y
s
2
,
k
конец
62
30.
Алгоритм получения зависимости от шага интегрирования
фактической ошибки вычисления интеграла по формуле Ньютона-Котеса.
Рис.7.5 – алгоритм получения зависимости фактической ошибки интегрирования функции по формуле Ньютона от шага интегрирования
начало
b
a,
)
(
)
(
a
F
b
F
s
20
)
1
(
1
k
0
,
3
y
k
a
b
h
h
x
h
b
x
a
x
3
;
2
;
)
3
(
)
2
(
3
)
(
3
)
(
h
x
f
h
x
f
h
x
f
x
f
y
x
8
/
3
*
h
y
h
y
s
3
,
k
конец
63
Алгоритмы для получения зависимостей фактической ошибки и временных затрат вычисления интеграла от задаваемой ошибки интегрирования
Рис.7.6 – алгоритм получения зависимости фактической ошибки интегрирования функции методом Ньютона-Котеса от шага интегрирования
начало
n
i
H
n
b
a
i
)
1
(
0
,
,
,
,
)
(
)
(
a
F
b
F
s
20
)
1
(
1
k
0
,
y
kn
a
b
h
nh
x
h
n
b
x
a
x
;
)
2
(
;
)
(
ih
x
f
H
y
i
x
i,
nh
y
*
nh
y
s
,
k
конец
n
i
)
1
(
0
64
Эти алгоритмы (рис.7.7…7.10) являются адаптивными, т.к. здесь шаг интегрирования автоматически устанавливается, начиная со значения
a
b
h
, таким, чтобы выполнилось условие:
y
ys
y
- относительная ошибка интегрирования меньше заданной
. Для изменения величины шага используется правило Рунге. Через
ys обозначено предыдущее значение интеграла; в начале вычислений оно может принято равным любому числу, здесь – нулю. Затраты машинного времени в алгоритме оцениваются косвенно – по числу вызовов подынтегральной функции
)
(x
f
, для чего использована переменная
t
Адаптивные алгоритмы для других частных случаев можно легко составить по аналогии с представленными алгоритмами.
31.
Адаптивный алгоритм получения зависимостей от задаваемой
ошибки затрат машинного времени и фактической ошибки вычисления
интеграла по формуле прямоугольников.
Хз че это
65
Рис.7.7 – адаптивный алгоритм получения зависимостей фактической ошибки и временных затрат при интегрировании функции по формуле прямоугольников от задаваемой ошибки
начало
b
a,
)
(
)
(
a
F
b
F
s
10
*
;
011 0
;
10 5
0
,
ys
a
b
h
h
x
b
x
t
y
h
a
x
;
;
0
,
0
,
2
/
t
x
f
y
),
(
x
h
y
*
t
y
s
,
,
конец
y
ys
y
2
/
h
h
y
ys
да
66
32.
Адаптивный алгоритм получения зависимостей от задаваемой
ошибки затрат машинного времени и фактической ошибки вычисления
интеграла по формуле Гаусса.
И это хз что
67
Рис.7.8 – адаптивный алгоритм получения зависимостей фактической ошибки и временных затрат при интегрировании функции методом Гаусса от задаваемой ошибки
начало
)
1
)(
1
(
0
,
,
,
,
,
n
i
X
A
n
b
a
i
i
)
(
)
(
a
F
b
F
s
10
*
;
011 0
;
10 5
0
,
ys
a
b
h
h
x
h
b
x
t
y
a
x
;
2
/
;
0
,
0
,
t
hX
x
f
A
y
i
i
),
(
x
i,
h
y
*
t
y
s
,
,
конец
y
ys
y
2
/
h
h
y
ys
да
)
1
)(
1
(
0
n
i
68
33.
Адаптивный алгоритм получения зависимостей от задаваемой
ошибки затрат машинного времени и фактической ошибки вычисления
интеграла в частных случаях формулы Ньютона-Котеса.
Опять хз что это
69
34.
Адаптивный алгоритм получения зависимостей от задаваемой
ошибки затрат машинного времени и фактической ошибки вычисления
интеграла по формуле Ньютона-Котеса.
Рис.7.9 – адаптивный алгоритм получения зависимостей фактической ошибки и временных затрат при интегрировании функции по формуле трапеций от задаваемой ошибки
начало
b
a,
)
(
)
(
a
F
b
F
s
10
*
;
011 0
;
10 5
0
,
ys
a
b
h
h
x
h
b
x
t
y
a
x
;
2
/
;
0
,
0
,
2
),
(
)
(
t
h
x
f
x
f
y
x
2
/
*
h
y
t
y
s
,
,
конец
y
ys
y
2
/
h
h
y
ys
да
70
Рис.7.10 – адаптивный алгоритм получения зависимостей фактической ошибки и временных затрат при интегрировании функции методом Ньютона-Котеса от задаваемой ошибки
начало
n
i
H
n
b
a
i
)
1
(
0
,
,
,
,
)
(
)
(
a
F
b
F
s
10
*
;
011 0
;
10 5
0
,
/
)
(
ys
n
a
b
h
nh
x
nh
b
x
t
y
a
x
;
2
/
;
0
,
0
,
t
ih
x
f
H
y
i
),
(
x
i,
nh
y
*
t
y
s
,
,
конец
y
ys
y
2
/
h
h
y
ys
да
n
i
)
1
(
0
71
35.
Численные методы интегрирования ОДУ и систем ОДУ
(перечислить). Метод интегрирования ОДУ с помощью ряда Тейлора.
Численные методы интегрирования:
•
Метод Эйлера
•
Эйлера с пересчетом
•
Эйлера Ломаных
•
Метод Рунге-Кутта третьего порядка РК3
•
Методы прогноза и коррекции РК4
Метод интегрирования ОДУ с помощью ряда Тейлора
Этот метод теоретически годится для решения любых ДУ но практического интереса не представляет Однако ценность его в том что он может быть использован в качестве эталона по которому сравниваются практически удобные методы Эйлера и
Рунге-Кутта
Разложение функции
)
(x
y
в ряд Тейлора в окрестности точки
m
x
x
имеет вид
0
)
(
3 2
!
)
(
!
3
)
(
!
2
)
(
)
(
)
(
k
k
m
k
m
m
m
m
m
m
m
m
y
k
x
x
y
x
x
y
x
x
y
x
x
y
x
y
Если значения x расположены на расстоянии h друг от друга то
!
!
3
!
2
)
(
0
''
'
3
''
2
'
1
k
m
k
k
m
m
m
m
m
y
k
h
y
h
y
h
hy
y
y
(2)
Ряд Тейлора позволяет получить решение в узле
1
x
по начальным условиям задачи
Коши:
)
,
(
0 0
0
'
0 0
1
y
x
hf
y
hy
y
y
а по решению в произвольном узле
m
x
- решение в узле
1
m
x
:
).
,
(
'
1
m
m
m
m
m
m
y
x
hf
y
hy
y
y
(3)
Сравнение (2) и (3) приводит к выводу что решение получено при отбрасывании членов ряда начиная с члена содержащего
2
h
Следовательно ошибка ограничения равна
)
(
2
h
O
Для уменьшения ошибки ограничения необходимо вычислять значения первой и последующих производных функции
)
,
(
y
x
f
в узле
m
x
что требует выполнения практически неприемлемого объема вычислений
72
36.
Метод Эйлера и модифицированные методы Эйлера интегрирования
ОДУ. Иллюстрации методов
Метод Эйлера основан на формуле (3) Иллюстрация метода Эйлера представлена на рис.8.1;
0
l
- касательная к
)
(x
y
в точке
)
,
(
0 0
y
x
и, следовательно, её наклон определяется значением производной
)
,
(
0 0
y
x
f
Погрешность метода пропорциональна h
2
и следовательно слишком велика при допустимой величине h
Для повышения точности метод Эйлера модифицируют путем использования значений функции
)
,
(
y
x
f
в дополнительных точках
Первый модифицированный метод Эйлера (“с пересчетом”)
Этот метод основан на формулах:
).
,
(
),
,
(
,
2
/
)
(
1 2
1 2
1 1
hk
y
h
x
f
k
y
x
f
k
где
k
k
h
y
y
m
m
m
m
m
m
(4)
Иллюстрация метода представлена на рис.8.2, где
1 0
, l
l
- касательные, наклон которых определяется значениями
2 1
, k
k
соответственно, наклон прямой
2
l
и параллельной ей
*
l определяется значением
2
/
)
(
2 1
k
k
. Отметим, что прямая
2
l
является биссектрисой угла, образованного прямыми
1 0
, l
l
Рис.8.1 – иллюстрация метода Эйлера
y
x
0
y
0
x
1
x
1
y
)
,
(
0 0
y
x
h
решение
вестное
неиз
но
точное
x
y
,
)
(
))
,
(
*
0
,
(
0 0
0
y
x
f
h
y
h
x
1 1
)
(
y
x
y
чения
ограни
ошибка
0
l
73
Иллюстрация метода наглядно доказывает уменьшение ошибки ограничения по сравнению с методом Эйлера.
Теперь докажем факт уменьшения ошибки ограничения, показав, как этот метод согласуется с разложением в ряд Тейлора. Итак, разложим
)
,
(
y
x
f
в окрестности точки
)
,
(
m
m
y
x
в ряд Тейлора:
)
(
)
(
)
,
(
)
,
(
y
m
x
m
m
m
f
y
y
f
x
x
y
x
f
y
x
f
, где
y
y
x
f
f
x
y
x
f
f
y
x
)
,
(
,
)
,
(
вычислены в точке
)
,
(
m
m
y
x
. Обозначив
)
,
(
m
m
y
x
f
через
f
, получим
)
(
)
,
(
2
'
h
O
hff
hf
f
hy
y
h
x
f
y
x
m
m
m
, где
)
(
2
h
O
- сумма остальных членов ряда. Подставим разложение в (4) и получим
)
(
)
(
2 3
2 1
h
O
ff
f
h
hf
y
y
y
x
m
m
Полученный результат позволяет сделать следующие выводы.
1)
Метод «с пересчетом» согласуется с разложением в ряд Тейлора вплоть до членов степени
2
h .
Рис.8.2 – иллюстрация первого модифицированного метода Эйлера (метода «с пересчетом»)
y
x
0
y
0
x
1
x
1
y
)
,
(
0 0
y
x
h
решение
вестное
неиз
но
точное
x
y
,
)
(
)
,
(
*
0
,
(
0 0
0
y
x
f
h
y
h
x
пересчета
после
чения
ограни
ошибка
0
l
1
l
2
l
2
*
|| l
l
74 2)
Порядок метода равен максимальной степени h в его разложении. По сравнению с методом Эйлера порядок увеличился на 1 и стал равным 2.
Поэтому метод «с пересчетом» называют методом Рунге-Кутта 2-го порядка, или РК2.
3)
Увеличение порядка на 1 явилось следствием вычисления производной
)
,
(
y
x
f
в дополнительной точке, т.е. вычисления
).
,
(
1 2
hk
y
h
x
f
k
m
m
Второй модифицированный метод Эйлера (метод “ломаных”)
Метод “ломаных” использует формулы:
).
,
(
),
2
/
,
2
/
(
1
m
m
m
m
m
m
y
x
f
k
hk
y
h
x
hf
y
y
(5)
Иллюстрация метода представлена на рис.8.3, где '
l - касательная к
)
(x
y
в точке
)
2
/
,
2
/
(
0 0
hk
y
h
x
и, следовательно, её наклон определяется значением производной
)
2
/
,
2
/
(
0 0
hk
y
h
x
f
Следовательно, в методе “ломаных” производная вычисляется в дополнительной точке
)
2
/
,
2
/
(
0 0
hk
y
h
x
, которая была получена после вычисления производной в точке
)
,
(
0 0
y
x
, как
)
,
(
0 0
y
x
f
k
. Можно показать, что этот метод также является методом РК2.
Итак, модифицированные методы Эйлера имеют погрешность порядка h
3
. Они относятся к семейству методов Рунге-Кутта второго порядка (РК2)
Рис.8.3 – иллюстрация второго модифицированного метода Эйлера (метод «ломаных»)
y
x
0
y
0
x
1
x
1
y
)
,
(
0 0
y
x
решение
вестное
неиз
но
точное
x
y
,
)
(
)
,
(
*
0
,
(
0 0
0
y
x
f
h
y
h
x
тата
резуль
чения
ограни
ошибка
0
l
'
l
'
||
*
l
l
2
h
2
h
75
37.
Методы Рунге-Кутта интегрирования ОДУ. Достоинства и
недостатки методов. Пример алгоритма.
Метод Рунге-Кутта третьего порядка РК3
Порядок погрешности метода -
4
h
Формулы метода:
2
,
,
2
,
2
/
,
,
,
4
*
6 2
1 3
1 2
1 3
2 1
1
hk
hk
y
h
x
f
k
hk
y
h
x
f
k
y
x
f
k
k
k
k
h
y
y
m
m
m
m
m
m
m
m
(6)
Метод Рунге-Кутта четвертого порядка РК4
Наиболее распространен и имеет порядок погрешности h
5
. Формулы метода:
).
,
(
),
2
/
,
2
/
(
),
2
/
,
2
/
(
),
,
(
,
2 2
*
6
/
3 4
2 3
1 2
1 4
3 2
1 1
hk
y
h
x
f
k
hk
y
h
x
f
k
hk
y
h
x
f
k
y
x
f
k
k
k
k
k
h
y
y
m
m
m
m
m
m
m
m
m
m
(7)
Все вышеприведенные методы называются одношаговыми так как для вычисления
1
m
y
достаточно знать лишь
m
y
- значение решения на предыдущем шаге
Это позволяет использовать их при переменном шаге интегрирования
76
38.
Интегрирование систем ОДУ и ОДУ высших порядков. Пример
алгоритма.
Прикладные задачи часто приводят к системам ОДУ и к ОДУ n -го порядка В нормальной форме система ОДУ n -го порядка имеет вид
),
,...,
,
,
(
2 1
'
n
i
i
y
y
y
x
f
y
,
...,
,
2
,
1
n
i
(9)
где
n
y
y
y
...,
,
,
2 1
- неизвестные функции от переменного
x , а
i
f
- заданные функции от
1
n
переменных
Рис.8.5 – алгоритм численного интегрирования системы ДУ методом РК3
начало
)
(
),
(
,
,
,
0 0
0
x
z
x
y
h
xk
x
)
(
1
),
(
0
,
0 0
0
x
z
y
x
y
y
x
x
1
,
0
,
y
y
x
h
x
l
l
l
h
y
y
k
k
k
h
y
y
hl
hl
y
hk
hk
y
hl
hl
y
l
hl
hl
y
hk
hk
y
hk
hk
y
k
l
h
y
k
h
y
l
h
y
l
l
h
y
k
h
y
k
h
y
k
y
y
y
l
y
y
y
k
),
3 2
4 1
(
6 1
1
),
3 2
4 1
(
6 0
0
),
2 2
1 1
2 2
1 0
/(
)
2 2
1 1
(
3
),
2 2
1 1
2 2
1 0
/(
)
2 2
1 0
(
3
),
1 2
1 1
2 0
/(
)
1 2
1
(
2
),
1 2
1 1
2 0
/(
)
1 2
0
(
2
),
1 0
/(
1 1
),
1 0
/(
0 1
2
/
h
xk
x
да
конец
77
Задача Коши для системы (9) состоит в отыскании решения удовлетворяющего начальным условиям
),
(
0
t
y
i
n
i
,...,
2
,
1
.
ОДУ n -го порядка разрешают относительно старшей производной
)
,...,
,
,
,
(
)
1
(
)
(
n
n
y
y
y
y
x
f
y
(10) и введением новых переменных по правилу
,
)
(
k
k
y
y
1
...,
,
1
,
0
n
k
приводят к нормальной системе ОДУ первого порядка
)
...,
,
,
(
,
,
,
1 1
0
'
1 1
'
2 2
'
1 1
'
0
n
n
n
n
y
y
y
x
f
y
y
y
y
y
y
y
(11)
с начальными условиями
)
(
)
(
...,
,
,
)
(
)
(
)
1
(
0 0
)
1
(
0 1
0 0
0 1
0 0
0 0
n
n
n
y
x
y
x
y
y
x
y
x
y
y
x
y
x
y
(12)
Решением ОДУ(10) является n раз дифференцируемая функция
,
),
(
x
x
x
y
которая обращает уравнение (10) в тождество и удовлетворяет начальным условиям
(12)
Например для решения системы ОДУ
)
,
,
(
),
,
,
(
z
y
x
g
z
z
y
x
f
y
с начальными условиями
0 0
0 0
0
)
(
,
)
(
,
z
x
z
y
x
y
x
x
формулы метода РК4 запишутся в виде
78
,
,
,
,
,
,
2
/
,
2
/
,
2
/
,
2
/
,
2
/
,
2
/
,
2
/
,
2
/
,
2
/
,
2
/
,
2
/
,
2
/
,
,
,
,
,
,
,
2 2
*
)
6
/
(
,
2 2
*
)
6
/
(
3 3
4 3
3 4
2 2
3 2
2 3
1 1
2 1
1 2
1 1
4 3
2 1
1 4
3 2
1 1
hl
z
hk
y
h
x
g
l
hl
z
hk
y
h
x
f
k
hl
z
hk
y
h
x
g
l
hl
z
hk
y
h
x
f
k
hl
z
hk
y
h
x
g
l
hl
z
hk
y
h
x
f
k
z
y
x
g
l
z
y
x
f
k
l
l
l
l
h
z
z
k
k
k
k
h
y
y
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
(13)
39.
Методы прогноза и коррекции (перечислить). Иллюстрация простого
метода прогноза и коррекции.
Методы прогноза и коррекции:
•
Метод Хемминга
•
Первый вариант метода Адамса
•
Второй вариант метода Адамса
•
Метод на основе методов Милна и Адамса-Башфорта
Для численного решения ОДУ
)
,
(
'
y
x
f
y
используют также многошаговые методы в которых вычисление
1
m
y
ведется не только по
,
,
,
m
m
m
y
x
f
y
но и по значениям
i
m
i
m
i
m
y
x
f
y
,
,
в нескольких предыдущих узлах Формулы k - шагового метода имеют вид
,
,
0 1
1 1
1 1
k
q
q
m
q
m
q
k
q
q
m
q
m
y
x
f
b
h
y
a
y
где a
-q
, b
-q
- постоянные коэффициенты Если
0 0
b
соответствующий метод называется экстраполяционным или явным если
0 0
b
- интерполяционным или неявным методом
Частным случаем многошаговых методов является метод Адамса:
,
0 1
1 1
k
q
q
m
q
m
q
m
m
y
x
f
b
h
y
y
79
Обычно вычисления ведут по паре формул одна из которых явная а другая - неявная Такие пары формул называются методами прогноза и коррекции Прогноз выполняемый один раз на шаге служит цели получения хорошего начального приближения для последующей коррекции Последняя может выполняться на каждом шаге заданное число раз (часто только один раз) или повторяться до сходимости
Применение многошаговых методов возможно лишь в том случае если известны решения в k первых узлах Для нахождения этих значений обычно пользуются одношаговыми методами что увеличивает объём программы Поэтому в настоящее время многошаговые методы употребляются значительно реже чем четырехточечный метод Рунге-Кутта РК4
Рассмотрим простейший вариант метода прогноза и коррекции.
Формулы метода:
,
,
*
2
,
,
2
)
1
(
1 1
)
(
1 1
)
0
(
1
k
m
m
m
m
m
k
m
m
m
m
m
y
x
f
y
x
f
h
y
y
y
x
hf
y
y
(14)
По первой формуле выполняется прогноз, по второй - k -я коррекция. На рис.8.4 представлена иллюстрация этого метода при k =1, т.е. для первой коррекции.
Коррекцию можно выполнять сколько угодно раз; для получения решения итерационный процесс должен быть сходящимся. Условие сходимости:
y
f
Z
где
Z
h
,
/
2
. Вычисления прекращаются, когда
,
1 1
)
(
1
i
m
i
m
y
y
где
- заданная ошибка.
80
40.
Оценки ошибок ограничения в методе прогноза и коррекции.
Уточнение решения.
Ошибки ограничения прогноза и коррекции
Представим решение в виде разложения в ряд:
m
m
m
m
m
m
m
x
x
y
x
x
y
x
x
y
x
x
y
x
y
,
),
(
6
)
(
2
)
(
)
(
)
(
''
'
3
''
2
'
При
h
x
x
x
m
m
1
:
1 1
1
''
'
3
''
2
'
1
,
),
(
6 2
m
m
m
m
m
m
x
x
y
h
y
h
hy
y
y
При
h
x
x
x
m
m
1
:
m
m
m
m
m
m
x
x
y
h
y
h
hy
y
y
,
),
(
6 2
1 2
2
''
'
3
''
2
'
1
Из разности
1 1
m
m
y
y
получим
Рис.8.4 – иллюстрация метода прогноза и коррекции
y
x
0
y
0
x
1
x
1
y
h
решение
вестное
неиз
но
точное
x
y
,
)
(
)
0
(
2 2
)
(
y
x
y
прогноза
чения
ограни
ошибка
1
l
)
(
||
1 2
прогноз
l
l
3
l
)
1
(
2 2
)
(
y
x
y
коррекции
раничения
ог
ошибка
4
l
)
(
||
4
*
коррекция
l
l
2
x
h
)
,
(
)
0
(
2 2
y
x
)
,
(
)
1
(
2 2
y
x
81
1 1
3 3
''
'
3
'
1 1
,
),
(
3 2
m
m
m
m
m
x
x
y
h
hy
y
y
Отсюда ошибка ограничения прогноза:
1 1
''
'
3
)
(
,
),
(
3
m
m
n
T
x
x
y
h
E
Подобным же образом можно получить ошибку ограничения коррекции:
1 1
''
'
3
)
(
,
),
(
12
m
m
к
T
x
x
y
h
E
Получение оценки ошибки ограничения в процессе получения решения
Пусть
m
Y
- точное значение решения при
m
x
x
. Тогда
).
(
12
),
(
3
''
'
3
)
(
''
'
3
)
0
(
y
h
y
Y
y
h
y
Y
i
m
m
m
m
Вычтя из одного выражения другое и учитывая, что в промежутке
1 1
,
m
m
x
x
третью производную можно считать практически постоянной, получим ''
'
3
)
(
)
0
(
12 5
y
h
y
y
i
m
m
, откуда
)
(
5 1
12
)
(
)
0
(
''
'
3
)
(
i
m
m
к
T
y
y
y
h
E
Решение на шаге получается в результате i -той коррекции, и здесь же как побочный продукт вычислений можно получить оценку ошибки ограничения коррекции по выражению
)
(
5 1
)
(
)
0
(
)
(
i
m
m
к
T
y
y
E
, которая позволяет уточнить решение:
)
(
5 1
)
(
)
0
(
)
(
i
m
m
i
m
m
y
y
y
y
В заключение приведем формулы для ряда методов прогноза и коррекции (здесь
x - функция,
t
- аргумент).
41.
Достоинства и недостатки методов прогноза и коррекции.
Интегрирование ОДУ и систем ОДУ в Matchad’е.
Главное преимущество итерационных методов заключается в том, что при наличии сходимости они позволяют получить высокую гарантированную точность решения при простой расчетной схеме метода. При этом необходимо учитывать, что
82 при малых значениях e в суммарной погрешности метода будет возрастать доля погрешности, вносимой округлениями.
Скачать 1.85 Mb.