§ 5.3. Численные методы и алгоритмы моделирования

Поскольку мы строим процессы моделирования на основе ре­шения задачи Коши, будем исходить из уравнения (3.10), пере­писав его здесь для простоты в виде

(5.24)

Далее достаточно будет рассмотреть скалярное уравнение

(5.25)

Численным решением уравнения (5.25) является таблица, ко­торая на каждом шаге h дает числовые значения переменной у’ (k, h) достаточно близкие к истинным значениям y(k, h).

Классические методы численного интегрирования обыкновен­ных дифференциальных уравнений восходят к Л. Эйлеру (рос­сийский математик и механик, 1707—1783), Д. Адамсу (англий­ский математик, 1812—1892), К. Рунге (немецкий математик, 1856—1927), А. Н. Крылову (математик и инженер, член АН СССР, 1863—1945). Появление ЦВМ и бурное их развитие для ре­шения широкого круга задач науки и техники потребовало пере­оценки классических методов численного интегрирования, их ма­шинной ориентации и разработки новых, собственно машинных ме­тодов, учитывающих специфику ЦВМ. Машинная ориентация клас­сических численных методов состоит в распространении этих ме­тодов на системы уравнений большой размерности, автоматизации операций реализации методов на ЭВМ, выбору наиболее подхо­дящих методов и алгоритмов для решения конкретной задачи и величины шага h в этих алгоритмах.

Численные методы, основанные на классическом разложении правых частей (5.25) в степенные ряды Тейлора в аналитической форме, обычно исключаются из практического применения ввиду необходимости вычислений производных функций f(y,-t) порядка л>1. Однако применение ППП «Аналитические преобразования» (см. § 4.2) позволяет для ряда ММ САУ «заготавливать» эти про­изводные в аналитическом виде, что принципиально меняет воз­можности применения численных методов на степенных разложе­ниях, делает их эффективными. В соответствии с отмеченными об­стоятельствами рассмотрим построение таких методов. Будем ис­кать решения (5.25) в форме ряда

сходящегося на промежутке 0<=t<=R. Тогда все дело сводится к на­хождению неопределенных коэффициентов a_k, k принадлежит множеству от 0, до бесконечности, в (5.26). Для этого подставим (5.26) в (5.25), предварительно разложив f(y, t) в ряд Тейлора относительно точки t=0, у = уо:

Приравнивая в этом разложении коэффициенты при одинаковых степенях t, найдем

(5.27)

Таким образом, можно последовательно найти а_k и тем самым построить решение (5.26).

Однако этот алгоритм обладает существенными недостатка­ми— необходимостью ограничиваться только промежутком сходи­мости ряда (5.26), оперировать с большим числом членов этого ряда и брать производные высокого порядка от правой части в (5.25). От первого недостатка мало избавиться воспользовавшись аналитическим продолжением решения y(t)

где t_j— переменный центр разложения, a (t—tj)<R, j=0, 1, 2, ..., t₀ = 0 (рис. 5.2, а).

В этом случае коэффициенты разложения a_k каждый раз пе­ресчитываются для нового центра разложения //, у$.

От второго недостатка избавляются применением асимптотиче­ского представления у (t):

Для асимптотического представления y(t) удобно воспользо­ваться формулой Тейлора в точке t = tj:

(5.28)

где - остаточный член, представляющий малую величину вышеn-ого порядка; y^(t)(t_j) – i-ая производная от

В теории рядов Тейлора доказано, что если для какой-либо функции (t), имеющей производные n-ого порядка включительно, в точке t=h выполняется условие

(5.29)

то ии представляет собой бесконечную малую порядка малости вышеn. Тогда, обозначив в (5.28)

и обеспечив выполнение условия (5.29) выбором п и h, получим

На промежутке асимптотическое разложениеy(t) примет вид

(5.30)

В этом выражении есть два регулируемых параметра: п и h . Выбором этих параметров можно добиваться точности порядка h^n+l при числе членов разложения п. Имеющиеся в составе САПР

САУ средства программирования позволяют автоматически полу­чать аналитические выражения для y⁽ⁱ⁾(t) в (5.30). Поэтому про­цедуру численного интегрирования уравнения (5.25) можно по­строить в соответствии с алгоритмом (5.27), возложив на ЭВМ оп­ределение производных, расчет коэффициентов а_i выбор шага h<=R. Распространение алгоритма (5.27), составленного для ска­лярного уравнения на систему уравнений (3.10), состоит в вычис­лении производных и в составлении в аналитическом виде коэф­фициентов для векторного уравнения

(5.31)

где dF/dt — вектор, полученный покомпонентным частным диффе­ренцированием, a dF/dY есть (n1xn1) матрица Якоби функции F(Y, t) системы (3.10), в которой для простоты положено X=λ=U=0.

Построение численных методов, в которых переход от y’(t_j) к y’(t_j+1) не требовал бы вычисления производных от правых час­тей в (5.25), основано на идее Рунге. Принцип построения таких формул численного интегрирования можно усмотреть из геомет­рических построений (рис. 5.2, б, в). Продолжая этот процесс и организуя комбинации k_i(h) для нескольких точек внутри проме­жутка t_j+1—t_j, запишем общую формулу этой группы численных методов:

(5.32)

где

Из (5.32) при различных значениях _i, _i получим как частные случаи известные формулы численного интегрирования Рунге и Кутта (венгерский математик, 1867-1944) разного порядка точности:

первого порядка

(5.33)

- метод Эйлера;

второго порядка

(5.34)

- метод Эйлера – Коши;

(5.35)

- усовершенствованный метод Эйлера;

четвертого порядка

Точность представления у на одном шаге в каждой формуле оценивается остаточным членом разложения функции φ(h) в ряд Тейлора:

Величина r в (5.37) называется порядком точности метода или степенью метода. Важно сразу же отметить, что выбор шага h ис­ходя из заданой точности на одном шаге в принципе не гаранти­рует требуемой точности конечного результата.

Следующая группа методов численного интегрирова­ния, основанная на идеях Адамса, строится на основе аппрокси­мации подынтегральной функции в формуле Ньютона—Лейбница

известным полиномом. Эта формула получена после интегрирова­ния уравнения (5.25) на промежутке t_j, t_j+1. Если в качестве тако­го полинома выбирать полином Ньютона, то у’ (t) можно предста­вить в виде

Задавая предыдущие значения y’(t_j), y’(t_j-1), y’(t_j-2)...в предпо­ложении, что они вычислены заранее, вычисляем неизвестные ко­эффициенты:

(5.40)

где . Подставляя (5.40) в (5.38), после интегрирования получим

или в общем виде

Методы типа (5.41) называют также многошаговыми в отличие от одношаговых (5.32), так как для получения y_j+1 необходимо предварительно вычислить п предыдущих значений y_j-k ,k = 0, ..n. Очевидно, недостатком методов этого типа является необходимость предварительного вычисления «разгонных» значений y’(t_j-1). Методы (5.32), (5.41) реализованы в виде стандартного програм­много обеспечения современных ЭВМ.

Применимость методов (5.32), (5.41) к моделированию САУ тесно связана с выбором шага h в формулах численного интегри­рования. В случае одного уравнения (5.25) значение h выбирается из соображений точности на одном шаге и оно же обеспечивает устойчивость вычислительного процесса. Когда же рассматривает­ся система уравнений (3.10), то именно устойчивость разностных уравнений (5.32), (5.41) ограничивает шаг. Наибольшие шаги, обеспечивающие устойчивость разностной схемы, или, как их обыч­но называют, критические шаги, различаются по величине для раз­ных методов. Так, например, если собственные значения вычислен­ной в данной точке t=t* матрицы Якоби системы (3.10) вещественны, то критический шаг выбирается из соотношения

(5.42)

где λ_mах — максимальное собственное значение; r — число, завися­щее от метода.

Для метода Рунге — Кутты четвертой степени оно равно 2.78, для методов Эйлера и Эйлера — Коши — 2, для ме­тода Адамса четвертой степени — 0,3. Причем даже незначительное увеличение значения h против допустимого приводит к потере вы­числительной устойчивости — «взрыву погрешности».

Подобные явления возникают при численном решении задачи Коши не только для САУ, но и для других объектов. Это обстоя­тельство повлияло на усиленное внимание к разработке специ­ально машинных методов численного интегри­рования, учитывающих особенности ЭВМ.

Одна из основных идей таких методов заключается в построе­нии формул численного интегрирования с «обратной связью» — так называемых неявных методов. В неявных методах можно дос­тичь увеличения шага h по сравнению с явными без потери устой­чивости вычислительного процесса. В простейшем случае неявный метод получается из формулы (5.38), если воспользоваться фор­мулой правых прямоугольников при взятии интегралов в (5.25):

Аналогично строятся неявные формулы для методов (5.32), (5.41). Например, неявная формула второго порядка будет иметь вид

Методы этого класса являются устойчивыми и не имеют ограничености может выбираться сообразно с требуемой длительностью ре­ний вида (5.42) на величину шага дискретности. Шаг дискретно шения.

При реализации неявных методов необходимо на каждом шаге решать нелинейные уравнения типа (5.45), что вызывает значи­тельные осложнения и является недостатком этих методов. Поэто­му применение неявных методов сочетают со специальными мето­дами решения нелинейных уравнений на каждом шаге.

Так, метод Эйлера (5.45) с модификацией его по формуле Нью­тона примет вид

Из других неявных методов приведем методы, предложен­ные в 1952 г. американскими математиками С. Ф. Куртисом и Н. О. Хиршфельдером, так как их применение оказалось эффек­тивным при моделировании САУ. Неявные алгоритмы Куртиса — Хиршфельдера строятся на принципе «дифференцирования назад» и имеют вид

т. е. в выполнении своеобразного «дифференцирования назад».

Применение неявных методов при моделировании CAY показы­вает их преимущества над явными в смысле устойчивости вычис­лительного процесса при увеличении шага, однако точность полу­чения решений существенно зависит от точности решения нелиней­ных алгебраических уравнений (5.45), (5.49). Причем если в яв­ных методах выбором шага по (5.44) удается получить, хотя и с погрешностью, решение (5.25) качественно верное, то в неявных при неправильном выборе шага можно получить устойчивые ре­шения, но качественно отличные от искомого. Значительные уве­личения шага при моделировании САУ неявными методами ока­зываются ограниченными из-за резкого снижения точности. Прак­тически по сравнению с явными методами удается увеличить шаг не более чем в 5—10 раз.

Эффективное применение неявные методы и алгоритмы нашли в процедуре, предложенной Гиром, с изменениями и дополнения­ми, внесенными Хиндмаршем (С. В. Гир и А. С. Хиндмарш — аме­риканские математики). В процедуре Гира применяются методы (5.49). Однако это не исключает возможности применения и дру­гих численных методов. Неявный характер употребляемых в про­цедуре Гира методов предполагает решение на каждом шаге сис­темы алгебраических уравнений, для чего предусматривается ис­пользование трех модификаций итерационного метода Ньютона.

Для этого каждое из уравнений (5.49) представляется в виде

где α_i, β_i— числовые коэффициенты; К₁, K₂ — количество соответ­ствующих слагаемых в (5.49).

Pешение нелинейных алгебраических уравнений (5.50) относи­тельно y_j+1 строится как следующий процесс:

В процедуре Гира автоматически выбираются численный метод и шаг по заданной точности вычислений.

Если в некоторый момент времени используется метод степени l (например, в (5.51) К₁ = l, К₂=l—1), то одновременно оценива­ется локальная погрешность методов степени l—1, l+1. В дальней­шем применяется тот метод, который при заданной точности поз­воляет интегрировать с большим значением шага h.

Еще одно преимущество процедуры Гира в том, что она не тре­бует стартового алгоритма. Интегрирование начинается с методов первой степени, а затем степень метода при необходимости изме­няется в соответствии с задаваемой пользователем точностью. Сле­дует отметить, что алгоритмы процедуры подвержены неустойчи­вости, которая проявляется в колебательном характере шага, при­нимающем то разумные значения, то неприемлемо малые. Возни­кающие трудности можно исключить уточняя соответствующие константы α_j, β_j эмпирическим путем.

Все рассмотренные методы являются скалярными, т. е. для каж­дого уравнения системы (3.10) составляется соответствующее раз­ностное уравнение, поэтому они обладают присущими скалярным методам недостатками при применении их к системам большой размерности. Качественно новым является создание системных (матричных) методов, которые позволяют учитывать свойства пра­вой части системы (3.10) и обеспечивают независимость выбора шага от устойчивости вычислительного процесса.

Системные методы были предложены Ю. В. Ракитским (со­ветский инженер и математик, 1932—1986). Эти методы построены на введении в разностное уравнение, эквивалентное исходному (3.10), матричного множителя Ф(А, H) который «подобен» шагу интегрирования h в методе Эйлера (5.33):

при этом Hи h связаны соотношением H=h*2^N, где h — началь­ный шаг интегрирования, а точность вычислений Ф₀ соответствует точности метода (5.34) при интегрировании системы Y=AY. На рис. 5.4, а приведена последовательность операций по реализации метода (5.52).

Верхняя граница для h удовлетворяет соотношению h <=1/[[A]] а достаточное порочности нижнее значение h определяется срав­нением Y(k, h), Y(kh/r), вычисленных при h и h/r.

Рассмотрим теперь системные методы для решения систем ли­нейных уравнений, получаемых из общего уравнения (3.10) после

где под √А понимается матрица, для которой (√A)²=A.

Формулы (5.59) — (5.65) определяют Y(t) для любого момента времени. Казалось бы, задача решена и с точки зрения машинной реализации сведена к наиболее простой операции — «счету по фор­мулам». Однако численная реализация этих формул упирается в серьезную проблему вычисления функциональных рядов (5.62) — (5.65). Непосредственное применение формул (5.62) — (5.65) для этой цели при достаточно больших t и плохой обусловленности матрицы А приводит к необходимости учета большого числа чле­нов, увеличивающихся с ростом t, поэтому для организации вычис­лений от (5.56), (5.57) осуществляется переход к эквивалентным разностным уравнениям, по которым и осуществляется вычисление решения шаг за шагом. Такой переход уже был сделан для систе­мы (3.15) (см. § 3.3, формула (3.26)).

Для систем (5.56), (5.57) эквивалентные разностные уравнения имеют вид

где Н — шаг дискретности, выбираемый из удобства наблюдения процесса Y(t) и не зависящий от свойств этого процесса.

По существу, в (5.66), (5.67) все дело сводится к вычислению матриц е^АН и вектора

Существует много способов вычисления матрицы е^АH , что указы­вает на сложность этой процедуры и зависимость ее от свойств матрицы А и величины Н. При малых шагах h матрица е^АН вы­числяется значительно проще. В практике численного интегриро­вания линейных уравнений, описывающих движение САУ, матри­ца е^Аh вычислялась двумя способами: по формуле (5.62) при t=h; путем решения вспомогательного матричного уравнения U = AU, U(0)=E с последовательным формированием е^Ah по столб­цам

Матрица Ф₀ вычисляется по формуле (5.55). На рис. 5.4, б при­ведена последовательность операций по реализации алгоритма (5.69).

Выбор методов и алгоритмов численного интегрирования, под­ходящих для решения конкретной задачи (3.10) или (5.24), можно проверить исходя из свойств правых частей уравнений (5.24). Для моделирования САУ с числом обусловленности ρ(Я)<10³ и вре­менем решения t_Н<<1 мин вполне приемлемым оказываются ска­лярные методы (5.32), (5.49).

В случае моделирования САУ с числом обусловленности их ММ ρ(Я)<10⁶ t_н<5 мин вполне приемлемыми оказываются не­явные методы (5.47), (5.41) с их реализацией в виде процедуры, аналогичной процедуре Гира (см. § 5.3).

Моделирование САУ с числом обусловленности р(Я)>Ю^б и временем наблюдения /_н>5 мин скалярными методами становит­ся малоэффективным, и необходимо применять системные методы (5.52), (5.60) — (5.69).

Практика моделирования САУ с помощью системных методов показывает, что эти методы позволяют существенно — в 20— 30 раз — уменьшить время получения решения по сравнению со скалярными. Это дает возможность эффективно решать задачи мо­делирования САУ на ЭВМ по полным их моделям. Недостатком системных методов (5.52) является необходимость вычисления матрицы Ф(А, Н) в каждом частном случае и составления специ­альной программы пересчета шагов.