Непосредственное создание языка моделирования следует начинать с разработки синтаксических анализаторов различных видов математических моделей. Далее, имея некоторые предположения о структуре объектного кода, необходимо разработать, соответствующие предложениям языка, семантические процедуры, которые обеспечат последующee генерирование объектного кода.

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

Для рассмотрения вопросов, связанных с семантикой математических моделей, воспользуемся стилем задания грамматик, принятым для входных спецификаций таких программ как Yacc, Bison, ZUBR [40] и близким к нормальной форме Бэкуса (BNF, НФБ) [21,49].

Например, грамматику целых чисел можно записать следующим образом:

%%
число:   чс
       ;
чс:      чс  цифра
       | цифра
       ;
цифра:   '0'
       | '1'
       | '2'
       | '3'
       | '4'
       | '5'
       | '6'
       | '7'
       | '8'
       | '9'
       ;
%%

Здесь символы число, чс и цифра являются нетерминальными; цифры заключены в одинарные кавычки, что говорит генератору синтаксических анализаторов о том, что они будут поступать в виде символьных переменных [33] непосредственно от лексического анализатора, то есть данные символы являются терминальными.

Кроме того, во входной спецификации можно определять код на языке C, который должен быть ограничен символами %{ и %}. Например, так:

%{

int number = 7;

%}

Данную возможность мы будем использовать для упрощения лексической привязки к контексту разбираемого входного потока и связанных с этим пояснений.

Системы обыкновенных дифференциальных уравнений

Рассмотрим простую гамматику:

%{
   int maxp = 0; /* счетчик порядка производных  */
   int vtmp = 0; /* счетчик временных переменных */
   int ctmp = 0; /* счетчик числовых констант    */
%}

%union
{
   SYMBOL *sym;
}
%token    CONST NUMBER
%token    VARIABLE BUILTIN DIFF
%type     expr
%right '='
%left '+' '-'   /* left associative, same precedence */
%left '*' '/'   /* left associative, higher precedence */
%left UNARYMINUS
%right '^'      /* power */
%%

deset: /* пусто */
       | deset
       | deset de
       ;

de:
         DIFF '=' expr ';'
       ;

expr:    CONST
       | VARIABLE
       | DIFF
       | expr '+' expr
       | expr '-' expr
       | expr '*' expr
       | expr '/' expr
       | expr '^' expr   /* возведение в степень */
       | '(' expr ')'
       | BUILTIN '(' expr ')'
       | '-' expr %prec UNARYMINUS
       ;

%%

Для синтаксического анализатора, построенного по данной грамматике, система ОДУ может быть задана в виде списка deset, состоящего из отдельных уравнений de. Дифференциальное уравнение, в свою очередь, должно представлять собой запись, состоящую из двух частей. Слева от знака равенства должна располагаться производная от переменной состояния (терминал DIFF). Справа – выражение, в состав которого могут входить: числовые константы (CONST), переменные (VARIABLE), производные (DIFF) и вызовы встроенных функций (BUILTIN), таких, например, как exp(), sin(), cos() и т.д. Для задания выражений (expr) в правой части допускаются операторы сложения (+), вычитания (-), умножения (*), деления (/), возведения в степень (^), неограниченной вложенности скобочные выражения и операции унарного минуса.

Допустим, что лексический анализатор при получении из входного потока, например, следующих символов:

x''

расценивает их как вторую производную от переменной x и возвращает синтаксическому анализатору указатель (*sym) на созданную им (для производной от x) запись в таблице символов. Данная запись кроме самого идентификатора переменной x содержит поле p, значением которого является порядок производной:

struct sym
{
   char *id;
   int    p;

   . . .
};

Для предотвращения появления в правой части уравнений производных равного или более высокого порядка, чем у производных стоящих в левой части, мы ввели глобальную переменную maxp. Начальное значение этой переменной равно нулю maxp = 0. Мы будем увеличивать это значение на единицу при каждом выводе по правилу

expr: DIFF

следующим образом:

expr:
       | DIFF
         {
            maxp = max( $1->p, maxp );
            . . .
         }
       ;

Здесь мы ввели семантическую процедуру, написанную на языке C и ограниченную, с обеих сторон, фигурными скобками. Псевдопеременная $1 служит для ссылки на значение первой компоненты правила DIFF, которое было возвращено лексическим анализатором.

После успешной проверки во время вывода по правилу:

de:
       | DIFF '=' expr ';'
         {
            if( maxp >= $1->p )
               error( "уравнение не разрешено относительно"
                      "старшей производной" );
            else
            {
               /*
                 Обнуляем maxp перед разбором очередного
                 уравнения
                */
               maxp = 0;
            }

            . . .

         }
       ;

мы можем обнулить переменную maxp и приготовиться к анализу следующего уравнения.

Справочный материал

Прежде чем рассматривать семантику ОДУ, необходимо определить структуру записи в таблице символов и способ внутреннего представления программ. Мы не будем подробно рассматривать все необходимые структуры и функции работы с ними. Более того, код на языке C, приводимый здесь, следует рассматривать как примерное описание на метаязыке, служащее для пояснения действий, связанных с семантическим анализом, а не как готовые программы. Следующий материал нужен лишь для того, чтобы создать основу, на базе которой мы сможем продолжить дальнейшее обсуждение.

Таблица символов

В таблице символов мы будем хранить записи, соответствующие символам программы. Такими символами могут быть: числовые константы, обычные переменные, производные и встроенные функции. О них мы говорили ранее. Следовательно, в структуре записи таблицы символов нам понадобится целочисленное поле (int type), определяющее тип переменной. Для хранения идентификатора символа будем использовать поле char *id, являющееся указателем на строковую переменную. Фактическое значение будем хранить в поле double val, представляющем собой вещественное число (для простоты ограничимся одним типом данных). Для производных введем целочисленное поле int p, которое будет содержать порядок производной. Способ описания и хранения встроенных функций, здесь, не имеет приципиального значения, поэтому рассматривать его мы не будем.

Кроме того, нам понадобятся специальные флаги, которые мы введем позднее.

Итак, структуру записи таблицы символов на языке C можно определить следующим образом:

typedef struct sym sym;
struct sym
{

   int     type;  /* DIFF, VARIABLE, CONST, BUILTIN */

   char   *id;    /* идентификатор */
   double  val;   /* значение */
   int     p;     /* порядок производной */

   int     flags; /* флаги */

   . . .

};

Представим, что с помощью функции install(char *id, int type, double val), принимающей в качестве аргументов идентификатор, тип и значение символа, мы сможем заносить в таблицу необходимые записи. Возвращаемым значением функции install() будет указатель *sym на запись в таблице, которая создана с ее помощью. Для поиска в таблице символов записи по идентификатору будет служить функция lockup(char *id), возвращающая указатель *sym на найденную запись.

Внутреннее представление

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

Поля up_tetrade и down_tetrade служат для организации двусвязного списка тетрад, который и будет представлять содержание программы.

typedef struct tetrade tetrade;
struct tetrade
{
   sym   *func;
   sym   *result;
   sym   *operand1;
   sym   *operand2;

   struct tetrade *up_tetrade;
   struct tetrade *down_tetrade;
};

Для более удобной работы с программами, как с едиными конструкциями, определим структуру program следующего вида:

typedef struct program program;
struct program
{
   struct tetrade *root_up;
   struct tetrade *root_down;

   struct tetrade *ptr_tetrade;
};

Теперь мы можем работать с программой посредством указателя. Это мы сделаем, определив переменную text_prog:

static program text_prog;

и глобальный указатель на нее

program *current_prog;

Размещение программы в памяти можно представить схемой, показанной на рис. 4.1.

Для создания и добавления тетрад (структур tetrade) в программу нам понадобится функция add_tetrade():

extern tetrade *add_tetrade( program *p,
                             sym     *fn,
                             sym     *rc,
                             sym     *op1,
                             sym     *op2 );

первым параметром которой является указатель на программу. Остальные параметры функции add_tetrade() представляют собой указатели на записи в таблице символов.

Вот все, что нам понадобится для перехода к началу семантического анализа ОДУ.

Семантический анализ

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

expr: expr '+' expr

можно поставить следующую семантическую процедуру

expr:
         . . .

       | expr '+' expr
         {
            char t_name[80], t_num[80];
            sym *tmp;
            sym *func;

            strcpy( t_name, "_$tmp$_" );
            sprintf( t_num, "%u", vtmp++ );
            strcat( t_name, t_num );

            tmp = install( t_name, VARIABLE, 0.0 );
            func = lookup( "add" );
            add_tetrade( current_proc, func, tmp, $1, $3 );
            $$ = tmp;
         }

         . . .

       ;

Здесь создается временная переменная с идентификатором _$tmp$_X, где X – текстовое представление значения глобальной переменной itmp, которая служит счетчиком создаваемых таким образом переменных. Тетрада формируется из указателей на записи в таблице символов для: функции add(), временной переменной, созданной для хранения результата сложения, и операндов. Указатели на записи таблицы символов, созданные для операндов, хранятся во внутреннем стеке синтаксического анализатора. Обращение к ним осуществляется посредством псевдопеременных $1 и $3. Данные указатели были созданы ранее при выводе нетерминальных символов expr, присутствующих в правой части правила. Более подробно об этом можно прочитать в [40]. Заметим здесь, что функция add(), как и другие подобные функции, реализующие коды операций виртуальной машины, может быть реализована на языке C и ее описатель должен быть добавлен в таблицу символов заранее [40].

Поскольку результат данной операции нам понадобится позже, при выводе других правил, для построения очередных тетрад программы, мы присвоили значение указателя *tmp на результирующую переменную символу expr, находящемуся в левой части правила.

Оператор задержки

Кроме обычных арифметических операций нам понадобится оператор задержки \Delta t. Договоримся о том, что операторы \Delta t будут срабатывать в программе только после того, как сработали операторы других видов. Все операторы \Delta t работают по, так называемым, переменным состояния модели и срабатывают одновременно, пропуская на выход значение входа [29], полученное путем выполнения всех операторов программы. Эти операторы позволяют описывать рекурентные вычисления. Последние, в свою очередь, дают возможность в дискретной форме описывать динамику моделируемого объекта. Значения выходов операторов \Delta t для первого шага моделирования всегда задаются (исходное состояние системы или начальные условия).

Семантика производных

При разборе правил, в которых присутствует терминальный символ DIFF, дело обстоит иначе.

Введем новую операцию нашей виртуальной машины – delay(), и договоримся представлять интеграторы в виде двух команд виртуальной машины так, как мы делали это в разделе Вопросы и цели при переходе от (1.4) к (1.5) (см. рис. 4.2):

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

Итак, на рис. 4.2 показан принцип формирования имен создаваемых переменных при семантическом анализе производных переменной x. Для формирования имени первой производной к имени переменной "x", слева, добавляется строка "_1p_", для второй – строка "_2p_" и так далее. Имя входной переменной для оператора \Delta t создается путем добавления к имени выходной переменной (оператора \Delta t ) строки "_plus1".

Рассмотрим процедуру семантического анализа второй производной от x. При встрече, во входном потоке, цепочки x'', как мы договорились ранее, лексический анализатор возвращает имя переменной "x" и порядок ее производной p = 2. В рамках семантической процедуры мы должны выполнить следующие действия:

• Создаем в таблице символов новую переменную с именем "_2p_x".
• Создаем переменные с именами "_1p_x_plus1", "_1p_x", "x_plus1", "x", и отмечаем их как входные и выходные переменные оператора \Delta t (см. элемент flags структуры записи (sym) в таблице символов).
• Формируем следующую последовательность команд во внутреннем представлении программы:
  
  

  	КОП	Результат	Операнд 1	Операнд 2
  	+	_1p_x_plus1	_1p_x	_2p_x
  	\Delta t	_1p_x	_1p_x_plus1

  (4.1)

  
  

  	+	x_plus1	x	_1p_x
  	\Delta t	x	x_plus1

  (4.2)

  
  Следует отметить важное обстоятельство. Если при анализе второй производной оказалось, что первая производная уже была разобрана ранее (это легко определить по наличию в таблице символов переменной с именем "_1p_x"), то семантическую процедуру для первой производной повторять нет необходимости.
• И последнее действие, которое надо совершить, это вернуть значение данной семантической процедуры синтаксическому анализатору. В случае второй производной этим значением будет указатель на запись в таблице символов для переменной с именем "_2p_x" (естественно, если мы анализировали только первую производную x', то возвращаемым значением будет указатель на "_1p_x"). На языке C для генераторов синтаксических анализатороа, таких как YACC, BISON, ZUBR, это можно записать с помощью следующих операторов:
  

  sym  *tmp;
  
  tmp = install( "_2p_x", VARIABLE, 0.0 );
  $$  = tmp;

Процедура семантического анализа ОДУ

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

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

x'' = k^2*(a*x' - x); 

Прежде чем анализатор дойдет до синтаксической единицы, которую он сможет разобрать согласно какого-либо правила грамматики ОДУ, ему придется прочитать следующую цепочку символов:

x'' = k^2*(a*x' 

При этом в таблице символов лексический анализатор создаст следующие переменные и числовые константы:

Разумеется, все значения прочитанных лексических единиц будут помещены в стек синтаксического анализатора и, будут доступны когда настанет их очередь в процессе синтаксического анализа [40].

Первой синтаксической единицей, к которой можно применить правило вывода, является последовательность "a*x'", однако, перед выполнением синтаксической процедуры, соответствующей операции умножения, необходимо разобрать производную "x'" по методике (4.2). В результате анализа первой производной в таблицу символов будут добавлены переменные с именами "_1p_x", "x_plus1", "x" и, кроме того, переменные "x_plus1", "x" будут отмечены как входная и выходная переменные оператора задержки \Delta t, соответственно. Внутреннее представление программы, после разбора первой производной, будет содержать следующие два оператора:

Результирующим значением семантической процедуры разбора первой производной будет указатель на переменную с именем "_1p_x". Данный указатель понадобится при формировании следующих тетрад программы.

Далее по правилу, описывающему операции умножения, синтаксический анализатор сможет разобрать цепочку символов "a*x'", что даст еще одну тетраду нашей программы:

Здесь "_$tmp$_1", – временная переменная, получившая уникальное имя согласно счетчику vtmp и являющаяся результатом семантической процедуры разбора последовательности "a*x'".

Поскольку на данный момент синтаксический анализатор не может более произвести никаких выводов, лексический анализатор прочитает еще пару символов последовательности (4.3):

x'' = k^2*(a*x' - x

Теперь наши анализаторы в состоянии сгенерировать тетраду для цепочки "a*x’ – x":

Далее, чтение цепочки (4.3) приводит нас к концу уравнения ";", где анализаторы оказываются в состоянии вывести значение правой части уравнения:

Наконец, по правилу вывода дифференциального уравнения:

de:   DIFF '=' expr ';'
         ;

и применяя процедуру синтаксического анализа второй производной, получим программу:

где, последний оператор, фактически, присваивает значение правой части уравнения (4.3) второй производной "x''".

Поскольку переменная _2p_x является производной наивысшего порядка (_2p_x не является выходом оператора задержки \Delta t) и ее значение приравнивается значению правой части уравнения, мы можем удалить из программы оператор mov и заменить переменную _2p_x на переменную _$tmp$_4. В результате получим заключительную версию программы:

где в скобках, справа от каждой тетрады, показан номер уравнения системы ОДУ в нормальной форме Коши.

Если теперь мы произведем следующую замену имен переменных

то увидим, что программа (4.10) дает схему вычислений

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

Более того, переменные и операторы, помеченные в программе (4.10) серым цветом, являются ни чем иным как процедурами вычисления значений правых частей уравнений

системы в нормальной форме Коши, получаемой ( x=x_1,\:\:\dot{x}=x_2 ) из уравнения (1.1):

В разделе, посвященном основам Динамической теории погрешностей, а также в работах [38,37,39] показано как знания [формулы: (4.10), (4.12), (4.13), (4.14)], полученные с помощью формализованного алгоритма семантического анализа ОДУ, можно использовать не только для решения уравнений по схеме Рунге-Кутты (2.8), но и для вычисления адекватности МОДЕЛЬ – МОДЕЛЬ.

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

Разностные уравнения

Положим, что лексический анализатор распознает, например, цепочку "x[n+1]" не как набор отдельных символов, а как значение переменной с именем "x" в момент времени [n+1]. Тогда, если принять стиль формирования имен, подобный тому, который мы рассмотрели для ОДУ, то есть добавлять к имени переменной, например, цепочки: "_minus2", "_minus1", "_0", "_plus1", и т.д., и если все анализируемые разностные уравнения (РУ) решены относительно младшего момента времени, легко переложить методику разбора ОДУ на системы разностных уравнений.

Графическое представление математических моделей

Легкость, с которой мы провели разбор уравнения (4.3), невозможна в случае анализа графических форм представления ММ. Здесь существует две проблемы. Первая состоит в том, что после анализа графических примитивов программа вычислений представляет собой набор операторов в хаотическом порядке, т.е. операторы не расположены в порядке их вычисления. Вторая состоит в том, что графическое представление систем ОДУ позволяет не изображать дважды операторы, относящиеся одновременно к двум или нескольким уравнениям системы. Поэтому мы будем рассматривать по крайней мере два прохода идентифицирующего планирования вычислений.

Суть идентифицирующего планирования вычислений поясним на примере следующей ММ.

Уравнение (4.15) в первом приближении можно считать, например, моделью бортовой качки корабля на регулярном волнении. Для представления данной модели в виде Алгоритмической Сетевой Модели (АСМ) [29], уравнение (4.15) необходимо преобразовать в систему двух ОДУ первого порядка. Введем обозначения \theta=x, \dot{\theta}=y и для простоты величину F\: cos\: \omega_k t обозначить как F\: cos\: \omega_k t=F. Тогда получим следующую систему ОДУ:

С помощью метода Эйлера данную модель можно преобразовать в разностную схему:

которой соответствует АСМ, представленная на рис. 4.3.

После анализа графического представления АСМ, в системе КОГНИТРОН [29], все операторы сети заносятся в список или матрицу операторов. При этом порядок следования операторов в данном списке произволен и зависит только от координат графических объектов на экране или последовательности занесения их пользователем в память ЭВМ. Для того, чтобы выстроить операторы сети в порядке необходимом для расчета в системах, подобных среде моделирования КОГНИТРОН, предусмотрен так называемый Планировщик вычислений, работа которого основана на том, что оператор может быть вычислен только тогда, когда получены значения переменных на его входах.

В результате планирования вычислений представление сети (рис.4.3) может выглядеть, например, следующим образом.

Здесь переменные b_{1}, b_{2}, b_{3} и b_{4} служат для передачи результатов между операторами и в случае АС называются внутренними переменными сети. Переменные 2\vartheta_{\theta}, {\eta}_{\theta}^2 и F_{n-1} являются для АС входными переменными.

Справа от списка (4.18) приведены номера уравнений системы (4.17), которым принадлежат операторы сети.

В случае применения в качестве оператора дискретизации (ОД) метода Эйлера нас не интересует взаимное расположение операторов в списке (4.18). Достаточно лишь соблюдения условия их вычислимости. Поскольку планировщик АСМ в принципе не продполагает, что исходные модели могут представлять собой ОДУ, планирование вычислений в системе КОГНИТРОН могло бы дать и другой результат, т.к. в случае (4.18) совершенно не важно взаимное расположение операторов, например, 7 и 8, 4 и 5, 5 и 6, 1 и 2. Возможный конечный результат планирования зависит только от последовательности “снятия” операторов с экрана графического редактора.

Алгоритм идентифицирующего планирования вычислений

Перепишем программу (4.18), перемешав операторы так, как они могли бы располагаться в списке до начала планирования вычислений.

Организуем сортировку операторов в списке (4.19) (планирование вычислений) не только в необходимом для их вычисления порядке, но и разделяя их на две группы, относящиеся к отдельным уравнениям системы (4.17).

Так как в списке (4.19) имеется два оператора \Delta t, создадим два пустых списка, в которые будем перемещать операторы, относящиеся к отдельным уравнениям системы (4.17). Поиск операторов будем осуществлять сверху вниз.

Для идентификации уравнений в первую очередь необходимо найти операторы \Delta t. В списке (4.19) первым оператором \Delta t является оператор под номером 2. Переместим его на “дно” первого списка, удалив его из списка (4.19). Входом этого оператора является переменная x_{n}. Найдем оператор, выходом которого является переменная x_{n} и переместим его из списка (4.19) в первый, созданный нами список над оператором \Delta t. Этим оператором является оператр под номером 3. В результате на местах 7 и 8 в новой программе окажутся операторы 2 и 3 списка (4.19).

Нашей задачей является выстроить над оператором \Delta t все операторы, работающие на его вход. Чтобы определить, найдены ли все такие операторы, проведем анализ входов оператора 7 в новом списке (4.20). Первым входом оператора 7 является переменная x_{n-1}, которая для данной АСМ является выходом оператора \Delta t, и ее значение определяется начальными условиями для модели (4.17). Следовательно, поиск операторов, влияющих на вход x_{n-1} оператора 7 можно прекратить. Проведя аналогичные рассуждения для входа y_{n-1}, можно сделать вывод о том, что составление группы операторов, идентифицирующих первое уравнение системы (4.17), закончено.

Перейдем к поиску следующих операторов \Delta t среди оставшихся в списке (4.19). Таковым является оператор 5. Переместим его на “дно” второго списка, удалив из списка (4.19). Найдем оператор, работающий на его вход. В старом списке это оператор 1. После перемещения его в новый список получим следующее состояние строящейся программы:

Здесь первым входом оператора 5 является переменная y_{n-1}, являющаяся выходом оператора \Delta t. Следовательно, поиск операторов, влияющих на данный вход прекращается.

Вторым входом оператора 5 служит переменная b_{4}, которая не является ни входной переменной сети, ни входом или выходом какого-либо оператора \Delta t. Следовательно, (поскольку список (4.19) еще не пуст) следующим шагом будет поиск оператора, выходом которого будет переменная b_{4}. Таковым является оператор 4 старого списка (4.19), переместим его на место 4 новой программы и проанализируем аналогичным образом его входы. Переменная F_{n-1} является входной для АС, а в формировании значения b_{3} участвует оператор 6 старого списка (4.19). Поэтому мы имеем следующее состояние строящейся программы:

Оператору 3 нового списка следует уделить большее внимание, т.к. при анализе первого входа (переменная b_{1}) программа планирования вычислений может “уйти” слишком далеко.

Для того, чтобы после прохождения всех операторов, влияющих на переменную b_{1} мы могли вернуться назад для анализа второго входа оператора 3 (переменная b_{2}), необходимо запомнить в стеке два числа: первое – номер разбираемого оператора, а второе – номер следующего входа, подлежащего анализу. В нашем случае в стеке мы должны сохранить числа 3 и 2.

Учитывая все сказанное приведем программу, полученную после выполнения идентифицирующего планирования вычислений по списку (4.19).

Графическое представление систем ОДУ в форме АСМ (рис.4.3) не требует изображать дважды операторы, относящиеся одновременно к двум или нескольким уравнениям системы и, поскольку к данному моменту список (4.19) оказывается пустым, может возникнуть необходимость копирования таких операторов по вновь созданным спискам из тех, в которые они попали первыми. Этот и ряд других нюансов говорят о том, что приведенный здесь алгоритм можно назвать лишь первым проходом полной процедуры идентифицирующего планирования вычислений в системах с графическим вводом ММ.

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

Второй проход идентифицирующего планирования

Второй проход процедуры идентифицирующего планирования вычислений рассмотрим на примере модели (1.7) [4], рассмотренной нами в части посвященной математической поддержке моделирования раздела Вопросы и цели. Модель (1.7), при использовании метода Эйлера, дает следующую разностную схему.

Модели (4.24) соответствует АСМ, представленная на рис.4.4.

После анализа графического представления ММ, на экране редактора, операторы программы могут быть перемешаны, например, так:

Тогда после первого прохода идентифицирующего планирования мы получим следующее внутримашинное представление программы.

Из программы (4.26) видно, что второе и третье уравнения системы (4.24) мы идентифицировали не полностью. Так значения переменных X5 и X12 (показанных серым цветом) в рамках одного уравнения вычислить невозможно, следовательно, невозможно составить функции расчета правых частей уравнений системы (1.7). Именно для разрешения таких ситуаций и служит второй проход идентифицирующего планирования.

Алгоритм второго прохода весьма прост. Необходимо проверить вычислимость списков, представляющих отдельные уравнения системы, и, вслучае необходимости, продублировать недостающие операторы. Так для программы (4.26) мы должны сделать следующее.

1. Скопировать оператор 9 из списка 1 в начало списка 2.
2. Скопировать затем подсписок 2, состоящий из операторов: 9, 20, 1, 5, 2, 6, 8, 19, в начало списка 3.

Тогда окончательный вариант внутреннего представления ММ будет выглядеть следующим образом.

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

Уравнения, неразрешенные относительно старшей производной

Рассмотрим модель вертикальной (\zeta) и килевой (\psi) качки корабля.

В практике проектирования и исследований принимается [44], что вертикальная и килевая качка корабля на волнении описывается следующей системой уравнений:

где F(t), M(t) – возмущающие силы и моменты, а коэффициенты модели определяются сложными, в том числе и нелинейными, зависимостями от многих параметров корабля и окружающей среды.

Запишем систему (4.28) в виде линейной последовательности байт, заменив \zeta и \psi на x и y, соответственно, и переписав ее в несколько ином виде:

  x'' = (F - B_xx*x' - C_xx*x - A_xy*y'' - B_xy*y' - C_xy*y) / A_xx;
  y'' = (M - B_yy*y' - C_yy*y - A_yx*x'' - B_yx*x' - C_yx*x) / A_yy;

Если разрешить наличие в правой части уравнений производных второго порядка и провести синтаксический и семантический анализ данных уравнений по грамматике ОДУ (рассмотренной нами ранее), то после первого прохода идентифицирующего планировщика мы получим следующее внутреннее представление программы:

Алгоритм второго прохода идентифицирующего планирования вычислений должен распознать сложившуюся ситуацию (см. программу (4.29)), чтобы не копировать операторы из уравнения 2 в уравнение 4 и наоборот до бесконечности. Распознать такую ситуацию довольно легко. Подсказкой здесь может служить то обстоятельство, что переменные _$tmp$_11 = _2p_x и _$tmp$_22 = _2p_y являются входными для интеграторов и символизируют старшую производную, значение которой не должно входить в список начальных условий задачи Коши.

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

Аналитически, приведение системы (4.28) к системе четырех уравнений первого порядка возможно по алгоритму, предложенному Абгаряном К.А. [13]. Поэтому задача автоматизации приведения систем уравнений к нормальной форме вполне разрешима. Однако если компилятор языка моделирования не имеет соотвестствующих средств, то в нем, по крайней мере, должна быть реализована возможность распознать такие уравнения и предупредить пользователя об ошибке.

Схемы для аналоговых вычислитеьных машин

Для того чтобы переложить все сказанное выше на язык схем АВМ [41], необходимо иметь возможность обеспечения адекватности преобразований аналоговых интеграторов в дискретные операторы ЦВМ (см. рис. 4.2). Поскольку задача идентификации исходных математических моделей, на этапе семантического анализа, практически разрешима для систем ОДУ, мы можем говорить о том, что Динамическая теория погрешностей численных методов может найти широкое применение не только для анализа результатов вычислительных экспериментов, но и при создании компиляторов языков моделирования.

Следует заметить, что алгоритм идентифицирующего планирования вычислений можно распространить не только на схемы для АВМ, но и схемы в переменных состояния [56], а также на структурные схемы систем автоматического управления [56,32].

Структурные схемы систем автоматического управления

Перепишем схему [32], показанную на рис. 1.6, в виде, более удобном для интерпретации:

Заменив идеальные интеграторы (1/s) на дискретные операторы согласно схеме, показанной на рис.4.2, а переменные \theta и \beta на x и y, соответственно, затем, "сняв" все операторы с изображения рис.4.5, и осуществив идентифицирующее планирование вычислений, получим программу:

Если теперь в программе (4.30) произвести замену имен переменных:

мы получим следующую разностную схему:

Тот же самый результат можно получить, приведя систему (1.9) к нормальной форме Коши (произведя замену: \theta=x_1, \dot{\theta}=x_2, \beta=x_3, \dot{\beta}=x_4) и применив к ней метод Эйлера.

Оптимизация кода

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

Выход здесь достаточно прост.

Если входными переменными какого-либо оператора не являются переменные состояния модели (вход, выход \Delta t) и переменные, значения которых вычисляются в блоке инициализации шага интегрирования в зависимости от текущего времени (расчет коэффициентов), то данный оператор можно отнести к списку однократно вычисляемых. Следовательно, выходную переменную данного оператора также можно отнести к разряду однократно вычисляемых. Двигаясь таким образом к корню дерева вычислений (и помечая переменные, в таблице символов, соответствующим образом) можно выделить поддерево (если, разумеется, оно существует), которое целиком можно вынести в процедуру инициализации модели. Подобные алгоритмы давно разработаны в теории алгоритмического сетевого моделирования [29], однако, не нашли еще практического применения для оптимизации кода.