Афанасьев А.А. Математические основы теории систем управления
Подождите немного. Документ загружается.
Таблица 4.1
t
0
t ht +
0
ht 2
0
+
K
hnt
+
0
K
Tt
Nht
+=
=+
0
0
x
)(
0
tx ) )2(
~
0
htx(
~
0
htx + +
K
)(
~
0
nhtx +
K
)(
~
0
Ttx +
Здесь h > 0 - некоторая малая величина, называемая шагом интегрирования
или шагом дискретности таблицы. Будем обозначать через
),,1,0,(,)(
~~
0
h
T
Nnhntttxx
nnn
==+==
K
сеточную (решетчатую) функцию.
Для вычисления значений )(
~
n
tx
(
nn
txx
, приближающихся с необходимой
точностью к истинным значениям )= , дифференциальное уравнение (4.4)
заменяют алгебраическим, называемым разностным уравнением. Разностное
уравнение оперирует с величинами, фиксированными на дискретной
равномерной сетке . Оно может иметь такой вид [15]:
n
tt=
)(
~~~
110110
nnmnmnn
bfbfbhxaxaxa
−−− mnm
f
−
+++=+++
KK
, (4.5)
где
)
~
,(
nnn
xtff =
,
h
T
Nnm =≤≤ ; числовые коэффициенты ),,1,0(, mkba
kk
K
=
не зависят от n, причем
.0
0
≠a
Уравнение (4.5) называют методом численного интегрирования. Его
следует рассматривать как рекуррентное соотношение, выражающее новое
значение через найденные ранее значения .
n
x
~
mnnn
xxx
−−−
~
,,
~
,
~
21
K
Расчет начинается со значения n=m, т.е с уравнения
)(
~~
0000
fbfbhxaxa
mmmm
++=++
KK
.
Отсюда видно, что для начала расчета необходимо задать m начальных
значений
.
~
,,
~
,
~
,
1210
−m
xxxx
K
Эту совокупность называют началом таблицы , а способ вычисления значений
- стартовым алгоритмом.
11
~
,
~
−
m
xx
K
При m=1 метод численого интегрирования называют одношаговым, а
при m > 1- многошаговым.
121
Построение разностного уравнения разностной схемы (4.5) может
производится на основе исходного уравнения (4.4) и формулы Ньютона-
Лейбница
(4.6)
∫
+
+=
+
1
))(,()()(
1
n
n
t
t
nn
dxftxtx
τττ
с последующей аппроксимацией интеграла.
Явный метод (ломаных) Эйлера получают при приближенном
вычислении интеграла по способу правых треугольников (рис.4.1, а) .
При этом уравнение (4.6) принимает вид
001
~
;,0,1, = ,
~~
xxNnfhxx
nnn
=+=
+
K
, (4.7)
где
}.
~
,{
nnn
xtff =
К неявному методу (ломаных) Эйлера приводит приближенное
вычисление интеграла по способу левых треугольников (рис 4.1,б)
11
~~
++
+=
nnn
fhxx
, (4.8)
где
)
~
,(
111
+++
=
nnn
xtff
.
а
б
Рис.4.1
Видим, что в этом уравнении присутствует как в левой, так и в правой
части, т.e. уравнение (4.8) не разрешено явно относительно . Неизвестное
1
~
+n
x
1
~
+n
x
122
1
~
+n
x
находят, решая уравнение (4.8) итерационным методом Ньютона, выбирая
начальное приближение
)0(
1
~
+n
x равным
n
x
~
.
В общем случае метод (4.5) называют явным, если 0, и, следовательно,
искомое значение
0
=b
n
x
~
выражается явным образом через предыдущие значения
.
mnn
xx
−−
~
,,
~
21
K
n
x
−
~
,
В противном случае
)0(
0
b ≠
метод называют неявным . Неявным,
например, является метод трапеций, получаемый при аппроксимации в (4.6)
по методу трапеций
)(
2
~
n
x +
~
11
nnn
ff
h
x +=
++
. (4.9)
Методы Эйлера и трапеций относятся к одношаговым методам (m=1) .
Из многошаговых методов в практике вычислений наибольшее
распространение получили методы Адамса, в которых производная
d
t
dx
аппроксимируется только по двум точкам т.e в формуле (методе) (4.5)
имеем
,
1
−nn
tиt
.),,2(,0,1
10
mkaa
K
,1 a
K
==−==
Таким образом, методы Адамса имеют вид:
kn
m
k
k
n −n
fb
h
xx
−
=
∑
=
~ ~
−
0
1
или
~
)(
~
110
mnmnnn
fbfbfbhx
−+
1
n +
1
x
+
= + + ++
K
. (4.10)
В случае методы Адамса называются явными, в случае 0
0
=b 0
0
≠b-
неявными.
При построении разностных схем наряду с формулой Ньютона-Лейбница
(4.6) используется представление решения (4.4) в виде ряда Тейлора
,
!!
1
1
11
∞
=
∑
k
k
h
)
1
−
−
∞
=
−
∑
==−
k
n
k
k
k
k
n
k
k
nn
dt
fd
k
h
dt
xd
xx (4.11)
взятого в окрестности точки
(
nn
txx =
. Здесь высшие производные находят
путем последовательного дифференцирования правой части уравнения (4.4).
Уравнение (4.11) можно переписать в разностной форме первого порядка:
123
1
1
1
1
1
!
~~
−
−
=
−
+
∑
+=
k
n
k
k
k
nn
d
t
fd
k
h
hxx
ν
.
При 1=
ν
это уравнение совпадает с уравнением явного метода Эйлера
nnn
fhxx +=
+
~~
1
.
4.2.2. Погрешность аппроксимации вычислительной схемы
В соответствии с (4.11) погрешность аппроксимации схемы Эйлера
K
&&
!3!2
~~
2
111
h
x
h
f
h
xx
h
xx
h
xx
nn
nnnnnn
n
+=−
−
=
−
−
−
=
+++
ψ
имеет первый порядок, т.е.
nn
илиh
ψ
ψ
)(0
1
=
<Mh , где M-константа.
Под погрешностью аппроксимации разностной схемы (4.5) понимают ее
невязку
)()(
)
~~
()(
00
00
mnmnmnmn
mnmnmnmnn
fbfbhxaxa
xaxaxaxa
−−
−−
++−++=
=
++−++=
KK
KK
ψ
на точном решении исходного уравнения состояния (4.4).
В общем случае говорят, что разностная схема имеет p-ый порядок
аппроксимации, если
P
n
Mh≤
ψ
и разностная схема аппроксимирует исходное дифференциальное уравнение,
если
0→→ hприo
n
ψ
.
4.2.3. Порядок точности вычислительной схемы
Известно также, что разностная схема (разностный метод) имеет p-ый
порядок точности, если
0
~
→≤− hприMhxx
P
nn
.
Можно доказать [15], что порядок точности разностного метода
совпадает с порядком аппроксимации.
124
Можно также показать [15], что 2m - неизвестных коэффициентов
(по условиям нормировки коэффициенты линейно
зависят от остальных: ) в разностном методе (4.5)
удовлетворяют системе из p -линейных уравнений. Следовательно, этот метод
обеспечивает порядок точности p=2m.
mm
bbaa ,,,,
1,1
KK
00
,ba
1;0
00
==
∑∑
==
m
k
k
m
k
k
ba
Итак, наивысший достижимый порядок точности неявных m-шаговых методов
равен 2m, а явных очевидно, 2m-1. ),0(
0
=b
Многошаговые методы Адамса (4.10) имеют наивысший порядок
точности p = m+1 для неявного варианта (имеем p уравнении при (m+1)
неизвестных коэффициентах ) и - p =m для явного варианта
m
bbb ,,,
10
K .)0(
0
=b
Для явного m-шагового метода Адамса условия m-го порядка точности
имеют вид:
pml
l
bk
k
m
k
l
===
∑
=
−
,,2,1,
1
1
1
K . (4.12)
Решая систему (4.12), можно найти коэффициенты для
конкретного m=p. Например, при m=1 (первом порядке точности:
),,2,1( mkb
k
K
=
1
==
p
m )
получаем из (4.12): Следовательно, из (4.10) имеем
.1
1
=b
.
~~
1
n
nn
f
h
xx
=
−
+
Этот алгоритм является явным методом Эйлера.
При 2 имеем из (4.12), принимая последовательно =m 2,1=
l
:
2
1
2
,1
21
21
=+
=+
bb
bb
.
Из этой системы уравнений получаем
2
1
,
2
3
21
−== bb . Вычислительный
алгоритм будет иметь вид:
1
1
2
1
2
3
~~
−
+
−=
−
nn
nn
ff
h
xx
, (4.13)
125
При m=3, 4 получаем соответственно следующие методы m-порядка
точности:
)51623(
12
1
~~
21
1
−−
+
+−=
−
nnn
nn
fff
h
xx
, (4.14)
)9375955(
24
1
~~
321
1
−−−
+
−+−=
−
nnnn
nn
ffff
h
xx
. (4.15)
Для m-шаговых неявных методов Адамса
mnmnn
nn
fbfbfb
h
xx
−++
+
+++=
−
1110
1
~~
K
наивысший порядок точности равен m+1.
При m =1 получаем метод 2-го порядка точности (p =2).
Для него
2
1
10
== bb
:
)(
2
1
~~
1
1
nn
nn
ff
h
xx
+=
−
+
+
.
Это упоминавшийся выше метод трапеций.
При m=2,3 получаем методы соответственно 3- и 4-го порядка точности
)85(
12
1
~~
11
1
−+
+
−+=
−
nnn
nn
fff
h
xx
, (4.16)
)5199(
24
1
~~
211
1
−−+
+
+−+=
−
nnnn
nn
ffff
h
xx
. (4.17)
Методы Рунге-Кутта
отличаются от разностных методов тем, что в них
допускается вычисление правых частей f(t,x) не только в точках сетки, но и
некоторых промежуточных точках. Все классические методы Рунге-Кутта
являются явными и одношаговыми [16]. Наиболее распространены на практике
схемы Рунге-Кутта второго и четвертого порядков точности.
Вычисления по схеме Рунге-Кутта второго порядка точности проводятся
в два этапа. На первом этапе находится промежуточное значение по явной
схеме Эйлера с шагом
x
n
∗
h
α
)
~
,(
~
nnnn
xtfhxx
α
+=
∗
.
На втором этапе находится значение
1
~
+n
x
по формуле
,),()
~
,()1(
~~
1
∗
+
++−+=
nnnnnn
xhthfxtfhxx
α
σ
σ
126
где
>0, σ>0 - параметры. α
Более кратко этот алгоритм запишем так
nnnn
hkhkxx
211
)1(
~~
σ
σ
+−+=
+
,
где
).
~
,(),
~
,(
21
••
==
nnnnnn
xtfkxtfk
Здесь обозначено:
htt
nn
α
+=
•
.
Исключая из последнего уравнения , получим для
∗
n
x
1
~
+n
x
схему:
))
~
,(
~
,()
~
,()1(
~~
1
nnnnnn
nn
xtfhxhtfxtf
h
xx
αασσ
+++−=
−
+
nn
kk
21
)1(
σ
σ
+−=
. (4.18)
Порядок точности зависит от параметров
α
σ
,
.
Найдем погрешность (невязку)
ψ
n
этой схемы:
h
xx
h
xx
nnnn
n
~~
11
−
−
−
=
++
ψ
.
Воспользовавшись разложением по формуле Тейлора (4.11) для
nn
xx
−
+
1
и
формулой (4.18), получим:
)()
2
1
(
22
)(
!3!2
)1()
!3!2
(
222
22
hMxhMhx
h
h
h
xx
Mhx
h
xx
f
h
f
h
fffff
h
f
h
f
nn
nn
nnn
nnnnnnnnnn
++−=++
−
−=++−
=+++−=−−−++=
•
•
••
&&&&
&&
&&&&
K
&&&
K
&&&
σαα
α
σσ
σσσσ
ψ
.
Отсюда видно, что схема (4.18) имеет второй порядок точности
2
Mh
n
≤
ψ
при
выполнении условия
2
1
=
σα
.
Рассмотрим частные случаи:
1)
1,
2
1
==
σα
(усовершенствованный метод ломаных)
)
~
,(
2
~
nnnn
xtf
h
xx +=
∗
,
nnnnnn
hkxx
h
tfhxx
21
~
),
2
(
~~
+=++=
∗
+
.
Это известная схема предиктор-корректор (или счет-пересчет). Ее можно
переписать, после исключения , в виде
∗
n
x
127
)
2
~
,
2
(;)
~
,(;
~~
1212
1
nnnnnnnn
nn
k
h
x
h
tfkxtfkk
h
xx
++===
−
+
. (4.19)
2)
α=
,
1
σ=
1
2
(метод Эйлера-Коши)
)
~
,(;)
~
,(;)(
2
1
~~
112121
1
nnnnnnnnn
nn
khxtfkxtfkkk
h
xx
+==+=
−
+
+
.
(4.20)
Приведем формулы для схемы Рунге-Кутта 4-го порядка точности
)22(
6
1
~~
4321
1
nnnn
nn
kkkk
h
xx
+++=
−
+
, (4.21)
где
),
~
,(
1 nnn
xtfk
=
)
2
~
,
2
(
12 nnnn
k
h
x
h
tfk
++=
, )
2
~
,
2
(
23 nnnn
k
h
x
h
tfk
++=
,
)
~
,(
4 nn
xhtfk
3nn
kh
++=
.
Отметим, что рассмотренные методы справедливы и в том случае, когда
(4.4) представляет собой систему дифференциальных уравнений
),(
xtf
d
t
xd
=
,
т.е когда переменные
fиx являются вектор-функциями:
.][,][
T
2121 n
T
n
ffffxxxx
KK
==
При интегрировании рассмотренными разностными схемами систем
уравнений состояния их запись в векторной форме остается такой же, как и в
случае одного уравнения, только в соответствующих выражениях вместо
скалярных величин появятся векторы и матрицы. Будем иметь:
;][
21
Е
nnnnn
xxxx
K
=
);
~
,(][
121111 nn
Е
nnnnn
xtfkkkk ==
K
128
);
2
~
,
2
(
12 nnnn
k
h
x
h
tfk
++=
)
2
~
,
2
(
23 nnnn
k
h
x
h
tfk ++=
;
)
~
,(
34 nnnn
khxhtfk
++=
.
Вопросы для самопроверки
1.
Напишите рекурсивную формулу вычисления матричной экспоненты по
методу Ю.В. Ракитского.
2.
Назовите одношаговые алгоритмы численного решения уравнений
состояния.
3.
Как выглядят алгоритмы явного и неявного метода Эйлера?
4.
Что такое порядок точности вычислительной схемы?
5.
Каков порядок точности метода трапеций?
6.
Каков наивысший достижимый порядок точности m-шаговых явных и
неявных методов?
7.
Скольки шаговыми являются методы Рунге-Кутта. Они являются явными
или неявными?
8.
Приведите формулу метода Рунге-Кутта 4-го порядка.
4.3. Устойчивость методов численного интегрирования
Выбор шага численного интегрирования уравнения состояния зависит не
только от необходимой точности его решения, но и от условий обеспечения
устойчивости соответствующей разностной схемы.
Посмотрим как изменится решение уравнения (4.4) при изменении
начального условия с на . Пусть
0
x
0
~
x )(
~
t
x
и x(t)- решения этого уравнения с
начальными условиями
0
~
)0(
~
xx
= и
0
(x )0 x
= . Для погрешности
)()(
~
)(
t
x
t
x
t
z
−=
получаем уравнение
,
~
)0(,0,)(),()
~
,(
000
xxzzTtztxtfxtf
d
t
dz
−==≤<=−=
α
(4.22)
где
),(
),()
~
,(
)( zxt
d
t
f
z
xtfxtf
t
θ
∂
α
+=
−
= , .10 ≤≤
θ
129
Решением (4.22) является экспоненциальная функция
∫
=
t
dss
eztz
0
)(
0
)(
α
.
Если
0
<
x
f
∂
∂
для всех t и x, то
[]
,,0
~
)()(
~
)0()(
00
Ttвсехдляxxtxtxилиztz
∈−<−<
т.e., решение уравнения (4.4) устойчиво по начальным данным.
Решение задачи (4.22) при
∞→
t
ведет себя аналогично решению
линейного однородного уравнения
,)0(,0,
0
zzTtz
d
t
dz
=≤<=
λ
которое можно рассматривать как модельное (тестовое) уравнение при
изучении устойчивости.
Для оценки устойчивости по начальным условиям нелинейного
уравнения (4.4) будем рассматривать модельное уравнение вида
.)0(,0,0Const,
0
xxtx
d
t
dx
=><==
λλ
(4.23)
Его решение
t
extx
λ
0
)( =
убывает при 0
<
λ
и
0
)( xtx
≤
при 0
≤
λ
для всех t >0, т.е уравнение (4.23)
устойчиво при
0≤
λ
, что соответствует условию 0
≤
x
f
∂
∂
(при уравнение
считают асимтотически устойчивым, при 0
λ<0
=
λ
- просто устойчивым, при 0
>
λ
-
неустойчивым).
Оценим устойчивость некоторых разностных схем.
1) Явная схема Эйлера.
Применим для решения модельного уравнения (4.23) явную схему
Эйлера. Имеем, учитывая, что в модельном уравнении
nn
xf
λ
=
,
130