Современные методы исследования, проектирования и создания сложных объектов и систем неразрывно связаны с разработкой, реализацией на ЭВМ и исследованием математических моделей. Исследование динамических свойств и характеристик таких объектов часто проводится с помощью моделей, представленных в форме систем обыкновенных дифференциальных уравнений (ОДУ), в общем случае, нелинейных и жестких. Для получения решений данные непрерывные [12] математические модели преобразуются в алгоритмические или машинные (дискретные) модели с применением какого-либо численного метода интегрирования ОДУ. Hедостаточно корректное применение того или иного метода дискретизации и выбора шага интегрирования может привести к неадекватности алгоритмической и математической моделей.

Адекватность МОДЕЛЬ – МОДЕЛЬ

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

Далее, в качестве непрерывной математической модели будем рассматриваеть модель в форме системы ОДУ, а в качестве операторов дискретизации (ОД) мы будем подразумевать методы Рунге-Кутты (РК-методы). Модель, получающуюся из непрерывной математической модели (ММ) с помощью ОД, то есть систему разностных уравнений, будем называть алгоритмической моделью (АМ), что вполне соответствует принятой нами терминологии [см. раздел Вопросы и цели].

Динамическая теория

Дискретная модель АМ, в процессе ее реализации на ЭВМ (в виде ПМ), дает возможность “наблюдать” протекающие в ММ процессы в дискретные моменты времени. Но дискретизация, как правило, может приводить к существенным искажениям воспроизведения дискретной моделью динамических свойств и характеристик процессов, описываемых ММ и протекающих в моделируемом объекте.

Качественно и количественно искажения динамических свойств и характеристик АМ, по сравнению со свойствами и характеристиками ММ, полностью определяются динамическими свойствами и характеристиками ОД, применяемыми для дискретизации ММ. Именно на их использовании и основана технология расчета степени адекватности МОДЕЛЬ – МОДЕЛЬ.

Для построения динамических характеристик ОД и пояснения их сути и назначения рассмотрим следующую последовательность трех моделей и их преобразований [36,38], представленную на рис. 3.1.

Рис. 3.1. Последовательность преобразований моделей.

Рассмотрим блоки, представленные на рис. 3.1.

Блок ММ – исходная (для АСМ[29] – возможно, виртуальная) математическая модель в форме системы ОДУ и собственные значения \lambda_{j} ее матрицы Якоби \textbf{L}.

Блок ОД – численный РК-метод дискретизации ОДУ и его символическое векторно-матричное представление. В развернутой форме РК-методы обычно записываются в виде таблицы Батчера [27,57]:

\left.\begin{matrix} 0 \\ c_{2} \\ c_{3} \\ \vdots \\ c_{s} \\ \hline \, \end{matrix}\right| \begin{matrix} \: \\ a_{21} \\ a_{31} & a_{32} \\ \vdots & \vdots \\ a_{s1} & a_{s2} & \ldots & a_{ss-1} & \: \\ \hline b_{1} & b_{2} & \ldots & b_{s-1} & b_{s} \\ \end{matrix}\:\:\:\:\:\rightarrow\:\:\:\:\: \begin{matrix} \textbf{c} & \vline & \textbf{A} \\ \hline \, & \vline & \textbf{b}^T \end{matrix}

(3.1)

где \textbf{A} – матрица коэффициентов, \textbf{c} – вектор абсцисс, \textbf{b} – вектор весов, s – количество этапов явных или стадий неявных РК-методов.

Блок АМ – алгоритмическая модель в форме системы разностных уравнений, получающаяся из исходной ММ при ее дискретизации выбранным РК-методом.

Подчеркнем, что здесь АМ – это разностная модель, получающаяся из исходной ММ путем дискретизации последней численным РК-методом. В данном случае АМ представляет собой обобщение алгоритма численного решения ММ, который в случае простейшей модели:

\dot{x} = f(t,x), \:\:\: x(t_0) = x_0,

(3.2)

при известном значении x_n, состоит, в соответствии с (3.1), из последовательности i = 1, ..., s этапов вычислений вычислений коэффициентов k_i:

\begin{aligned} k_1 & = hf(t,y), \\ k_2 & = hf(t+c_2h, \:\: y+a_{21}k_1), \\ k_3 & = hf(t+c_3h, \:\: y+a_{31}k_1 + a_{32}k_2), \\ \cdots & \:\:\:\:\:\:\:\:\:\:\:\:\:\: \cdots \:\:\:\:\:\:\:\:\:\:\:\: \cdots \\ k_s & = hf(t+c_sh, \:\: y+a_{s1}k_1 + a_{s2}k2 + \cdots + a_{ss-1}k_{s-1}); \end{aligned}

(3.3)

а затем, – вычисления окончательного результата n+1-го шага интегрирования:

x_{n+1} = x_n + \sum_{i=1}^{s} b_ik_i\:.

(3.4)

Собственное значение \rho матрицы \textbf{R} алгоритмической модели – есть функция соответствующего значения \lambda матрицы \textbf{L} исходной ММ, параметров РК-метода (3.1) и шага интегрирования h:

\rho = f(h\lambda,\textbf{c},\textbf{A},\textbf{b}^T).

(3.5)

В теории численных методов дискретизации ОДУ функция, совпадающая с (3.5), называется функцией устойчивости [27], которая применяется для исследования свойств и построения областей абсолютной (вычислительной) устойчивости методов на комплексной плоскости параметров h\lambda.

Четыре известных и наиболее употребительных РК-метода с s = p = 1, ..., 4, независимо от величин их коэффициентов в таблице (3.1), имеют следующие функции устойчивости:

\rho_1 = 1 + h\lambda,

(3.6)

\rho_2 = 1 + h\lambda + \frac{1}{2!} h^2 \lambda^2,

(3.7)

\rho_3 = 1 + h\lambda + \frac{1}{2!} h^2 \lambda^2 + \frac{1}{3!} h^3 \lambda^3,

(3.8)

\rho_4 = 1 + h\lambda + \frac{1}{2!} h^2 \lambda^2 + \frac{1}{3!} h^3 \lambda^3 + \frac{1}{4!} h^4 \lambda^4,

(3.9)

которые являются полиномами h\lambda и совпадают с соответствующими отрезками ряда Тейлора для e^{h\lambda}.

Параметр \lambda, здесь, – это любой элемент \lambda_{j} матрицы \textbf{L} исходной ММ. Иначе говоря, функция устойчивости РК-метода всегда определена на спектре матрицы Якоби системы ОДУ, для численного решения которой применяется этот метод. Но тогда подстановкой в функцию устойчивости метода самой матрицы \textbf{L} исходной ММ, вместо \lambda, корректно определяется [45,58] функция этой матрицы:

\textbf{R} = f(h\textbf{L},\textbf{c},\textbf{A},\textbf{b}^T),

(3.10)

которая, в данном случае, непосредственно представляет собой матрицу \textbf{R} алгоритмической модели. Далее, если известны все собственные значения \lambda_{j} матрицы \textbf{L} исходной ММ, то соответствующие им собственные значения \rho_{j} матрицы \textbf{R} алгоритмической модели, также непосредственно вычисляются путем подстановки значений \lambda_{j} в формулу (3.5), или одну из формул (3.6) – (3.9), в зависимости от применяемого к исходной ММ численного метода.

Блок МД – модель-двойник исходной MM. Модель-двойник (MД), так же как и исходная MM, представляет собой систему ОДУ, кроме того, размерности MД и MМ совпадают. В действительности MД, в отличие от АM, не существует, однако, мы будем работать с этой моделью как с реально существующей. Нас будут интересовать собственные значения \eta_{j} ее матрицы \textbf{G}, а акже ее точное решение, которое мы будем сравнивать с точным, пусть даже нам не известным, решением исходной MМ.

Сопоставление моделей АМ и МД подсказано методом дискретной фазовой плоскости [23] исследования импульсных систем управления и оказалось весьма полезно для исследования динамических свойств РК-методов [35].

Воссоздание МД будем осуществлять, предполагая выполнение следующих условий:

Размерность МД в точности соответствует размерности АМ.
Параметры МД рассчитываются, исходя из условия совпадения решений МД и АМ в дискретные моменты времени t_{n} = nh, то есть мы подразумеваем точное выполнение граничных условий (ГУ), представленных в Блоке ГУ на рис. 3.1.

Рассмотрим решение этой задачи для случая однородных моделей:

\begin{aligned} \textbf{x}_n & = \textbf{R}\textbf{x}_{n-1}, \\ \dot{\textbf{y}} & = \textbf{G}\textbf{y}, \end{aligned}

(3.11)

с заданными граничными условиями:

\left.\begin{aligned} \textbf{y}(t_0) & = \textbf{y}_0 = \textbf{x}_0, \\ \textbf{y}(t_m) & = \textbf{x}_m\:. \end{aligned}\right\}

(3.12)

Здесь t_{m} = t_0 + mh, при t_0 = 0, t_m = mh, h – шаг интегрирования.

Решение \textbf{x}_{m} алгоритмической модели в точке t_m = mh, при известном решении \textbf{x}_{0} в точке t_0 = 0, выражается формулой:

\textbf{x}_{m} = \textbf{R}^m \textbf{x}_{0}\:.

(3.13)

В свою очередь, решение \textbf{y}(t_m) модели-двойника в точке t_m = mh, при известном решении \textbf{y}(t_0) = \textbf{y}_0, имеет вид:

\textbf{y}(t_m) = e^{\textbf{G}mh} \textbf{y}_{0}\:.

(3.14)

Сравнивая решения (3.13) и (3.14), с учетом (3.12), получим:

\textbf{R}^m = e^{\textbf{G}mh}

(3.15)

и, после логарифмирования, а также сокращения на m

\textbf{G} = \frac{1}{h} ln \textbf{R}\: .

(3.16)

На основании известных [57,58] свойств функций матриц, для собственных значений \gamma_{j} матрицы \textbf{G} модели-двойника (3.16) получим:

\gamma_{j} = \frac{1}{h} ln \rho_{j}

(3.17)

или

ln \rho_{j} = h \gamma_{j}\: ,

(3.18)

где \rho_{j} – любое из собственных значений матрицы \textbf{R} алгоритмической модели.

Однако \rho_{j}, в свою очередь, является функцией собственного значения \lambda_{j} матрицы \textbf{L} исходной ММ и эта функция, для любого \lambda_{j}, есть ни что иное, как функция устойчивости (3.5) того численного метода, который был применен для дискретизации ММ. Таким образом,

ln\:\rho = ln\:\rho(h\lambda,\textbf{c},\textbf{A},\textbf{b}^T).

(3.19)

Так как \lambda – любое из собственных значений \lambda_{j} матрицы \textbf{L} исходной ММ, то функция (3.19), через(3.17), определяет преобразование собственных значений \lambda_{j} (спектра) матрицы \textbf{L} исходной ММ в соответствующие собственные значения \gamma_{j} (спектр) матрицы \textbf{G} виртуальной модели-двойника.

Решения МД и АМ совпадают в дискретные моменты времени. Поэтому АМ является точным дискретным имитатором непрерывной МД, которая, в свою очередь, является непрерывным двойником исходной ММ. При анализе свойств численных методов и исследовании, на их основе, адекватности МОДЕЛЬ – МОДЕЛЬ мы будем сопоставлять ММ и АМ. Но при этом всегда будем иметь в виду, что через последнюю отображается структура, параметры и динамические свойства МД, как непрерывного (виртуального) двойника исходной непрерывной модели.

Свойства и характер решений МД и ММ полностью определяются собственными значениями (с учетом физического смысла их компонентов) и элементарными делителями их матриц \textbf{G} и \textbf{L}. Поскольку компоненты \lambda_{j} и \gamma_{j}, в основном, определяют динамические свойства моделей, функцию (3.19) преобразования \lambda_{j} в \gamma_{j} можно рассматривать как оператор преобразования динамических свойств исходной ММ в динамические свойства АМ. В то же время, эти преобразования определяются через функцию устойчивости применяемого для дискретизации ММ численного метода. Следовательно, функция (3.19) – есть некая характеристика свойств численных методов дискретизации ОДУ. Итак.

Определение 1. Динамической характеристикой \textbf{D} численного метода дискретизации ОДУ называется логарифм натуральный его функции устойчивости

\textbf{D} = ln\:\rho .

(3.20)

Так как, в общем случае, \lambda_{j}, \rho_{j} и \gamma_{j} – комплексные величины, то функцию \textbf{D} можно представить в виде комплексной динамической характеристики

\textbf{D} = ln\:\rho = ln(\sigma \pm i\tau) = h\gamma = h(\eta \pm i\omega),

(3.21)

а также в виде вещественной

\textbf{D}_{R} = Re\:\textbf{D} = \frac{1}{2}ln(\sigma^2 + \tau^2) = h\eta,

(3.22)

и мнимой

\textbf{D}_{I} = Im\:\textbf{D} = arctg \frac{\tau}{\sigma} + 2k\pi = h\omega, \:\:\:\: k = 0,\:\pm1,\:\pm2,...

(3.23)

динамических характеристик численного метода.

Вещественную \textbf{D}_{R} = h\eta и мнимую \textbf{D}_{I} = h\omega динамические характеристики можно изобразить в виде соответствующих поверхностей в трехмерном пространстве над плоскостью параметров h\lambda\:\:(Re\:h\lambda,Im\:h\lambda).

Укажем основные особенности динамических характеристик численных методов, которые могут быть использованы при реализации алгоритмов оценки и обеспечения адекватности МОДЕЛЬ – МОДЕЛЬ.

Графики \textbf{D}_{R} = h\eta, как функций определенных на плоскости параметров h\lambda пересекаются с этой плоскостью строго по границам \partial \Re областей устойчивости соответствующих численных методов, а внутри областей устойчивости, функции \textbf{D}_{R} = h\eta имеют отрицательные значения.
Поверхности \textbf{D}_{R} = h\eta, и \textbf{D}_{I} = h\omega динамических характеристик над некоторой окрестностью начала координат плоскости h\lambda, или практически полностью совпадают с наклонными плоскостями, изображающими соответственно вещественную h\mu и мнимую h\nu компоненты, умноженного на h собственного значения \lambda ММ, или незначительно отклоняются от них. Поэтому на той же плоскости параметров h\lambda можно выделить области, где относительные отклонения их друг от друга не превышают определенных, заранее заданных, величин:

\delta_R = \frac{h\eta - h\mu}{h\mu} = \frac{\eta - \mu}{\mu},

(3.24)

\delta_{I} = \frac{h\omega - h\nu}{h\nu} = \frac{\omega - \nu}{\nu}.

(3.25)

Пересечение этих двух областей на комплексной плоскости h\lambda образует непустую область \Re_{\delta} с границей \partial \delta, где одновременно ни одна из величин \delta_{R} и \delta_{I} не превышают заданной величины \delta. Из этого вытекает

Определение 2. Областью динамической \delta-согласованности численного метода дискретизации ОДУ, на комплексной плоскости h\lambda, называется множество всех h\lambda \in \Re_{\delta}, для которых ни одна из величин \delta_{R} (3.24) и \delta_{I} (3.25) не превышает заданного значения \delta.

Введем теперь понятие рабочих точек алгоритмической модели.

Определение 3. Рабочими точками АМ называются точки на комплексной плоскости h\lambda, соответствующие произведениям h\lambda_{j} всех собственных значений \lambda_{j} ММ на шаг интегрирования h.

Отметим, что рабочие точки h\lambda_{j} АМ и ее собственные значения \rho_{j} являются по сути разными объектами. В то же время при h=1 рабочие точки АМ являются собственными значениями \lambda_{j} ММ на комплексной плоскости h\lambda.

Теперь становится ясным основной смысл и практическое назначение областей динамической \delta-согласованности численных методов дискретизации ОДУ, которые заключаются в следующем.

Если все рабочие точки АМ, при выбранном шаге интегрирования h, находятся внутри заданной области \delta-согласованности применяемого РК-метода, то АМ воспроизводит процессы происходящие в ММ и их динамические свойства с погрешностью, не превышающей заданной величины \delta. В этом смысле можно говорить о \delta-адекватности АМ и ММ с точностью \delta.

Как следует из изложенного, каждый конкретный метод дискретизации ОДУ имеет собственную динамическую характеристику, которая, по сути, представляет собой некий оператор преобразования спектра матрицы \textbf{L} исходной непрерывной ММ в спектр матрицы \textbf{G}, также непрерывной, МД. На точном решении последней лежит дискретная последовательность точек точного решения АМ. Но точное решение АМ является, в свою очередь, тем приближенным численным решением исходной ММ, которое мы наблюдаем. В традиционных (классических) исследованиях свойств и погрешностей численных методов обычно сравниваются именно решения непрерывной ММ и дискретной АМ. Использование же динамических характеристик методов дает возможность сравнивать между собой спектры, а через них и решения двух непрерывных моделей ММ и МД, и на этой основе строить алгоритмы расчета точных (не оценочных, в отличие от классики) характеристик погрешностей моделирования, а также алгоритмы целенаправленной коррекции динамических свойств и характеристик дискретных моделей.

Алгоритмы расчета

Дифференциальные модели реальных объектов, в общем случае, нелинейны. Поэтому алгоритмы оценки адекватности МОДЕЛЬ – МОДЕЛЬ связаны, прежде всего, с необходимостью формирования матрицы Якоби \textbf{L} непрерывной ММ. Поскольку в распоряжении компилятора языка моделирования, после синтаксического и семантического анализа, имеются функции расчета правых частей системы ОДУ, представляющих исходную ММ, матрицу Якоби можно формировать численно, например, следующим образом.

Алгоритм численного формирования матрицы Якоби нелинейной системы ОДУ

Для работы нам понадобятся формулы численного дифференцирования. Приведем здесь несколько формул из таблицы 7.12 (см. [19], с.510):

Формула `y_{j} = y(jh),\:\:\: {{y}_{j}}' = {y}'(jh), ...`	Оценка погрешности	Номер формулы
`{y_0}' = \frac{1}{2h} (-y_{-1} + y_1)`	`- \frac{1}{6} h^2{y_0}''' - ...`	(F1)
`{y_0}' = \frac{1}{12h} (y_{-2} - 8y_{-1} + 8y_{1} - y_2)`	`+ \frac{1}{30} h^4{y_0}^{(5)} + ...`	(F2)
`{y_0}' = \frac{1}{h} (-y_0 + y_1)`	`- \frac{1}{2} h{y_0}'' - ...`	(F3)
`{y_0}' = \frac{1}{2h} (-3y_0 + 4y_1 - y_2)`	`+ \frac{1}{3} h^2{y_0}''' + ...`	(F4)
`{y_0}' = \frac{1}{12h} (-3y_{-1} -10y_0 + 18y_1 - 6y_2 +y_3)`	`- \frac{1}{20} h^4{y_0}^{(5)} + ...`	(F5)

(3.26)

Запишем систему нелинейных дифференциальных уравнений в виде

\left.\begin{aligned} \frac{dx_1}{dt} & = f_1(t,x_1,x_2,...,x_n,u_1),\\ \frac{dx_2}{dt} & = f_2(t,x_1,x_2,...,x_n,u_2),\\ \:...\:\: & = \:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\: ...\\ \frac{dx_n}{dt} & = f_n(t,x_1,x_2,...,x_n,u_n), \end{aligned}\right\}

(3.27)

и будем формировать матрицу Якоби

\textbf{W}(\textbf{x})=\begin{pmatrix} \frac{\partial f_1}{\partial x_1} & \frac{\partial f_1}{\partial x_2} & \cdots & \frac{\partial f_1}{\partial x_n} \\ \frac{\partial f_2}{\partial x_1} & \frac{\partial f_2}{\partial x_2} & \cdots & \frac{\partial f_2}{\partial x_n} \\ \vdots & \vdots & & \vdots \\ \frac{\partial f_n}{\partial x_1} & \frac{\partial f_n}{\partial x_2} & \cdots & \frac{\partial f_n}{\partial x_n} \end{pmatrix} = \begin{pmatrix} \frac{\partial f_i}{\partial x_j} \end{pmatrix}

(3.28)

приближенно в виде

\textbf{J}_c = \begin{pmatrix} {y_{11}}' & {y_{12}}' & \cdots & {y_{1n}}' \\ {y_{21}}' & {y_{22}}' & \cdots & {y_{2n}}' \\ \vdots & \vdots & & \vdots \\ {y_{n1}}' & {y_{n2}}' & \cdots & {y_{nn}}' \end{pmatrix} \approx \textbf{W}(\textbf{x}),

(3.29)

где {y_{ii}}', – подлежащие вычислению приближенные значения соответствующих частных производных.

Зададимся приращениями компонентов вектора состояния модели (3.27)

h_i = \Delta x_i, \:\:\:\: i=1, ...,n.

(3.30)

На каждом шаге численного интегрирования системы (3.27) выбор \Delta x_i можно осуществить следующим образом:

\Delta x_i = \delta(x_i[n]-x_i[n-1]), \:\:\:\: \delta = 0.1;\:0.01;\:...

(3.31)

где, \delta выбирается исходя из требуемой точности вычислений, а x_i[n], x_i[n-1] – значения компонента решения x_i на текущем и предыдущем шагах интегрирования системы, соответственно.

Для вычисления h_i в начале счета, необходимо сделать один шаг интегрирования системы и вычислить

\Delta x_i = \delta(x_i[1]-x_i[0]).

(3.32)

Это, в какой-то мере, избавит от проблем вычисления h_i при нулевых начальных условиях, т.е. при x_i[0] = 0. Правда в этом случае может оказаться, что (x_i[1] - x_i[0]) = 0, но это легко обойти, сделав, например, еще один шаг интегрирования или объявить \Delta x_i = \delta.

Будем формировать значение матрицы \textbf{J}_{c} по столбцам, используя (F2) из набора (3.26). Для первого столбца зададимся конкретной величиной h_1 и вычислим необходимые значения функции f_1() системы (3.27), при соответствующих приращениях компонента x_i:

\begin{aligned} y_{-2} &= f_1(t,x_1-2h_1,x_2,...,x_n,u_1),\\ y_{-1} &= f_1(t,x_1-\:\:h_1,x_2,...,x_n,u_1),\\ y_{1} &= f_1(t,x_1+\:\:h_1,x_2,...,x_n,u_1),\\ y_{2} &= f_1(t,x_1+2h_1,x_2,...,x_n,u_1). \end{aligned}

(3.33)

Естественно, что y_0=f_1(t,x_1,x_2,...,x_n,u_1).

Тогда величину {y_{11}}' можно вычислить по формуле (F2) следующим образом:

{y_{11}}' = \frac{1}{12h} (y_{-2} - 8y_{-1} + 8y_{1} - y_2).

(3.34)

Для вычисления {y_{21}}' необходимо сначала вычислить значения функции f_2() при тех же приращениях компонента x_1 и, аналогично (3.34), получить значение {y_{21}}'.

Данные действия надо повторить для всех функций f_3(), ..., f_n(), вычисляя затем значения {y_{31}}', ..., {y_{n1}}' матрицы \textbf{J}_{c}.

Повторяя действия, описанные выше, для h_2,h_3,...,h_n, давая приращения компонентам x_2,x_3,...,x_n вектора решения системы (3.27), мы получим значения остальных столбцов матрицы \textbf{J}_{c}.

Перспективное оценивание адекватности МОДЕЛЬ – МОДЕЛЬ

В обобщенной форме, алгоритм оценивания \delta-адекватности МОДЕЛЬ – МОДЕЛЬ представляется последовательностью следующих этапов:

1 этап. Формирование матрицы Якоби сходной ММ. Если ММ приведена к нормальному виду и является моделью с постоянными коэффициентами, то формирование матрицы Якоби сводится к формированию матрицы коэффициентов \textbf{L}, соответствующей однородной модели.

2 этап. Расчет (оценка) наибольшего (по модулю) собственного значения \lambda_{max}, сформированной матрицы Якоби \textbf{L}. Такой расчет может быть выполнен по какой-либо программе из библиотек прикладного математического обеспечения языков программирования.

3 этап. Расчет по динамической характеристике \textbf{D} = ln\:\rho, (3.20), применяемого численного метода дискретизации ОДУ, собственному значению \lambda_{max}, полученному на 2-м этапе, и выбранному шагу интегрирования h, собственного значения \gamma матрицы \textbf{G} модели-двойника для исходной модели ММ.

4 этап. Расчет, по компонентам, полученных на 2-м и 3-м этапах собственных значений \lambda_{max} = \mu \pm i\nu, \gamma = \eta \pm i\omega и, по формулам (3.24) и (3.25), их относительных рассогласований \delta_{R} и \delta_{I}. Максимальное из этих значений принимается за величину \delta оценки степени \delta-адекватности МОДЕЛЬ – МОДЕЛЬ.

На этом работа алгоритма оценки \delta-адекватности МОДЕЛЬ – МОДЕЛЬ. заканчивается.

Примеры таких расчетов можно найти в [38].

В данном алгоритме рассчитывается только одно, максимальное по величине, собственное значение \lambda_{max} матрицы \textbf{L} исходной ММ. Этого вполне достаточно, т.к. при этом относительные отличия остальных \lambda_{j}, меньших \lambda_{max}, от соответствующих \gamma_{j}, будут заведомо меньше \delta, полученного на 4-м этапе. Иначе говоря, если рабочая точка h\lambda_{max} АМ лежит в области заданной \delta-согласованности, выбранного численного метода, то остальные рабочие точки h\lambda_{j} будут тем более находиться внутри этой области.

Ретроспективное оценивание адекватности МОДЕЛЬ – МОДЕЛЬ

Второй способ оценивания адекватности МОДЕЛЬ – МОДЕЛЬ основан на анализе результатов вычислительных экспериментов и решении обратных задач теории динамических свойств численных методов дискретизации ОДУ [36,38]. Его целесообразно применять для оценивания адекватности моделей вновь (или впервые) разрабатываемых сложных объектов, когда непрерывная модель оказывается слишком сложной, нелинейной и не позволяет даже ориентировочно оценить свойства моделируемого объекта, не прибегая к непосредственной предварительной реализации на ЭВМ.

Данный способ предполагает использование алгоритмических аналогов анализаторов частоты и амплитуды некоторых наиболее динамичных или подозрительных процессов на выходе ПМ (анализ следа модели). В результате работы таких анализаторов фактически получаются оценки числовых величин собственных значений \gamma = \eta \pm i\omega непрерывной модели, точному решению которой принадлежит анализируемая дискретная последовательность значений на выходе ПМ. По полученным, таким образом, приближенным значениям \gamma непрерывной модели-двойника, и динамическим характеристикам ОД вычисляются собственные значения \lambda, непрерывной модели, по которой построена данная ПМ.

Поясним суть данного метода на простом примере.

Пусть в результате работы анализаторов приближенно вычислены собственные значения

y_{1,2} = 0.0972 \pm i 0.5689

(3.35)

непрерывной модели, определяющие некоторое ее точное решение, которому принадлежит проанализированная дискретная последовательность значений, соответствующего процесса в ПМ. Пусть также в АМ (ПМ) используются ОД, соответствующий методу Эйлера (РК-1) дискретизации ОДУ и шаг равный единице h=1.

Обратим внимание на то, что ПМ, в данном случае, неустойчива.

По динамическим характеристикам данного ОД (3.6) находим

1 + \lambda h = e^{\gamma h}

откуда, при h=1, получим

\lambda = e^{\gamma} - 1.

Подстановка сюда числового значения \gamma (3.35) дает

\lambda_{1,2} = e^{0.0972}(cos\:0.5689 \pm i\:sin\:0.5689) - 1 = -0.0715 \pm i\:0.5937.

(3.36)

Сравнивая между собой полученные результаты (3.35) и (3.36), можно строго количественно оценить степень динамической адекватности ПМ непрерывной модели. В данном случае, если сравнивать собственные значения по модулю, то получим \delta < 4\%. Однако здесь ПМ оказалась неустойчивой, в отличие от исходной непрерывной модели.