Для того чтобы численное решение (9.7) было устойчивым, необходимо, чтобы выполнялись неравенства jl-hXj\>l (;=1,..., п). Очевидно, что эти неравенства выполняются при любом положительном шаге h, если исходная система (9.3) устойчива. Таким образом, в неявном методе Эйлера соображения устойчивости численного решения не накладывают никаких ограничений на шаг интегрирования.

Анализ устойчивости методов интегрирования принято проводить на линейном уравнении

х=7,х, x{to)=Xo.

Из численного решения этого уравнения можно получить функцию роста R{hX) =Xi+i/xi. Область устойчивости определяется неравенством \R{hX)

Если область устойчивости включает всю левую полуплоскость комплексной плоскости hX, метод называется А-устойчивым. Неявный метод Эйлера является примером А-устойчивого метода. А-устойчивость является желательным свойством метода решени» жестких систем ОДУ. Однако даже А-устойчивые методы не всегда обеспечивают хорошую точность при решении жестких систем. Действительно, рассмотрим А-устойчивый метод второго порядка

Xi+i=xi-f 0,5/i[f (xi) -f f (Xi+i) ]. (9.8)

Функция роста для этого метода R(hX) = {l-\-0,5hX)/(1-0,5hX). Пусть I отрицательно, a /гЯ,>1. Тогда R(hl)!v - 1. Мы получили колебательное медленно затухающее решение, хотя на самом деле решение должно быстро затухать. Для того чтобы численное решение быстро затухало, требуется выполнение условия

J?(oo)=0. (9.9>

А-устойчивый метод, для которого выполняется условие (9.9), называется L-устойчивым. Неявный метод Эйлера является L-устойчивым, а метод (9.8) не является L-устойчивым. Как правило, для решения жестких систем ОДУ применяют неявные L-yc-тойчивые методы.

Методы численного интегрирования принято подразделять на многошаговые, использующие информацию, полученную на предыдущих шагах интегрирования, и одношаговые, использующие только информацию, полученную на последнем шаге. Многошаговые методы рекомендуется использовать в тех случаях, когда правые части ОДУ достаточно гладкие и решение необходимо получить с высокой точностью. В то же время многошаговые методы критичны к изменению шага: частые изменения могут привести к ухудшению точности и даже потере устойчивости. Поэтому в тех случаях, когда правые части могут резко изменяться и необходимы частые изменения шага, следует использовать одно-шаговые методы. Одношаговые методы имеют лучшие характеристики устойчивости, чем многошаговые, что видно из сравнения их областей устойчивости [68].

Из одношаговых методов наибольшее распространение получили методы Рунге-Кутта (теория методов Рунге - Кутта подробно изложена в [78, 89]). Общая форма представления методов Рунге - Кутта:

Xj+i = Xi-f- 2 bjkj, (9.10)

где векторы kj определяются из уравнений

kj=hi(xi+ i р,кЛ,/=1, ...,q. (9.11)

Такой метод называют (/-стадийным или (/-точечным. Если рл= = 0 для всех , то получаем обычные, явные методы Рунге - Кутта. Вычисление векторов kj по формулам (9.11) в явных методах производится последовательной подстановкой без решения алгебраических уравнений. Если хотя бы один из коэффициентов pji (Ij) не равен нулю, метод является неявным.

Можно показать, что любой явный метод Рунге - Кутта имеет ограниченную область устойчивости и не может быть А-устой-чивым. Поэтому для интегрирования жестких систем целесообразно использовать неявные методы, которые могут быть А- и L-устойчивыми. При реализации неявного метода Рунге-Кутта .необходимо на каждом шаге интегрирования решать систему нелинейных алгебраических уравнений, в общем случае порядка дп. -Практический интерес среди неявных методов Рунге - Кутта представляют полунеявные методы [78, 91] (иногда их называют полуявными или диагональными), в которых pj!=0, если />/ и система алгебраических уравнений (9.11) порядка дп распадается на д последовательно решаемых систем порядка п:

kjM (xi+ 2 hiki \ i=\ 7

/=1.

Для решения алгебраических уравнений в неявных методах обычно используют метод Ньютона или его модификации. В от-дичие от явных методов, которые однозначно определяются формулами (9.10), (9.11), неявные методы могут быть реализованы различными способами, отличающимися выбором начальных при-•ближений и числом итераций при решении алгебраических уравнений.

Для уменьшения времени счета при реализации неявного метода было предложено ограничить решение алгебраических уравнений одной ньютоновской итерацией. Такие методы получили название линейно-неявных [78]. Как и неявные методы, они также могут быть А- и L-устойчивыми. При реализации линейно-неявного метода необходимо на одном шаге интегрирования выполнить одно вычисление матрицы Якоби \=df{x)/dx и одно обращение (или LU-разложение) матрицы I-ah\, где а - параметр метода.

Применительно к полунеявным методам Рунге - Кутта методы такого типа были предложены Розенброком. Вычислительные формулы методов Розенброка имеют вид

(9.12)

х,+1 -х,- +

2 bjkj, /=1

где к;-=-(1 -a/iAj)- /if(Xj +

+ 2 Рл к,). / = 1. А, = {di (x)/ax)Ux. .

Принимая во внимание то, что автоматические системы часто содержат элементы, которые имеют разрывные характеристики, •ограничимся рассмотрением одношаговых методов интегрирования невысокого (не выше второго) порядка. Выбор методов невысокого порядка объясняется тем, что при умеренных требованиях к точности и возможных разрывах правых частей ОДУ они оказываются более эффективными. Кроме того, при совместном использовании явного и неявного методов необходимо, чтобы промежуточные точки, в которых вычисляются правые части ОДУ, •совпадали для обоих методов. Это накладывает дополнительные •ограничения в виде равенств на параметры методов, в результате задача выбора этих параметров для методов высокого порядка является весьма трудной.

Используемые ниже при построении алгоритмов моделирования методы Рунге - Кутта (явные, неявные и линейно-неявные) представлены в табл. 9.1. Явный и неявный методы первого порядка представляют собой не что иное, как явный и неявный методы Эйлера. Методы второго порядка реализуются в две стадии с одной промежуточной точкой, равной ti+ph для явного метода и ti + ah для неявного и линейно-неявного методов. Значения вектора переменных состояния в этих точках обозначены соответственно х,- 1-р и Xj-+a. Явные методы второго порядка образуют од-

Таблица 9.1