2.1.5. Устойчивость численного решения
Выше было показано, что численное решение системы дифференциальных уравнений (1) с матрицей L собственных значений эквивалентно точному решению системы уравнений со смещенной матрицей собственных значений L+dL, определяемой выражением (19). Так как входящие в него матрицы диагональны, то легко написать соответствующие скалярные соотношения:
eZkh = e(Zk+dZk)h = ebkh , bk = Zk+dZk , k=1…n (2.1.23)
где bk – смещенные собственные значения.
Устойчивость исходной системы уравнений (1) или, что то же самое, ее точного решения, определяется следующими условиями:
Re(Zk) < 0 , k = 1…n.
При численном решении дифференциального уравнения корни смещаются на dZk . Это может привести к тому, что какой-то смещенный корень bk переместится в правую полуплоскость, что будет соответствовать неустойчивости численного решения.
Условием устойчивости численного решения является, очевидно, соотношение:
Re(bkh) = Re[(Zk +dZk)h] < 0 , k=1…n. (2.1.24)
Из (23) с учетом (18) найдем значение смещенного корня:
bkh = ln(eZkh) = ln[1+Zkh+…+(Zkh)r/r!] (2.1.25)
Представим eZkh в виде комплексного числа:
eZkh = u + iv = Neiф , (2.1.26)
N
=
=
,
ф=
arctg(v/u)
Учитывая (25) и (26) , условия устойчивости численного решения (24) можно записать:
Re(bk
) = ln
< 0 (2.1.27)
В координатах u и v из (26) областью устойчивости является круг единичного радиуса с центром в начале координат:
u2 + v2 < 1 (2.1.28)
Соотношения (27) и (28) позволяют избежать грубых ошибок при выборе шага дискретизации h и порядка
метода r.
2.1.6. Повышение точности численного решения методом коррекции уравнений движения
Приближенное численное решение X(t+h) = eAh X(t) дифференциального уравнения (1) эквивалентно точному решению Х(t+h) = e(A+dA)h X(t) измененного уравнения
,
где dA определяется из соотношения
eAh = e(A+dA)h
Это позволяет поставить задачу о предварительной коррекции исходного уравнения таким образом, чтобы численное решение скорректированного уравнения стало равным точному решению исходного уравнения. Исходное уравнение (1) имеет точное решение на шаге:
X(t+h) = eAh X(t) (2.1.29)
Введем скорректированное уравнение
(2.1.30)
Численное решение уравнения (30) на шаге можно записать так:
X(t+h) = e(A+dA)h X(t) (2.1.31)
Здесь
e(A+dA)h = [I + (A+dA)h + …(A+dA)rhr /r! ] (2.1.32)
- усеченная экспонента матрицы (A+dA) .
Приравнивая решения (29) и (31), можно получить соотношения:
eAh = e(A+dA)h ; e(A+dA)h – eAh = deAh (2.1.33)
Отсюда, используя соотношения (6) и (32)и пренебрегая старшими степенями dA , получим :
dAh [I + Ah + …+ (Ah)(r-1) /(r-1)! ] = deAh (2.1.34)
В первом приближении оценку dA получим из (5.34) с учетом (6):
dAh (Ah)(r+1) /(r+1)! ( 2.1.35)
Сравнивая (35) с (14) можно видеть, что в первом приближении
dA -dA
Уравнения (5.33)…(5.35) с разной степенью точности определяют корректирующую матрицу dA , которую нужно прибавить к матрице исходного уравнения (5.1), чтобы численное решение скорректированного уравнения (5.30) совпало с точным решением на шаге исходного уравнения (1). Диагонализируя соотношения (33)…(35), получим:
eLh = e(L+dL)h ; e(L+dL)h – eLh = deLh
dLh [I + Lh + …+ (Lh)(r-1) /(r-1)! ] = deLh (2.1.36)
dLh = (Lh)(r+1) /(r+1)!
Все выражения в (36) являются диагональными матрицами, поэтому могут быть переписаны в скалярном виде:
eZkh = e(Zk+dZk)h ; e(Zk+dZk)h – eZkh = deZkh
dZkh [I + Zkh + …+ (Zkh)(r-1) /(r-1)! ] = deZkh (2.1.37)
dZkh = (Zkh)(r+1) /(r+1)! k= 1…n
Уравнения (37) с разной степенью точности определяют
приращения dZk , на которые нужно скорректировать
собственные значений Zk матрицы А, чтобы численное решение измененного таким образом уравнения совпало с точным решением на шаге исходного уравнения (1).
Приведенные в разделе 2.1 результаты получены для решения линейных однородных дифференциальных уравнений методом разложения в ряд Тейлора. Однако, во многоих случаях, полученные оценки могут быть распростронены на:
- неоднородные уравнения [2.5]
- линеаризуемые уравнения
- другие одношаговые методы решения дифференциальных уравнений, основанные на приближении локальной ошибки к ошибке разложения в ряд Тейлора.