нопараметрическое семейство, определяемое параметром р. Оптимальное значение р, обеспечивающее минимум главного члена погрешности, 2/3. При совместном использовании явного и неявного

методов второго порядка следует задавать р = а=1-12/2, в этом случае промежуточные точки для обоих методов совпадают. Неявный метод второго порядка при выбранных параметрах является не только L-устойчивым, но и S-устойчивым [91], что весьма ценно при интегрировании жестких систем. Отметим, что форма записи неявного и линейно-неявного методов в табл. 9.1 несколько отличается от общепринятой (см. формулы (9.10) -(9.12)), однако она оказалась более удобной при алгоритмической реализации.

Для выбора шага интегрирования можно использовать следующую процедуру. Интегрирование от точки ti до точки ti+i = ti + hi выполняется дважды: с шагом hi н с шагом hil2. За текущее значение вектора переменных состояния принимается решение хг+ь

полученное с шагом hi/2. Решение Xj+i, полученное с шагом hi, используется для оценки локальной ошибки, которая вычисляется по формуле

е=(хг+,-Xi+i)/(2-l),

где k - порядок метода. Оценку относительной локальной ошибки получаем в виде

е=е1/хг+1-х,-.

Если оценка ошибки больше допустимой (8>8доп), то следует уменьшить шаг hi в 2 раза и повторить расчет из точки ti. Если е<;едой/2, то шаг следует увеличить в 2 раза, т. е. принять hi+i = = 2/Zi. И наконец, если 8доп/288доп, то шаг интегрирования остается без изменений.

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

9.2. МОДЕЛИРОВАНИЕ СИСТЕМ ПО УРАВНЕНИЯМ СОСТОЯНИЯ

Рассмотрим автоматическую систему, уравнения которой представлены в нормальной форме

x=f (X, t), X{to) =Хо.

Для моделирования такой системы непосредственно применимы численные методы интегрирования ОДУ. В приведенной ниже подпрограмме RKE использован явный метод Рунге - Кутта второго порядка из табл. 9.1 при Р = 2/3, Расчетные формулы этого метода:

Xi-i-i = Xi+ hf {Xi, ti) + hf (x,-+2/g, tt + h

4 4 \ 3

Подпрограмма RKE

Назначение; решение системы дифференциальных уравнений явным методом Рунге -Кутта с автоматическим выбором шага.

Обращение: CALL ,RKE(TO, Tl, DT, HMAX, HMIN, EPS, X, N, IER, FUN, OUT).

Параметры:

TO - начальное значение независимой переменной, Tl - конечное значение независимой переменной,

DT - шаг вывода результатов, при DT=0 вывод через каждые два шага Интегрирования, . НМАХ - максимальный шаг .интегрирования, рекомендуется HMAX=DT,. HMIN - минимальный шаг интегрировании, рекомендуется НМШ= (Т1- -Т0)>1<1.Е-5,

EPS - допустимая относительная ошибка, рекомендуется EPS=0.01, не следует задавать EPS больше 0.1 или меньше 1.Е-4,

X - начальное значение вектора переменных состояния (на входе),

конечное значение вектора переменных состояния (яа выходе),

N - размерность вектора переменных состояния,

1ER - код завершения: 0 - интегрирование завершено успешно, 1 - ошибка во входных данных, й - процесс расходится (одна яз переменных по абсолютной величине превысила 1.Е 15), 3 - заданная точность не обеспечивается, вероятно, система является жесткой, и для ее решения следует использовать неявный метод интегрирования,

FUN - имя подпрограммы вычисления производных, ее параметры: Т - значение независимой переменной, X - вектор переменных состояния, РХ - вектор производных,

OUT - имя подпрограммы вывода или запоминания результатов, ее параметры: Т - значение независимой переменной, X - вектор переменных состояния, РХ - вектор производных, N - размерность вектора переменных состояния, IER - код завершения, при необходимости досрочно закоичить интегрирование параметру 1ER следует присвоить положительное значение.

Для удобства использования подпрограммы iRKE была написана программа, обеспечивающая ввод исходных данных, интегрирование уравнений состояния и вывод результатов моделирования в виде графиков на АЦПУ или на экран дисплея. В этой программе кроме подпрограммы RKE используется подпрограмма GRAF1K, приведенная в § 12.4.

Для решения конкретной задачи пользователь должен в основной программе задать параметры моделирования, а в подпрограммах FUN н OUT записать выражения для производных переменных состояния и выводимых переменных.

В основной программе необходимо задать следующие параметры: ТО - начальное время моделирования; Т1 - конечное время моделирования; DT - Шаг печати; EPS - допустимая относительная ошибка; N - размерность вектора переменных состояния; Х(1), X(N) - начальные значения вектора переменных состояния; М - количество выводимых переменных; TITLE(l), ... TITLE (М) - 8-байтовые текстовые константы, используемые как заголовки

для выводимых переменных, первый из восьми символов каждой константы используется для печати графика.

В подпрограмме FUN необходимо записать выражения для производных РХ(1), PX(N) в зависимости от переменных состояния Х(1), X(N) и времени Т. Необходимо также присвоить значения выводимым переменным Y(l), .... Y(M). Выводимыми переменными могут быть переменные состояния, их производные, а также функции от переменных состояния и их производных.

Пример. Система описывается нелинейным диффереицнальным уравнением

Произведем моделирование системы при u{t) = l(t). Для этого предварительно составим уравненпя в нормальной форме;

Х\=Х2, X2 = -Xi-0,8X2+U(t].

Результаты моделирования представлены на рис. 9.4.

O.Ot-Ol 4 -4.0t-01 -2

ВРЕМЯ

O.OOt-01 ».........

2.00t4>l » 4.(ЮЬ-01 . » 6.00t-01 . » 8.00Е-01 . I

1.00t*<X) . 1.20t+«> . 1.40t+00 . 1.60fc+<X) . 1.60t+00 .

2.00fct00 ........

2.2tt+00 . 2.4ОЫО0 .

г.&оеюо .

2.80e+00 . 3.00fc*00 . !« 3.20fc+00 .K 3.40t*00 . к 3.40t*00 . я 3.80E100 . 4.00fc+00........

*.2oetoo .

4.40tt00 . 4.60t+00 . 4.8ОЫО0 . 5.00t400 . 5.30E+0O . 5.40E+00 . ъ.мыоо . 5.80t+00 .

A.OOfc+00........

A.20fc+00 . A.40t+00 . 6.А0ЫОО . A.80E+00 . 7.00t+00 . "/.20t+00 . /.40+00 . /.60*00 . /.80fc+00 .

e.oot+oo ........

Ot-01 8.0fc-01 1 Ofc-01 O.Ot-01 2

« .4

2t+00 1

AE+OO CE-01

Ofc+OO OE-01

2.4E->00 » 8.0e-01 » » - X(T)

... O.OOE-01

. 1.90e-OE

. 7.24E-02

. 1.55E-01

. 2.ffiE-01

. 3.B9E-01

. 5.32E-01

«. 6.85E-01

к . 8.40e-01

. 9.87t-01

... I.IIE+OO

. 1.21E+00

. 1.24t+0O

. 1.27E+0O

. 1.23ttoo

. 1.17t+00

. l.lOt+00

. l.Oi+00

. 9.«>£-01

. 9.0BE-O1

... 8.72t-01

. 8.&4E-01

. 8.53t-01

. 8.6(S£-01

. 8.90e-Ol

. 9.21E-01

. 9.56E-01

. 9.90E-01

. 1.02E+00

. 1.04E+00

... l.OAEtOO

. 1.06E+00

. 1.06E+00

. l.CBE+00

. 1.04E+00

. I.026+00

. l.Olt+OO

. 9.91E-01

- 9.79t-01

. 9.72t-01

... 9.i9E-01

» -IiX/PT O.OOE-01 1.85t-0I 3.42E-01 4.76E-01 5.в8Е-01 6.80£-0l 7.t5E-01 7.77E-01 7.63E-01 6.91E-01 5.55E-01 3.64E-01 l.*3E-01 -6.74t-02 -2.3At-0l -3.*0t-0l -3.77E-01 -3.59E-01 -3.01E-01 -2.21E-01 -1.32t-01 -4.60L-02 3.13E-02 9.46e-02 1,40E-01 1.A7E-01 1.74E-01 1.626-01 1.356-01 9.59E-02 5.13E-02 A.94t-03 -3.16E-02 -A.05t-02 -7.74E-02 -8.29b-02 -7.80t-02 -A.50t-02 -4.A9E-02 -2.A5t-02 -6.3bt-03