авторефераты диссертаций БЕСПЛАТНАЯ БИБЛИОТЕКА РОССИИ

КОНФЕРЕНЦИИ, КНИГИ, ПОСОБИЯ, НАУЧНЫЕ ИЗДАНИЯ

<< ГЛАВНАЯ
АГРОИНЖЕНЕРИЯ
АСТРОНОМИЯ
БЕЗОПАСНОСТЬ
БИОЛОГИЯ
ЗЕМЛЯ
ИНФОРМАТИКА
ИСКУССТВОВЕДЕНИЕ
ИСТОРИЯ
КУЛЬТУРОЛОГИЯ
МАШИНОСТРОЕНИЕ
МЕДИЦИНА
МЕТАЛЛУРГИЯ
МЕХАНИКА
ПЕДАГОГИКА
ПОЛИТИКА
ПРИБОРОСТРОЕНИЕ
ПРОДОВОЛЬСТВИЕ
ПСИХОЛОГИЯ
РАДИОТЕХНИКА
СЕЛЬСКОЕ ХОЗЯЙСТВО
СОЦИОЛОГИЯ
СТРОИТЕЛЬСТВО
ТЕХНИЧЕСКИЕ НАУКИ
ТРАНСПОРТ
ФАРМАЦЕВТИКА
ФИЗИКА
ФИЗИОЛОГИЯ
ФИЛОЛОГИЯ
ФИЛОСОФИЯ
ХИМИЯ
ЭКОНОМИКА
ЭЛЕКТРОТЕХНИКА
ЭНЕРГЕТИКА
ЮРИСПРУДЕНЦИЯ
ЯЗЫКОЗНАНИЕ
РАЗНОЕ
КОНТАКТЫ


Pages:     | 1 | 2 || 4 | 5 |   ...   | 6 |

«МОСКОВСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ ИМЕНИ М.В. ЛОМОНОСОВА Факультет вычислительной математики и кибернетики ...»

-- [ Страница 3 ] --

Сборник статей молодых ученых факультета ВМК МГУ, № 10, 120 А. В. Мальцева В общем случае при поиске областей устойчивости для аф финных полиномов необходимо учитывать соотношения меж ду размерностью пространства параметров, размерностью про странства коэффициентов аффинного полинома и рангом мат рицы A линейного преобразования.

a q a.

. = A..

.

..

qk an 4.1 Построение внутренней аппроксимации.

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

Теорема 7. Пусть A : Rk Rn+1 линейное преобразование, удовлетворяющее условиям:

1. k n + 1;

2. rg A = n + 1.

параллелотоп из Rk ) Тогда Q = A (Q) (Q невырожденный зонотоп.

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

1. k n + 1;

2. матрица A имеет максимальный ранг.

где k размерность пространства параметров, n + 1 размер ность пространства коэффициентов аффинного полинома.

Способ внутренней аппроксимации параллелотопа Q пря мыми параллелотопами состоит в следующем:

Сборник статей молодых ученых факультета ВМК МГУ, № 10, Метод внутренней аппроксимации 1. Преобразованием A параллелотоп Q переводится в зоно топ Q, dim Q = n + 1. При помощи V -преобразования получаем параллелотоп P такой, что Q P. Далее бе рем покрытие {Pi } параллелотопа P (P = Pi, Pi па раллелотопы с заданной шириной ).

2. Из покрытия {Pi } выбираем максимальное подпокрытие {Pi } Q. Для реализации представим п. 2 представим зонотоп Q как множество решений некоторой системы линейных неравенств.

Алгоритм построения системы линейных неравенств, за дающих зонотоп Q (в Rn ):

(a) Перебор всевозможных наборов по n точек из мно жества образов вершинных точек параллелотопа Q при преобразовании A.

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

лежат ли они по одну сторону от нее).

(c) Отбрасывание гиперплоскостей, для которых точки лежат по обе стороны.

(d) Определение знаков неравенств.

В результате получаем систему линейных нера венств F1 (x) F2 (x) (8),.

.

.

Fn (x) задающую зонотоп Q.

По вершинным точкам каждого параллелотопа из {Pi } определяется принадлежность множеству Q (подстановкой в систему (8)).

Сборник статей молодых ученых факультета ВМК МГУ, № 10, 122 А. В. Мальцева 4.2 Сходимость метода внутренней аппроксима ции.

Каждый параллелотоп множества {Pi } проверяется на устойчивость/неустойчивость и, в зависимости от этого, либо заносится в множество S или N, либо в U и к нему применя ется описанный выше алгоритм SIVIA.

Применение алгоритма SIVIA к параллелотопам {Pi } поз воляет получить разбиение {Pi } = S U N, которое, в свою очередь, порождает соответствующее разбиение множества па раметров Q, т. е. Q = S U N. При этом верна следующая теорема.

Теорема 8. Пусть A : Rk Rn+1 линейное преобразование, удовлетворяющее условиям:

1. k n + 1;

2. rg A = n + 1.

параллелотоп в Rk и Q = A (Q). Тогда, если для Пусть Q последовательности вложенных множеств Q U1 U... Um... µ (Um ) 0, то µ (Um ) 0, где m m A (Um ) = Um.

Доказательство теоремы 8 основывается на свойствах ли нейных преобразований в евклидовых пространствах.

Теорема 8 фактически утверждает сходимость метода внут ренней аппроксимации.

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

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

Пусть q S, а V Q множество решений системы Ax =. Тогда любой q V стабилизирующий набор параметров.

q Сборник статей молодых ученых факультета ВМК МГУ, № 10, Метод внутренней аппроксимации 5 Применимость метода внутренней ап проксимации для задачи одновременной стабилизации При решении задачи одновременной стабилизации [7] воз никает необходимость поиска областей устойчивости аффин ных полиномов (s;

q) = 0 (q) + 1 (q) s +... + n+k1 (q) sn+k1, задающих преобразование пространства параметров в про странство коэффициентов qk.

... b1...

1 0 0 0 0 0 0 0 0 0.

n+k q. a1... b0 b1...

1 0 0 0 0 0 0 0 a. a1... b0 b1...

1 0 0 0 0 0 0 =... pk,....

.....................................

.

... a0 a1... b0 b 0 0 0 0 0 0 0.

.

... a0... b 0 0 0 0 0 0 0 0 0 p где ai, bi, pi, qi коэффициенты полиномов (s) = sn + an1 sn1 +... + a0, (9) n (10) (s) = bn1 s +... + b0, p (s) = pk sk + pk1 sk1 +... + p0, q (s) = sk +... + q0.

При этом (s), (s) взаимно простые, заранее заданные по линомы и k n + 1. Тогда верна следующая теорема.

Теорема 9. Пусть (s), (s) взаимно простые полино Сборник статей молодых ученых факультета ВМК МГУ, № 10, 124 А. В. Мальцева n + 1. Тогда матрица мы (9, 10), k 0... an an1 an2...... a 0... an an1 an2... a1 a.........

..................

0 a1 a... an an1 an2...

A= bm... bm1... b0......

0... bm... b1 b0...

.........

..................

0......... bm... b1 b имеет полный ранг.

Доказательство. Пусть даны многочлены n ai xi = a0 + a1 x +... + an xn, (x) = i= m bi x i = b0 + b 1 x +... + bm x m.

(x) = i= Тогда матрицей Сильвестра для этих многочленов будет квад ратная матрица (n + m) (n + m) вида 0... an an1 an2...... a 0... an an1 an2... a1 a............

...............

0 a1 a... an an1 an2...

S, = bm bm1...... b0......

0... bm... b1 b0...

............

...............

0......... bm... b1 b Количество строк матрицы, содержащих коэффициенты мно гочлена (x) равно m, а многочлена (x) n.

Сборник статей молодых ученых факультета ВМК МГУ, № 10, Метод внутренней аппроксимации bm1 sm1 +...+b1 s+b, W (s) = sm +am1 sm1 +...+a1 s+a pl sl +pl1 sl1 +...+p1 s+p R (s) = sl +q sl1 +...+q s+q 1 l (s) = sm + am1 sm1 +... + a0 sl + ql1 sl1 +... + q0 + + bm1 sm1 + bm2 sm2 +... + b0 pl sl + pl1 sl1 +... + p sm+l1 : ql1 + am1 + bm1 pl sm+l2 : ql2 + am2 + am1 ql1 + bm1 pl1 + bm1 pl...

s2 : a2 q0 + a0 q2 + b2 p0 + b0 p2 + a1 q1 + b1 p s1 : a1 q0 + a0 q1 + b1 p0 + b0 p s 0 : a 0 q 0 + b0 p m+l1... b1...

1 0 0 0 0 0 0 0 0 0 0 m+l a1... b0 b1...

1 0 0 0 0 0 0 0.

......

.................................................

= b.... a0 a1 a2... b0 b 0 0 0 0 0 0 0 0 2... a0 a1... b0 b 0 0 0 0 0 0 0 0 0 0 0... a0... b 0 0 0 0 0 0 0 0 0 0 0 0 0 ql.

.

q. pl, A : R2l+1 Rl+m, rg A = l + m.

.

.

.

p Вычеркиваем l m + 1 столбцов в зависимости от b0 и a0.

1. Если b0 = 0, 1 b1 det a1 b0 b1.

lm+ det (...) = b a 0 0 b 2. Если b0 = 0 a0 = 0, 1 b1 det a1 b0 b1.

lm+ det (...) = a a 0 0 b Таким образом, при решении задачи одновременной стаби лизации выполнены основные предположения применимости метода внутренней аппроксимации для поиска областей устой чивости аффинных полиномов.

Сборник статей молодых ученых факультета ВМК МГУ, № 10, 126 А. В. Мальцева Заключение Таким образом, в представленной работе 1. Вскрыты механизмы отсутствия сходимости метода внешней аппроксимации.

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

Список литературы [1] Коровин С. К., Кудрицкий А. В., Фурсов А. С. Конструк тивный алгоритм поиска регулятора, одновременно ста билизирующего семейство объектов // Нелинейная дина мика и управление: Сборник статей, Вып. 7 / Под ред.

С. В. Емельянова, С. К. Коровина, - М.: ФИЗМАТЛИТ, 2010, с. 5–16.

[2] Коровин С. К., Фурсов А. С. Одновременная стабилиза ция: синтез универсального регулятора // Автоматика и телемеханика, 2011, № 9, с. 61–73.

[3] Г. Алефельд, Ю. Херцбергер Введение в интервальные вы числения // М.: Мир, 1987. 360 с.

[4] Л. Жолен, М. Кифер, О. Дидри, Э. Вальтер Прикладной интервальный анализ // М. Ижевск: Институт ком пьютерных исследований, 2007 468 с. ISBN 5–93972– 585–6.

[5] Б. Поляк, П. Щербаков Робастная устойчивость и управ ляемость // М. Наука, 2002.

Сборник статей молодых ученых факультета ВМК МГУ, № 10, Метод внутренней аппроксимации [6] Емельянов С. В., Коровин С. К., Фомичев В. В., Фур сов А. С. Задачи и теоремы по теории линейной обрат ной связи // М.: Издательский отдел ф-та ВМиК МГУ им. М. В. Ломоносова, 2004.

[7] А. С. Фурсов Одновременная стабилизация: теория по строения универсального регулятора для семейства ди намических объектов // Диссертация. Москва. 2012.

Сборник статей молодых ученых факультета ВМК МГУ, № 10, 128 Сборник статей молодых ученых факультета ВМК МГУ, № 10, УДК 004. Библиотечная реализация стековой вычислительной модели на примере языка Joy © 2013 г. А. В. Мошкина a.v.moshkina@gmail.com Кафедра алгоритмических языков 1 Введение Перед началом решения некоторой задачи, программист естественным образом сталкивается с вопросом выбора языка программирования, на котором будет написан проект. Часто этот вопрос состоит не столько в выборе языка, сколько в вы боре парадигмы программирования, в рамках которой удобнее всего решать поставленную задачу.

Понятие парадигмы программирования не имеет однознач ного определения. В рамках данной работы под парадигмой программирования понимается система идей и понятий, кото рые определяют ход мышления программиста и стиль написа ния программ. Язык программирования не всегда однозначно определяет парадигму, используемую при работе с этим язы ком. Согласно [1], язык программирования может иметь раз личное отношение к конкретной парадигме: от навязывания применения парадигмы в рамках языка, до полного запреще ния применения парадигмы. Поэтому языку могут соответ ствовать сразу несколько парадигм, в разной степени харак терных ему.

Однако, парадигма, которую удобно использовать в одной задаче, не всегда может подойти для других, в то время как в рамках одного проекта могут решаться множество разнород ных задач. Например, любое полноценное приложение должно Библиотечная реализация языка Joy иметь графический интерфейс пользователя, который удобно проектировать в рамках объектно-ориентированной парадиг мы. Почти каждое приложение осуществляет взаимодействие с базами данных, запросы к которым носят декларативный ха рактер в запросе описываются свойства требуемых данных, и никак не оговариваются способы их получения. Кроме того, в процессе обработки и манипулирования этими данными мо жет быть удобно применение императивной, функциональной или, например, сценарной парадигм программирования.

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

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

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

Метод непосредственной интеграции успешно применялся для расширения языка С++ возможностями различных аль тернативных языков [3–5], что позволяет говорить о перспек тивности этого подхода.

Работа посвящена реализации метода непосредственной ин теграции для стекового функционального языка Joy.

Сборник статей молодых ученых факультета ВМК МГУ, № 10, 130 А. В. Мошкина 2 Стековые языки программирования Существует множество классификаций языков программи рования, и одна из них разделение языков на конкатенатив ные и аппликативные. Идея аппликативных языков заключа ется в том, что вычисление программы производится посред ством применений функций к их аргументам. К этой группе можно отнести большое число языков общего назначения, на пример, C, Java, Python, Haskell и другие. В конкатенативных языках программирования вычисление программы состоит в композиции нескольких функций (то есть последовательного применения одной функции к результату вычисления преды дущей функции), причем все они оперируют с одним и тем же объектом данных, который передается от функции к функ ции. Чаще всего этим объектом данных становится стек. Также важно заметить, что описанная композиция функций реализу ется простой конкатенацией программного кода этих функций.

К семейству конкатенативных языков можно отнести такие языки программирования, как Forth [6], Joy [7], PostScript [8], Cat [9], Factor [10] и другие. Все они имеют сильное семейное сходство, однако значительно различаются по дизайну, реали зации и назначению.

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

Особенностью стековых языков программирования являет ся то, что в них присутствуют сразу несколько стеков. Один Сборник статей молодых ученых факультета ВМК МГУ, № 10, Библиотечная реализация языка Joy из них это стек вызовов (стек возвратов или действий), ко торый используется для реализации рекурсии, но кроме него существует также стек данных (стек значений), предназначен ный для передачи данных между функциями. Стек данных используется при программировании чаще, поэтому в случаях, когда это не вызывает недоразумений, он называется просто стеком.

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

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

считать, что любая функция имеет ровно один аргумент состояние стека, и возвращает в качестве значения новое состояние стека.

В стековых языках программирования проще осуществля ется возврат нескольких значений из функции, чем в апплика тивных языках. В последних это осуществляется посредством объединения отдельных элементов данных в кортежи (напри мер, списки для языка Lisp или структуры для языка C). В то время как в стековых языках это осуществляется последова тельным помещением на вершину стека нескольких элементов Сборник статей молодых ученых факультета ВМК МГУ, № 10, 132 А. В. Мошкина данных. Это исключает накладные расходы по предваритель ной упаковке данных и последующей распаковке перед их ис пользованием.

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

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

3 Язык Joy Стековый язык программирования Joy был разработан в 2001 году Манфредом фон Таном в университете Мельбур на, Австралия. Joy представляет собой чисто функциональный язык программирования: в программе на Joy отсутствует по нятие состояния, и, как следствие, отсутствуют переменные и присваивание.

Типы данных, поддерживаемые языком, можно разделить на две группы атомарные и агрегатные. Атомарные типы включат символы, целые числа, числа с плавающей точкой, логические величины. К агрегатным типам относятся строки, множества, списки, файлы и квотированные программы. Под держиваются только числовые множества, содержащие целые значения в диапазоне от 0 до 31. Список является гетероген ной структурой данных, его элементами могуть быть данные Сборник статей молодых ученых факультета ВМК МГУ, № 10, Библиотечная реализация языка Joy любых типов, как атомарных, так и агрегатных, в том числе и другие списки. Квотированная программа это расширение понятие списка: наряду с данными, квотированная программа может также содержать функциональные вызовы. Элементы списков и квотированных программ записываются в квадрат ных скобках и разделяются пробелами, например:

[2 ["string" literal] ’c’ {7 5 3} []] [[1 2 3 4] [dup *] map] Программа записывается в обратной польской нотации, в которой знак операции записывается после ее операндов.

Например, программа, вычисляющая арифметическое выра жение (2 + 3) · 15, будет выглядеть следующим образом:

2 3 + 15 *.

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

Деквотирование совершает обратное действие активизирует квотированный программный код, делая его готовым к немед ленному исполнению. Таким образом, выявляется концепция единства данных и программного кода, которая впервые была реализована Джоном Маккарти в языке Lisp [11].

4 Примеры некоторых комбинаторов язы ка Joy Функции, осуществляющие в каком–либо виде деквоти рование называются комбинаторами. С помощью комбинато ров осуществляется управление потоком программы (ветвле ние, рекурсия), а также с помощью комбинаторов реализуются функции высших порядков.

Сборник статей молодых ученых факультета ВМК МГУ, № 10, 134 А. В. Мошкина Для начала рассмотрим каким образом в Joy записываются определения функций.

DEFINE square == dup * ;

cube == dup dup * * END Здесь представлены определения двух функций. Функция square вычисляет квадрат своего числового аргумента. Для этого она сначала удваивает аргумент, лежащий на вершине стека функцией dup, а затем функция * достает с вершины стека два верхних элемента и возвращает на стек их произ ведение. Аналогично, функция cube возводит свой аргумент в третью степень.

Блок определений функций начинается с ключевого слова DEFINE и заканчивается словом END. Внутри блока может на ходится одно или несколько определений. Между названием и телом функции ставится разделитель ==. В случае нескольких определений, они отделяются друг от друга точкой с запятой.

Язык Joy не предоставляет специального синтаксиса управ ляющих конструкций программы. Как уже было сказано, управление программой производится с помощью комбинато ров. Например, комбинатор ifte реализует ветвление. Функ ция ожидает на вершине стека 3 параметра, каждый из кото рых является квотированной программой. Рассмотрим пример:

DEFINE factorial == [0 =] [pop 1] [dup 1 - factorial *] ifte END Здесь представлено рекурсивное определение функции, счита ющей факториал целого числа. Третий аргумент комбинатора Здесь под словом понимается лексическая единица, а не специфич ное для стековых языков название функции. Ключевые слова не являются функциями, это всего лишь метки.

Сборник статей молодых ученых факультета ВМК МГУ, № 10, Библиотечная реализация языка Joy (if-part) проверяет аргумент на равенство нулю. Если усло вие выполнено, то деквотируется второй агрумент (then-part) и на вершину стека попадает ответ (0! = 1). Если условие не выполнено, то деквотируется первый аргумент (else-part), в котором происходит рекурсивный вызов для аргумента, умень шенного на единицу (n! = (n 1)! · n).

Для реализации рекурсии Joy предоставляет большое число рекурсивных комбинаторов. Рассмотрим два из них комбинатор линейной рекурсии linrec и комбинатор бинарной рекурсии binrec. Комбинатор linrec ожидает на вершине сте ка 4 параметра, каждый из них квотированная программа:

[if-part] [then-part] [rec1-part] [rec2-part] linrec.

Четвертый параметр проверяет условие выхода из рекурсии и, если оно выполняется, то выполняется нерекурсивная часть определения. Если условие не выполняется, то сначала выполняется rec1-part, затем происходит рекурсивный вызов для только что вычисленного значения, после возврата из которого выполняется rec2-part. В следующем примере дается определение факториала с помощью linrec:

DEFINE factorial1 == [null] [succ] [dup pred] [*] linrec END Полиморфная функция null возвращает true, если ее ар гументом является нуль (целого типа или с плавающей точкой) или пустой объект агрегатного типа (строка, список или мно жество), а каждая из функций pred и succ уменьшает или, со ответственно, увеличивает значение целого аргумента на еди ницу.

Комбинатор бинарной рекурсии совершает два рекурсив ных вызова. Комбинатор binrec ожидает на вершине стека 4 аргумента-квотированные программы. Почти все они име ют тот же смысл, что и для комбинатора linrec кроме rec1-part. Теперь результатом выполнения этой квотирован ной программы должны быть два значения, для каждого из которых происходит рекурсивный вызов.

Сборник статей молодых ученых факультета ВМК МГУ, № 10, 136 А. В. Мошкина Комбинаторы можно использовать не только в определе ниях пользовательских функций. Следующий фрагмент про граммного кода вычисляет восьмой элемент последовательно сти Фибоначчи: 8 [small] [] [pred dup pred] [+] binrec.

Полиморфная функция small возвращает true, если ее целый числовой аргумент не больше единицы или, если объект агре гатного типа содержит не более одного элемента.

Рассмотрим несколько комбинаторов, которые будут ис пользоваться в примерах (схема вида A1...An - B1...Bn обозначает стековый эффект функции Ai обозначают аргу менты функции, а Bi обозначают ее возвращаемые значения):

• uncons: list - head tail разделяет список на голо ву и хвост, первым кладет в стек голову, поверх нее хвост;

• split: qprog list - listX listY разбивает список на два, которые возвращает на стек. В первый (listX) попадают элементы, для которых выполнение квотиро ванной программы (qprog) дает ненулевой результат (или true), из всех остальных формируется второй список (listY);

• cons: elem list - list1 помещает элемент в начало списка;

• dip: qprog param -... первым параметром являет ся квотированнная программа. Второй параметр убира ется со стека и сохраняется, затем выполняется квоти рованная программа, после чего сохраненный параметр помещается на вершину стека.

Более интересным примером использования комбинатора бинарной рекурсии является описание функции быстрой сор тировки:

DEFINE Сборник статей молодых ученых факультета ВМК МГУ, № 10, Библиотечная реализация языка Joy qsort == [small] [] [[uncons] [] split] [[swap] dip cons concat] binrec END 5 Библиотека InteLib и ее структура Проект InteLib [12] представляет собой библиотеку, реали зующую метод непосредственной интеграции, где в качестве базового выбран объектно–ориентированный язык C++ [13].

Язык C++ обладает механизмом переопределения символов стандартных операций для пользовательских типов данных, что является ключевым моментом в выборе базового языка программирования для реализации метода непосредственной интеграции.

Рассмотрим понятие символьного выражения (или S выражения). S-выражение представляет собой либо атомар ные данные (символ, числовая константа, строковый литерал и т. п.) либо точечную пару вида (A. B), где A и B также являют ся S-выражениями. Наряду с понятием точечной пары вводят ся понятия списка и пустого списка. Список это точечная па ра вида (A. (B. (C. NIL))), где NIL специальное обозна чение для пустого списка. Таким образом, на основе типизиро ванных атомарных S-выражений можно строить нетипизиро ванные контейнеры, являющиеся составными S-выражениями, и содержащими в себе разнородные данные.

Над S-выражениями определены операции композиции (со ставление точечной пары из двух заданных S-выражений), де композиции (получение каждого из элементов точечной пары) и вычисления S-выражения (некоторое отображение множе ства S-выражений на само себя).

Библиотека InteLib реализована в виде нескольких логи ческих слоев. Нижний слой это S-выражения, не связан ные ни с какими вычислительными моделями. Символьны Сборник статей молодых ученых факультета ВМК МГУ, № 10, 138 А. В. Мошкина ми выражениями в библиотеке InteLib представляются дан ные и программный код. S-выражения реализованны в ви де иерархичного дерева наследования с общим предком. Кор нем иерархии является класс SExpession, представляющий наиболее общие свойства всех возможных S-выражений, от него наследуются классы, представляющие конкретные S выражения, такие как SExpressionInt, SExpressionString, SExpressionChar, SExpressionCons и другие.

Для удобства использования класса SExpession был разра ботан класс SReference, являющийся интерфейсом для эффек тивной работы с S-выражениями. Этот класс инкапсулирует S выражение, его поведение аналогично поведению указателя на S-выражение, кроме того, SReference может реализовывать и некий дополнительный функционал (например, счетчик ссы лок на S-выражение, позволяющий организовать сборку мусо ра).

Верхний слой библиотеки реализует вычислительные мо дели S-выражений для различных языков, таких как Lisp, Scheme [14], Refal [15]. Кроме того, в библиотеке представлен средний слой, предоставляющий некоторое инструментальные средства для работы с S-выражениями (средства ввода-вывода, синтаксический анализ текстового представления и т. п.), а так же инкапсулирующий некоторые общие части вычислитель ных моделей для S-выражений. В ходе реализации библиотеки активно используются возможности определения абстрактных типов данных в языке С++, наследования и полиморфности объектов.

Идея реализации общей части для вычислительных моде лей в библиотеке InteLib была навеяна механизмом продолже ний, используемым в языке Scheme. Продолжение представ ляет собой последовательность действий, которые необходимо выполнить до завершения вычисления [3]. Результатом выпол нения этих действий являются некоторые значения, которые могут использоваться для совершения последующих действий.

Сборник статей молодых ученых факультета ВМК МГУ, № 10, Библиотечная реализация языка Joy Поэтому речь идет о двух основных сущностях: о совершаемых действиях и результатах их значений. При этом чем позднее было получено значение выражения, тем быстрее оно понадо бится для вычислений, что естественно приводит нас к концеп ции двух стеков стека значений и стека действий.

Библиотека InteLib предоставляет все необходимые инстру менты для библиотечной реализации языка Joy, а именно S выражения для реализации гетерогенных списков и вычисли тельную модель на основе продолжений для реализации кон катенативности.

6 Реализация вычислительной модели Вычислительная модель языка очень проста. После разби ения на лексемы текст программы представляет собой после довательность слов-функций и значений произвольного типа (атомарного или агрегатного), разделенных пробелами. Текст программы просматривается последовательно. Если текущий элемент является значением, то это значение просто кладется в стек, а если на очередном шаге встречается имя функции, то из стека извлекаются аргументы функции (в соответствии с ее стековым эффектом), далее происходит вычисление функции, после чего результат попадает на вершину стека.

Единственная конструкция, вычисление которой не соот ветствует вышеописанному алгоритму, это объявление поль зовательской функции. При вычислении блока определений, имя и тело определяемой функции не вычисляются, а заносят ся в специальную таблицу функций, откуда они будут доступ ны для дальнейшего использования. Такая таблица реализует ся по аналогии со словарем языка Forth.

Таким образом, вычисление программы не всегда проис ходит по одному сценарию. Вычислитель может находиться в двух состояниях: определение пользовательской функции и интерпретация программного кода. Реализуется такое поведе Сборник статей молодых ученых факультета ВМК МГУ, № 10, 140 А. В. Мошкина ние с помощью конечного автомата. Для удобства определению пользовательской функции соответствует не одно, а два состо яния. В первом из них происходит считывание и запоминание имени функции, а во втором последовательное считывание и запись элементов тела функции. Диаграмма состояний вычис лителя представлена на рис. 1. Переход по дугам, помеченным звездочками, происходит по имени функции (стандартной или уже определенной пользовательской) или по значению. Если в некотором состоянии вычислителю встречается лексема, по которой нет переходов в другие состояния, то такая ситуация является исключительной.

Рис. 1. Диаграмма состяний вычислителя Для реализации вычислительной модели был использован класс IntelibContinuation библиотеки InteLib. Этот класс ре ализует коцепцию продолжений (континуаций) и инкапсулиру ет два стека стек значений и стек действий и предоставля ет инструменты для работы с ними. На очередном шаге работы континуации из стека действий извлекается верхний элемент инструкция, предписывающая континуации выполнение опре Сборник статей молодых ученых факультета ВМК МГУ, № 10, Библиотечная реализация языка Joy деленного функционала. Данные, получаемые в ходе этих дей ствий, помещаются в стек результатов.

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

Вычислительная модель языка Joy реализуется клас сом JoyContinuation, являющимся наследником классa IntelibContinuation. В континуацию были добавлены две инструкции: joy_evaluate и joy_dequote. По инструкции joy_evaluate, примененной к очередной встретившейся лек семе, континуация производит вычисление этой лексемы в соответствии с предложенной выше вычислительной моде лью. Инструкция joy_dequote активизирует квотированный программный код, такая иструкция может быть применена только к квотированной программе. Технически эта инструк ция разбирает список на элементы и для каждого элемента списка применяет инструкцию joy_evaluate.

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

(A. (B. (C. NIL))), причем пустой список не является то чечной парой, а представляет собой атомарное S-выражение.

Таким образом, списки и квотированные программы языка Joy реализуются классом JExpressionQProgram, являющимся на следником класса SExpressionCons, а пустой список реализует ся отдельным объектом типа SExpressionLabel. Указатель на квотированную программу реализуется специальным классом JReferenceQProgram, унаследованным от класса SReference.

Этот класс реализует операцию деквотирования списка (дру гими словами, операцию его вычисления), которая представле на методом Evaluate. Метод Evaluate производит вычисление списка с помощью континуации, используя описанную выше инструкцию joy_dequote. Операция деквотирования списков используется при реализации библиотечных комбинаторов, а также при написании кода на языке Joy непосредственно в ко Сборник статей молодых ученых факультета ВМК МГУ, № 10, 142 А. В. Мошкина де C++.

Далее вводится класс JoyQProgramConstructor. В нем пере гружается операция побитового или |, которая конструирует и возвращает указатель на список из одного элемента (то есть, объект класса JReferenceQProgram), этот элемент единствен ный аргумент перегруженной внутри класса операции |. Для создания пустого списка перегружается операция побитового отрицания. Для добавления нового элемента в уже существу ющий список, в классе JReferenceQProgram переопределяется операция запятая, позволяющая максимально приблизить за пись списка целевого языка констукциями языка C++. Ниже представлен пример создания списка [1 2 3 4] и пустого спис ка []:

JoyQProgramConstructor J;

(J| 1, 2, 3, 4);

~J;

Множества языка Joy могут содержать только целые чис ловые значения в диапазоне от 0 до 31. Множество реализуется классом JExpressionSet, наследником SExpression. Для кон струирования множеств, так же, как и для списков вводятся классы JoySetConstructor и JReferenceSet. Ниже приводят ся примеры создания множеств {1 2 3 4} и {}:

JoySetConstructor S;

(S| 1, 2, 3, 4);

~S;

Для реализации специальных лексем в библиотеке InteLib существует класс SExpressionSpecialToken, являющийся на следником класса SExpressionLabel. Этот класс имеет аб страктный метод Evaluate, который переопределяется в на следнике класса, задавая специальную логику вычисления ре ализуемой лексемы.

Класс SExpressionSpecialToken используется для реали зации библиотечных функций языка Joy. От него наследуется класс JFunctionExpression, являющийся абстрактным пред ставлением функции. В нем реализуется некоторый дополни тельный функционал, однако метод Evaluate в нем остается Сборник статей молодых ученых факультета ВМК МГУ, № 10, Библиотечная реализация языка Joy абстрактным, он переопределяется в наследниках класса, реа лизующих конкретные библиотечные функции.

Поскольку при программировании непосредственно S выражения не используются, по аналогии со списками и множествами, для использования функций вводится класс JReferenceFunction, представляющий указатель на функцию.

Конструктор класса принимает в качестве аргумента указа тель на экземпляр класса функции. Например, следующий код определяет указатель на функцию MUL, реализуемую классом JFunctionMul:

JFunctionReference MUL(new JFunctionMul());

Для записи ключевых слов (DEFINE, END) и разделителей (;

, ==) используется класс JToken, конструктор которого при нимает на вход текстовое представление слова (разделителя).

Из-за требований к идентификаторам языка C++, а также для единообразия, вместо разделителей ;

и == используются соот ветственно ключевые слова SEMICOLON и EQUALS. Для записи имен функций при определении (и при последующем вызове) используется класс JLabel, конструктор которого так же при нимает на вход строку с именем функции.

Вспомним уже встречавшийся пример определения функ ций square и cube на языке Joy и посмотрим как аналогичный код будет выглядеть на языке C++.

// Required declarations JoyQProgramConstructor J;

// Declarations of keywords JToken DEFINE("DEFINE");

JToken EQUALS("EQUALS");

JToken SEMICOLON("SEMICOLON");

JToken END("END");

// Declarations of names of user’s functions JLabel square("square");

JLabel cube("cube");

Сборник статей молодых ученых факультета ВМК МГУ, № 10, 144 А. В. Мошкина // Declarations of library functions JFunctionReference DUP(new JFunctionDup());

JFunctionReference MUL(new JFunctionMul());

// Code of interest (J| DEFINE, square, EQUALS, DUP, MUL, SEMICOLON, cube, EQUALS, DUP, DUP, MUL, END ).Evaluate();

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

7 Пример реализации функции сопостав ления строки с образцом Ниже приводится пример полной программы на языке С++, определяющей в синтаксисе языка Joy функцию match сопоставления строки с образцом. Знак * в образце обознача ет произвольную подстроку, а знак ? один произвольный символ. В примере используется комбинатор cond, семантика которого унаследована из языка Lisp. Функция ожидает на сте ке 2 параметра: первый шаблон, второй сопоставляемая строка. После определения функции, средставами языка C за прашиваются у пользователя шаблон и строка, которую нужно сопоставить. И, наконец, определенная функция применяется к введенным параметрам. Функция возвращает 1, если строка Сборник статей молодых ученых факультета ВМК МГУ, № 10, Библиотечная реализация языка Joy соответствует шаблону, и 0 в случае неуспеха. Чтобы не загро мождать пример, ответ функции не анализируется, а просто выводится на экран (функцией print языка Joy).

#include stdio.h #include "joycont.hpp" #include "joynames.hpp" #include "joyqprogram.hpp" #include "jfunction.hpp" int main() { JoyQProgramConstructor J;

JoyContinuation stacks;

JFunctionReference PRINT(new JFunctionPrint());

JFunctionReference EQ(new JFunctionEq());

JFunctionReference IFTE(new JFunctionIfTE());

JFunctionReference FIRST(new JFunctionFirst());

JFunctionReference REST(new JFunctionRest());

JFunctionReference DIP(new JFunctionDip());

JFunctionReference COND(new JFunctionCond());

JFunctionReference DROP(new JFunctionDrop());

JToken DEFINE("DEFINE");

JToken EQUALS("EQUALS");

JToken SEMICOLON("SEMICOLON");

JToken END("END");

JLabel match("match");

JLabel starmatch("starmatch");

JLabel drop2("drop2");

(J| DEFINE, match, EQUALS, (J| (J| (J| FIRST, ’\0’, EQ), DROP, FIRST, ’\0’, EQ), (J| (J| FIRST, ’?’, EQ), (J| DROP, "", EQ), Сборник статей молодых ученых факультета ВМК МГУ, № 10, 146 А. В. Мошкина (J| drop2, false), (J| (J| REST), DIP, REST, match), IFTE), (J| (J| FIRST, ’*’, EQ), REST, starmatch), (J| (J| true), (J| (J| FIRST), DIP, FIRST, EQ), (J| (J| REST), DIP, REST, match), (J| drop2, 0), IFTE) ), COND, SEMICOLON, starmatch, EQUALS, (J| (J| (J| match), drop2, true), (J| (J| DROP, "", EQ), drop2, false), (J| (J| true), (J| REST), DIP, starmatch) ), COND, SEMICOLON, drop2, EQUALS, DROP, DROP, END ).Evaluate(stacks);

char pattern[200];

char string[200];

printf("%s", "Please, enter the pattern: ");

scanf("%s", pattern);

printf("%s", "Please, enter the string: ");

scanf("%s", string);

(J| string, pattern, match, PRINT).Evaluate(stacks);

return 0;

} Сборник статей молодых ученых факультета ВМК МГУ, № 10, Библиотечная реализация языка Joy 8 Заключение Целью данной работы являлась реализация библиотечного расширения языка C++ для работы со стековым функциональ ным языком Joy. В ходе работы применялся метод непосред ственной интеграции и использовалась библиотека InteLib. В результате работы была реализована вычислительная модель языка Joy и библиотека стандартных функций языка.

Несмотря на то, что код, моделируемый конструкциями С++, похож по синтаксису на оригинальный код Joy, он оста ется достаточно громоздким и требует написания дополни тельных инициализирующих конструкций. Поэтому реализа ция была дополнена транслятором, переводящим модули в ори гинальном синтаксисе Joy, в модули на языке С++, которые можно непосредственно включить в проект. Это позволяет не только избавить разработчика от необходимости написания ко да Joy в синтаксисе С++, но и дает возможность без больших трудозатрат преобразовать уже существующий код на языке Joy для использования в проекте на С++.

Список литературы [1] А. В. Столяров. Об одном подходе к построению универ сальных языков программирования. Сборник статей моло дых учёных факультета ВМиК МГУ, N 4. М.: Издатель ский отдел факультета ВМиК МГУ, 2007, стр. 135–146.

[2] И. Г. Головин, А. В. Столяров. Объектно ориентированный подход к мультипарадигмальному программированию. Вестник Московского Университета, сер. 15 вычисл. матем. и киберн., 2002, № 1, стр. 46–50.

[3] А. В. Столяров. Импорт вычислительной модели языка Scheme в объектно–ориенторованное окружение. Сборник статей молодых учёных факультета ВМиК МГУ, № 5. М.:

Сборник статей молодых ученых факультета ВМК МГУ, № 10, 148 А. В. Мошкина Издательский отдел факультета ВМиК МГУ, 2008, стр.

119–130.

[4] И. Е. Бронштейн, А. В. Столяров. Библиотечная поддерж ка вычислительной модели языка Рефал. Сборник статей молодых учёных факультета ВМиК МГУ, № 6. М.: Изда тельский отдел факультета ВМиК МГУ, 2009, стр. 36–46.

[5] И. В. Струков. Библиотечная реализация Пролог– решателя. Сборник статей молодых ученых факультета ВМК МГУ, № 9, 2012. М.: Издательский отдел факультета ВМиК МГУ, 2012, стр. 235–250.

[6] Leo Brodie. Thinking FORTH. A Language and Philosophy for Solving Problems. Englewood Clis, N. J., Prentice–Hall, Inc., 1984.

[7] Manfred von Thun. Joy: Forth’s Functional Cousin.

Proceedings of the 17th EuroForth Conference, 9 October 2001.

[8] Glenn C. Reid. Thinking in PostScript. Addison–Wesley, 1990.

[9] Christopher Diggins. Cat: A Functional Stack-Based Little Language. Doctor Dobbs Journal, April 15th, 2008.

[10] Slava Pestov, Daniel Ehrenberg, Joe Gro. Factor: a dynamic stack-based programming language. Dynamic Languages Symposium 2010.

[11] Э. Хювёнен, И. Сеппянен. Мир Лиспа. В 2-х т. Мир, 1990.

[12] Официальный сайт проекта InteLib. http://www.intelib.org [13] Бьерн Страуструп. Дизайн и эволюция языка C++. ДМК Пресс, Питер, 2006.

Сборник статей молодых ученых факультета ВМК МГУ, № 10, Библиотечная реализация языка Joy [14] R. Kesley, W. Clinger and J. Rees. Revised5 report on Algorithmic Language Scheme. ACM SIGPLAN Notices, 1998.

[15] В. Ф. Турчин Алгоритмический язык рекурсивных функ ций (РЕФАЛ). М.: ИПМ АН СССР, 1968.

Сборник статей молодых ученых факультета ВМК МГУ, № 10, 150 Сборник статей молодых ученых факультета ВМК МГУ, № 10, УДК 519. Задача оптимального управления в модели распространения вируса гриппа A(H1N1) © 2013 г. С. М. Орлов sergey.orlov@cs.msu.su Кафедра оптимального управления 1 Введение В работе рассматривается нелинейная задача оптимального управления x1 x x1 = a (u + b1 )x1 c, x1 (0) = x10, x1 + x2 + x x1 x x2 = b2 x2 + c x2 (0) = x20,, x1 + x2 + x (1) x3 = ux1 + b3 x2 b1 x3, x3 (0) = x30, T J(u) = Qx2 + u2 dt min, u U.

2 u(·) Здесь x = (x1, x2, x3 ) трёхмерная фазовая переменная, u скалярное управление, подчинённое геометрическому ограни чению u U [0, 0.9], T 0 заданная длительность про цесса управления, параметры a, b1, b2, b3, c, x10, x20, x30, Q заданные положительные константы. Переменная x1 характе ризует количество людей, склонных к заболеванию, x2 ко личество инфицированных, x3 количество людей, имеющих иммунитет к заболеванию. Управление u характеризует интен сивность вакцинации. Функционал взвешенная характери стика количества заболевших и затрат на лечение подлежит минимизации. Модель (1) предложена в статьях [2,3]. В ориги нальных работах используются другие обозначения для пере менных задачи.

Модель распространения вируса гриппа A(H1N1) Для решения задачи привлекается принцип максимума Понтрягина [1]. Проведена серия численных экспериментов по нахождению решения краевой задачи ПМП. Вычисления про водились в среде Maple. При выполнении численных экспери ментов использовались идеи широко известного метода про должения по параметру (см., например, [4]).

Кроме того, изучен ряд теоретических вопросов, возника ющих при анализе задачи (1) сложной нелинейной задачи оптимального управления.

2 Модель В статье [2] исследуется следующая задача оптимального управления dS = µS S I uS, S(0) = S0 0, dt N dI I dt = S N (µ + d + r)I, I(0) = I0 0, dR (2) dt = rI µR + uS, R(0) = R0 0, tf A J(u) = I(t) + u2 (t) dt min, u U.

2 u(·) Здесь N = S(0) + I(0) + R(0) общая численность населения.

Область управления есть отрезок U [0, 0.9]. Обозначения для начальных условий модели (2) изменены по сравнению со ста тьями [2, 3] для отличия начального значения R(0) от базовой скорости репродукции R0, обсуждаемой ниже.

Похожие модели описаны в книге [5] (см. [5, с. 70, с. 134, с. 638]).

Сборник статей молодых ученых факультета ВМК МГУ, № 10, 152 С. М. Орлов 2.1 Интерпретация модели (2) с точки зрения ин фекционных процессов Эта часть работы написана с использованием русскоязыч ной терминологии из монографии [5] и данных из статей [2, 3].

В апреле 2009 года в Мексике и США началась пандемия нового штамма вируса гриппа A(H1N1), которая затем распро странилась почти по всему миру. В СМИ этот вирус был попу лярен под названием свиного гриппа. В связи с этим собы тием для моделирования распространения этого заболевания в статье [3] была предложена модель (2) (без управляющих сла гаемых ±uS в системе ОДУ и функционала).

Для борьбы с эпидемией правительства нескольких госу дарств решили провести массовую вакцинацию населения. В этом контексте в работе [2] была предложена и исследована задача оптимального управления (2).

Рассмотрим модель (2). Здесь общая популяция хозяев N разделена на три категории:

количество восприимчивых индивидов (susceptible, S(t) см. [2, 3]) в момент времени t;

количество инфицированных индивидов (infective) в I(t) момент времени t;

количество иммунных индивидов (removed) в момент R(t) времени t.

В дифференциальных уравнениях модели (2) имеются следу ющие эпидемиологические параметры: скорость притока людей в популяцию, µ удельная скорость естественной гибе ли и гибели от причин, не связанных с болезнью, d удельная скорость гибели от болезни, r удельная скорость выздоров ления, параметр трансмиссивности, объединяющий боль шое количество эпидемиологических, природных и социальных факторов, влияющих на скорость передачи инфекции. Слага I емое S N показывает, что новые случаи инфекции возникают пропорционально численности восприимчивых и доли инфици рованных хозяев.

Сборник статей молодых ученых факультета ВМК МГУ, № 10, Модель распространения вируса гриппа A(H1N1) Для модели (2) имеет важное значение параметр R0 = = µ+d+r базовая скорость репродукции инфекции. Это сред нее число вторичных инфекций, вызываемых внедрением од ного инфицированного индивида в популяцию, состоящую ис ключительно из восприимчивых хозяев. Ясно, что для внед рения и закрепления в популяции хозяев паразитические ви ды должны иметь базовую скорость репродукции инфекции R0 1.

Параметр u характеризует скорость вакцинации, описыва ющую переход восприимчивых индивидов S(t) в класс иммун ных R(t) (долю восприимчивых индивидов, подвергающихся вакцинации, в единицу времени). Задача оптимального управ ления состоит в минимизации количества инфицированных ин дивидов и стоимости вакцинации. Математически это форму лируется в виде интеграла в модели (2), подлежащего миними зации, а параметр A характеризует соотношение важности между количеством инфицированных и стоимостью вакцина ции (чем больше параметр A, тем важнее минимизировать рас ходы на вакцинацию по сравнению с количеством инфициро ванных).

2.2 Модификация модели (2) и новые обозначения В данной работе исследована модификация модели (2), ко гда общая численность населения N = N (t) меняется со вре менем t. Для удобства были введены новые обозначения для переменных и параметров задачи x1 = S, x10 = S0, b1 = µ, a =, x2 = I, x20 = I0, b2 = µ + d + r, c =, x3 = R, x30 = R0, b3 = r, Q = 1/A, T = tf.


Сборник статей молодых ученых факультета ВМК МГУ, № 10, 154 С. М. Орлов 3 Исследование системы дифференциаль ных уравнений модели (1) 3.1 Положительность фазовых переменных Рассмотрим систему дифференциальных уравнений для за дачи (1). Так как xi0 0, i = 1, 2, 3, из этого следует, что су ществует такое время, что на всём отрезке [0, ] решение x(t) лежит в R+ (все компоненты xi (t) 0). Справедливы следую щие оценки, верные на отрезке [0, ]:

x1 (t) (u + b1 + c)x1 (t), x2 (t) b2 x2 (t), x3 (t) b1 x3 (t).

Тогда в силу теоремы сравнения Чаплыгина каждая компонен та xi (t), i = 1, 2, 3, ограничена снизу положительной экспонен той. Легко проверить, что это ограничение верно для любого отрезка времени [0, s], на котором решение существует (можно доказать от противного).

3.2 Продолжимость решения Обозначим X = x1 + x2 + x3. Сложив три дифференциаль ных уравнения задачи (1), получим следующее дифференци альное уравнение:

X = a b1 X (b2 b3 b1 )x2.

Заметим, что b2 b3 b1 = d 0, и тогда справедлива оценка X a b1 X. Тогда для любого t 0 верно, что a + X(0)eb1 t.

0 X(t) b Тогда в силу положительности переменных x1 (t), x2 (t) и x3 (t) оказывается, что каждая из этих компонент ограничена для Сборник статей молодых ученых факультета ВМК МГУ, № 10, Модель распространения вируса гриппа A(H1N1) любого t 0:

a 0 xi (t) + X(0), i = 1, 2, 3.

b Тогда любое решение системы дифференциальных уравнений задачи (1) продолжимо на луч [0, +).

3.3 Устойчивость положений равновесия Рассмотрим систему дифференциальных уравнений зада чи (1) при u = 0. Запишем характеристическое уравнение для матрицы Якоби системы дифференциальных уравнений зада чи (1):

b1 c x2 (x2 +x3 ) c x1 (x1 +x3 ) c xXx X2 X2 x2 (x2 +x3 ) x1 (x1 +x3 ) x1 x2 = 0 (3) b2 + c X 2 c X c X 0 b b Легко показать, что у системы уравнений задачи (1) все ( ba, 0, 0) и второе го два положения равновесия: первое, x, x ), определяемое по формулам (x1 2 x = X, R b1 (R0 1) x = X, R0 (b1 + b3 ) b3 (R0 1) x = X, R0 (b1 + b3 ) aR0 (b1 + b3 ) X =, b1 ((R0 1)b2 + b1 + b3 ) c R0 =.

b Параметр R0 здесь есть базовая скорость репродукции инфек ции. Можно показать, что когда R0 1, то второго положения Сборник статей молодых ученых факультета ВМК МГУ, № 10, 156 С. М. Орлов равновесия нет (какие-то компоненты становятся отрицатель ными), а когда R0 = 1, то второе положение равновесия совпа дает с первым.

Рассмотрим характеристическое уравнение (3) для первого положения равновесия ( + b1 )2 (c b2 ) = 0, или ( + b1 )2 (b2 (R0 1) ) = 0.

Собственные значения 1,2 = b1 0, а 3 = b2 (R0 1). Зна чит, при R0 1 положение равновесия ( ba, 0, 0) асимптотиче ски устойчиво. Это ситуация, когда инфекции нет.

Теперь запишем уравнение (3) для (x, x, x ). Будем пред полагать, что здесь R0 1. После упрощений оно принимает вид ( + b1 )(2 + 1 + 2 ) = 0, (4) где b1 c(R0 1) 1 = b1 +, R0 (b1 + b3 ) b1 c(R0 1) 2 = (b2 (R0 1) + b1 b3 ).

R0 (b1 + b3 ) Видно, что 1 0, а при выполнении соотношения b2 (R 1) + b1 b3 0 верно, что 2 0. Тогда все решения урав нения (4) имеют отрицательные вещественные части, а значит, положение равновесия (x, x, x ) асимптотически устойчиво.

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

Сборник статей молодых ученых факультета ВМК МГУ, № 10, Модель распространения вируса гриппа A(H1N1) 4 Краевая задача принципа максимума Составим функцию Гамильтона–Понтрягина для задачи (1) 1 x1 x K = Qx2 + u2 + 1 a (u + b1 )x1 c + 2 x1 + x2 + x x1 x + 3 [ux1 + b3 x2 b1 x3 ], (5) + 2 b2 x2 + c x1 + x2 + x полагая в ней 0 = 1. Принимая во внимание соотношения Ku = u + (3 1 )x1, Ku = 1 0, из условия максимума K max u[0,0.9] находим максимизатор (6) u (t, x, ) = argmax K(t, x,, u) u[0,0.9] функции Гамильтона–Понтрягина (5). Максимизатор (6) опре деляется однозначно и может быть описан конструктивно. Вве дём обозначение = (3 1 )x для корня уравнения Ku = 0 относительно аргумента u.

Несложный анализ даёт для максимизатора (6) следующее вы ражение 0, 0, u () =, 0 0.9, 0.9, 0. или (7) u () = max(0, min(0.9, )).

Сборник статей молодых ученых факультета ВМК МГУ, № 10, 158 С. М. Орлов Запишем сопряжённые уравнения (1 2 )x2 (x2 + x3 ) 1 = b1 1 u(3 1 ) + c, (x1 + x2 + x3 ) (1 2 )x1 (x1 + x3 ) 2 = Q + b2 2 b3 3 + c, (x1 + x2 + x3 ) (1 2 )x1 x 3 = b1 3 c (x1 + x2 + x3 ) и условия трансверсальности 1 (T ) = 0, 2 (T ) = 0, 3 (T ) = 0.

Таким образом, краевая задача принципа максимума принима ет вид x1 x x1 = a (u () + b1 )x1 c, x1 + x2 + x x1 x x2 = b2 x2 + c, x1 + x2 + x x3 = u ()x1 + b3 x2 b1 x3, 1 = b1 1 u ()(3 1 ) + c (1 2 )x2 (x2 + x3 ), (x1 + x2 + x3 )2 (8) 2 = Q + b2 2 b3 3 + c (1 2 )x1 (x1 + x3 ), (x1 + x2 + x3 ) 3 = b1 3 c (1 2 )x1 x2, (x1 + x2 + x3 ) x (0) = x, x (0) = x, x (0) = x, 10 2 20 3 1 (T ) = 0, 2 (T ) = 0, 3 (T ) = 0.

5 Численные эксперименты 5.1 Значения параметров модели Данные о значениях параметров для численных экспери ментов были взяты из работ [2, 3]. Эти данные основаны на Сборник статей молодых ученых факультета ВМК МГУ, № 10, Модель распространения вируса гриппа A(H1N1) анализе ситуации в Морокко (2009 г.). Одна временная едини ца обозначает один день, удельная скорость естественной ги бели была взята равной µ = 3.9139 105 в день (она посчи тана как обратная к величине ожидаемой продолжительности жизни от рождения в Морокко, которая составляет пример но 70 лет). Параметр взят равным = µN (0) = 1174.17.

Удельная скорость гибели от болезни положена равной d = = 0.63 %, так как по данным от 10 июня 2009 года в Морок ко было зафиксировано всего 2208 случаев заболевания и смертей от него. Количество дней на лечение взято равным 5-ти, поэтому удельная скорость выздоровления взята равной r = 5 = 0.2 индивидуума в день. Базовая скорость репродук ции инфекции R0 = 1.5, поэтому можно посчитать параметр = R0 (µ + r + d) 0.3095. Начальные значения переменных были взяты равными S(0) = 30 106, I(0) = 28, R(0) = 30. В статье [2] значение параметра A не приведено, и из некоторых эмпирических соображений он был взят равным A = 50 000.

Период вакцинации был взят равным tf = 250 дней.

5.2 Численное решение краевой задачи Краевая задача (8) исследовалась численно в среде Maple с помощью стандартной функции dsolve и встроенной в неё реализации метода продолжения по параметру. Было установ лено, что краевая задача (8) легко решается в среде Maple, когда краевые условия на левой границе отрезка времени [0, T ] выглядят следующим образом x1 (0) = x10, x2 (0) = 0, x3 (0) = x10.

Решение задачи (8) с такими краевыми условиями было взято в качестве начального приближения в методе продолжения по параметру. Для облегчения счёта были сделаны следующие за мены переменных: x1 = x1 /x10, x2 = 2x2 /x20, x3 = x3 /x30. Для контроля точности решения краевой задачи (8) проводилось исследование соответствующей задачи Коши.

Сборник статей молодых ученых факультета ВМК МГУ, № 10, 160 С. М. Орлов Ниже на рисунках 1–4 показаны графики экстремального управления, зависимости фазовых переменных от времени и график функции t Qx2 (t) + u2 (t) x0 (t) = dt.

Эта функция была посчитана как решение задачи Коши x0 = Qx2 + u2, (9) x0 (0) = с использованием найденного численного решения краевой за дачи (8). Она показывает характер изменения значения функ ционала J(u) при увеличении отрезка времени, по которому бе рётся интеграл, и верно равенство x0 (T ) = J(u ). Кроме того, x0 (t) u (t) Jextr t t Рис. 2. Управление u (t) Рис. 1. Функция x0 (t) в среде Maple вычислено значение функционала J задачи (1) для найденного решения краевой задачи принципа максимума методом Clenshaw–Curtis quadrature с точностью = 104 :

(10) Jextr = 0.0372,  Сборник статей молодых ученых факультета ВМК МГУ, № 10, Модель распространения вируса гриппа A(H1N1) x2 (t) x3 (t) x1 (t) t t Рис. 3. Траектории x1 (t), x3 (t) Рис. 4. Траектория x2 (t) причём имеет место следующее соотношение T T Qx2 dt 3/4, u dt/ 0 то есть значителен вклад обоих слагаемых в значение функци онала. Отметим, что x0 (T ), посчитанное как решение задачи Коши (9) в конечный момент времени T, совпадает до 4-го зна ка после запятой со значением функционала (10), посчитанным численным методом интегрирования.

Были посчитаны значения функционала задачи (1) для раз личных управлений:

u = 0, J(u) = 1721.9613, u = 0.9, J(u) = 101.2536, 0.9, t u=, J(u) = 0.4125.

0, t Видно, что значение функционала на экстремальном управле нии меньше, чем на других посчитанных управлениях. После Сборник статей молодых ученых факультета ВМК МГУ, № 10, 162 С. М. Орлов получения значений функционала J при различных управле ниях было решено провести следующее исследование.

5.3 Численное исследование значений функцио нала в двух классах управлений Был рассмотрен класс управляющих функций с одной точ кой переключения:

0.9, t, (11) u(t;

) = 0, t, и для него введена функция I1 ( ) = J(u(t;

)) значение функционала J на управлении u(t;

). С помощью стандарт ной функции Minimize из пакета Optimization, эта функция была численно минимизирована на отрезке времени [0, T ]. Най денное время переключения равно 1 = 0. и соответствующее ему значение функционала равно I1 = I1 (1 ) = 0.2545 Jextr.

На рисунках 5–7 показаны графики функции I1 ( ) на отрезках времени [0, 0.4], [0.4, 1], [1, 250].

После этого был рассмотрен более общий класс управлений с одной точкой переключения u0, t, (12) u(t;

, u0 ) = 0, t, и для него введена функция I2 (, u0 ) = J(u(t;

, u0 )) значение функционала J на управлении u(t;


, u0 ). Для подсчёта мини мума этой функции на квадрате [0, T ] [0, 0.9] минимизиро валась функция I3 (u0 ) = min [0,T ] I2 (, u0 ) на отрезке [0, 0.9].

Сборник статей молодых ученых факультета ВМК МГУ, № 10, Модель распространения вируса гриппа A(H1N1) I1 ( ) I1 ( ) Рис. 5. График функции I1 ( ) на Рис. 6. График функции I1 ( ) на отрезке [0, 0.4] отрезке [0.4, 1] I1 ( ) Рис. 7. График функции I1 ( ) на отрезке [1, 250] Функция I3 (u0 ) строилась в Maple с использованием функции Minimize и её минимизация также проходила при помощи этой стандартной функции методом Branch-and-Bound. Найден ные значения параметров равны u = 0.0469, 2 = 16.4694, Сборник статей молодых ученых факультета ВМК МГУ, № 10, 164 С. М. Орлов а значение функционала равно I2 = I2 (2, u ) = 0.041 Jextr.

В среде Maple были построены графики функции I2 ( ) = ) на отрезках [0, 16], [16, 17], [17, 250], и оказалось, что = I2 (, u качественное поведение этих графиков совпадает с поведением графиков, представленных на рисунках 5–7. Аналогично были исследованы графики функции I2 (u0 ) = I2 (2, u0 ) на отрезках [0, 0.04], [0.04, 0.05], [0.05, 0.9], и их поведение тоже соответству ет поведению графиков на рисунках 5–7.

Ниже на рисунках 8–10 показаны графики функции I3 (u0 ) на отрезках [0, 0.03], [0.03, 0.06] и [0.06, 0.9], а на рисунке показана поверхность I2 (, u0 ) в окрестности точки (2, u ).

I3 (u0 ) I3 (u0 ) u0 u Рис. 8. График функции I3 (u0 ) Рис. 9. График функции I3 (u0 ) на отрезке [0.03, 0.06] на отрезке [0, 0.03] Из проведённого исследования можно заключить, что экс тремальное управление, полученное из решения краевой зада чи принципа максимума (8), лучше, чем полученные оптималь ные управления в классах (11), (12).

Сборник статей молодых ученых факультета ВМК МГУ, № 10, Модель распространения вируса гриппа A(H1N1) I3 (u0 ) u Рис. 10. График функции I3 (u0 ) на отрезке [0.06, 0.9] Рис. 11. График функции I2 (, u0 ) в окрестности точки (2, u ) 6 Заключение В статье исследована управляемая модель (1) на основе принципа максимума Понтрягина [1], найден максимизатор (7) функции Гамильтона–Понтрягина (5), составлена краевая за дача принципа максимума (8), проведены численные экспери Сборник статей молодых ученых факультета ВМК МГУ, № 10, 166 С. М. Орлов менты по нахождению экстремального управления и числен ный поиск оптимальных управлений в классах управляющих функций (11) и (12). Установлено, что численно полученное экстремальное управление доставляет меньшее значение функ ционалу J(u), чем соответствующие оптимальные управления в классах (11), (12).

Кроме того, доказана положительность фазовых перемен ных в модели (1), продолжимость решения задачи Коши при фиксированной управляющей функции и исследованы на устойчивость положения равновесия системы (1) при u = 0.

Список литературы [1] Понтрягин Л. С., Болтянский В. Г., Гамкрелидзе Р. В., Мищенко Е. Ф. Математическая теория оптимальных про цессов. М.: Физматлит, 1961. 391 с.

[2] El hia M., Balatif O., Bouyaghroumni J., Labriji E., Rachik M. Optimal Control Applied to the Spread of Inuenza A(H1N1) // Applied Mathematical Sciences. 2012.

Vol. 6, no. 82. P. 4057–4065.

[3] Hattaf K., Yous N. Mathematical Model of the Inuenza A(H1N1) Infection // Advanced Studies in Biology. 2009.

Vol. 1, no. 8. P. 383–390.

[4] Аввакумов С. Н., Киселёв Ю. Н. Некоторые алгоритмы оп тимального управления // Управление, устойчивость и об ратные задачи динамики, Сборник научных трудов, Тр.

ИММ УрО РАН, 12. 2006. № 2. С. 3 17.

[5] Андерсон Р., Мэй Р. Инфекционные болезни человека.

Динамика и контроль: Пер. с англ. А. А. Романюхи и С. Г. Руднева под редакцией Г. И. Марчука. М.: Мир, Научный мир, 2004. 784 c. ISBN 5–03–003552– ( Мир ). ISBN 5–589176–285–4 ( Научный мир ).

Сборник статей молодых ученых факультета ВМК МГУ, № 10, Сборник статей молодых ученых факультета ВМК МГУ, № 10, 2013 УДК 004. Численное решение задачи о свободной конвекции в квадратной области в двумерной постановке конечно-разностным методом © 2013 г. В. О. Пиманов, И. С. Калачинская pimanov-vladimir@yandex.ru, kalach@cs.msu.su Кафедра вычислительных методов Введение Работа посвящена конечно-разностному методу для реше ния системы уравнений Навье-Стокса в приближении Буси неска, но не его описанию, а тестированию. В работе [2] пред ставлено несколько численных методов, в том числе метод ко нечных элементов и метод сингулярной свертки, решающих ту же задачу. Для их сравнения в статье [2] была выбрана зада ча о свободной конвекции в ее двумерной постановке в квад ратной каверне. На ней же был опробован и представленный численный метод, а результаты сведены в общие таблицы для сравнения. Качественно решение также освещено в графиках и plot-построениях.

В основе представленного метода лежит метод решения си стемы уравнений Навье-Стокса для несжимаемой жидкости без учета тепловых процессов [1]. Не представляет сложности реализовать этот метод и для трехмерной задачи.

Задача о свободной конвекции в прямоугольной области бы ла предложена еще в 80-х годах прошлого века как тестовая за дача для проверки и отладки программного кода и численных методов, решающих широкий класс задач, связанных с тепло выми процессами. Например, теплоизоляция реактора, охла ждение процессора, вентиляция помещений, пожарная безопас 168 В. О. Пиманов и И. С. Калачинская ность, получение солнечной энергии, отвод и распределение тепла или рост кристаллов в жидкости.

Система уравнений Навье-Стокса в при ближении Буссинеска Уравнения тепловой конвекции в приближении Буссинеска Обербека наиболее популярная модель для описания конвек ции в жидкостях и газах.

Модель включает уравнение Навье-Стокса, уравнение теп лопроводности и уравнение несжимаемости. Основная идея приближения состоит в особенности учёта зависимости плотно сти от температуры. Именно, в системе уравнений конвекции данная зависимость учитывается только при массовых силах:

v + (v )v = p + v + (T )g, t T + v · T = T, t div v = 0.

Здесь v скорость течения, T отклонение температуры, p давление, динамическая вязкость, коэффициент тем пературопроводности, g ускорение свободного падения.

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

(T ) = 0 (1 T ).

Здесь коэффициент объёмного расширения жидкости, плотность жидкости при некоторой равновесной температуре, T отклонение температуры от равновесия. Поскольку и от клонение температуры обычно относительно невелико, то ли нейное приближение обладает приемлемой точностью в боль шинстве исследуемых задач.

Сборник статей молодых ученых факультета ВМК МГУ, № 10, Задача о свободной конвекции Подстановка линейной зависимости плотности и перенор мировка давления позволяют исключить слагаемое 0 g. Окон чательно задача конвекции несжимаемой жидкости в прибли жении Буссинеска принимает следующий вид:

v + (v · )v = p + v T g, t T + v · T = T, t div v = 0.

Здесь кинематическая вязкость.

Приведённая задача конвекции в различных постановках неоднократно исследовалась. Наиболее широко известна зада ча Рэлея-Бенара о конвекции в плоском слое жидкости. При определённых условиях возможно точное решение задачи, на пример, для ламинарной конвекции в вертикальном слое при подогреве сбоку (иногда встречается под названием задача Гершуни ). Материал из Википедии.

Постановка задачи Форма каверны. Каверна имеет вытянутую форму с квад ратным поперечным сечением, так что ее длина в продольном направлении несравнимо больше и ширины, и высоты. Таким образом, постановка задачи двумерная. В приведенных пере менных в поперечном сечении квадрат размером 1.01.0. Пред полагаем, что решение также двумерное. От времени форма ка верны не зависит и стенки неподвижны. Изображение каверны на рисунке 1.

Тепловой режим. Верхняя и нижняя стенки теплоизо лированные. В приведенных переменных температура на левой стенке 1.0, на правой 0.0 постоянная.

Физические параметры. Число Пранделя (Pr) взято рав ным 0.71, что соответствует конвекции в воздухе. Число Ре Сборник статей молодых ученых факультета ВМК МГУ, № 10, 170 В. О. Пиманов и И. С. Калачинская Теплоизоляция, dT= Горячяя Холодная g стенка, стенка, T=1 T= Теплоизоляция, dT= Рис. 1. Постановка задачи лея (Ra) варьируется в диапазоне 103 Ra 108.

Интересующее решение. В работе исследуется стаци онарное решение поставленной задачи в диапазоне парамет ров 103 Ra 108. Результаты сопоставляются с результа тами известных тестов [2], и при Ra = 108 также реализуется стационарный режим.

Система уравнений Ставится внутренняя задача в квадратной области разме ром 1.01.0. Для моделирования процессов тепломассоперено са использована система уравнений Навье-Стокса в приближе нии Буссинеска. В безразмерных переменных уравнение имеет вид:

v = 0, v = (v · )v + Pr 2 v p + f, t T = (v · )T + 2 T.

t Здесь f = (0, 0, Ra Pr T ) подъемная сила, число Пранделя Pr = 0.71, вектор v поле скорости, скаляр p поле давление, Сборник статей молодых ученых факультета ВМК МГУ, № 10, Задача о свободной конвекции скаляр T поле температуры, число Релея Ra варьируемый параметр.

Граничные условия:

на границе, u= на левой стенке, T = на правой стенке, T = T на верхней и нижней стенках.

= n Численный метод Для решения системы уравнений предлагается использо вать метод, аналогичный [1], разработанный для решения си стемы уравнений Навье-Стокса без учета тепловых процессов.

Чтобы коротко охарактеризовать его, можно сказать, что это конечно-разностный метод решения уравнения Навье-Стокса для вязкой несжимаемой жидкости с полунеявным методом Рунге-Кутта третьего порядка для интегрирования по време ни. Метод построен на разнесенной неравномерной сетке.

Модификация метода [1] с учетом тепловых про цессов Температура, как и давление, определяется в центрах яче ек. Приближение Буссинеска представляет себой уравнение Навье-Стокса с учетом температуры только как внешней си лы. Поэтому, модифицируя метод [1], необходимо реализовать лишь интегрирование температуры по времени, так как в ме тоде [1] учет внешних сил реализован.

Уравнение эволюции тепла в векторной и скалярной фор Сборник статей молодых ученых факультета ВМК МГУ, № 10, 172 В. О. Пиманов и И. С. Калачинская мах:

T = (v · )T + 2 T, t T T T + 2 T.

= u v t x y Здесь T температура, v = (u, v) скорость.

Определим разностную аппроксимацию:

T x y = udx T vdy T + [dx dx + dy dy ]T.

t Здесь верхнее подчеркивание обозначает интерполяцию по со ответствующей переменной, d разностную производную по соответствующей переменной так, как это определено в ста тье [1].

Граничные условия:

x T = 1, при x= x T = 0, при x= dy T = 0, при y=0 или Этого достаточно, чтобы в тот момент, когда вычисляются скорости на новом слое по времени, вычислить еще и темпера туру на новом слое по времени, применив такой же полунеяв ный метод Рунге-Кутта третьего порядка, что и в статье [1].

Делая шаг по времени, мы используем при интегрировании температуры значение скорости с предыдущего слоя и при ин тегрировании скорости значение температуры с предыдуще го слоя. Однако есть возможность использовать при вычисле нии скорости, например, температуру, вычисленную для того же слоя по времени, что и сама скорость. Или наоборот.

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

Сборник статей молодых ученых факультета ВМК МГУ, № 10, Задача о свободной конвекции Тестовый расчет Для проверки правильности работы написанной програм мы и численного метода результаты расчета сравниваются с аналогичными результатами, представленными в работе [2]. Во всех расчетах коэффициент растяжения выбран таким обра зом, что отношения линейного размера центральной ячейки к крайней равно 10.

Исследование зависимости точности решение от сетки В таблице 1 приведено несколько параметров стационар ного течения при различных числах Релея на различных сет ках. Сетки квадратные, кол-во ячеек по одной и по второй оси равны друг другу. Приведенные величины: максимальное зна чение x-компоненты скорости umax, максимальное значение y компоненты скорости vmax, и изменение среднего значения чис ла Нуссельта вдоль горячей стенки N ucp.

Анализируя таблицу, можно заметить, что при числах Ре лея порядка 107, 108 на сетке 128128 погрешность еще ве лика и может достигать порядка одного-двух процентов. Есть смысл уменьшить расчетную сетку для достижения большей точности.

Эволюция во времени и выход на стационар Для визуализации была выбрана функция E(t), вычислен ная как максимальная разность в расчетной области между значением завихренности в настоящий момент t и в момент времени, отстающий на 10dt, где dt шаг интегрирования по времени. Если расчетная область, то:

Сборник статей молодых ученых факультета ВМК МГУ, № 10, 174 В. О. Пиманов и И. С. Калачинская Таблица 1. Зависимость решения от размера сетки N и числа Релея Ra = 103, 104, 105, 106, 107, 2107, 4107, 108. umax макси мальное значение x-компоненты, vmax максимально значение y компоненты, Nucp среднее значение числа Нуссельта вдоль горя чей стенки. По вертикали и по горизонтали число узлов совпадает.

N Ra 103 104 105 106 107 2107 4107 umax 64 3.64675 16.1859 43.7988 127.812 385.351 535.085 739.406 1125. 96 3.65131 16.1917 43.8117 128.426 386.984 535.950 741.288 1125. 128 3.65293 16.0214 43.8711 128.497 387.922 537.127 742.532 1125. 160 3.65319 - - - - - - 192 3.65319 - - - - - - vmax 64 3.70227 19.6867 68.7262 220.899 687.043 984.398 1378.23 2232. 96 3.70036 19.6864 68.6058 221.158 701.871 990.193 1392.87 2228. 128 3.70091 19.6788 68.6241 221.170 701.852 991.994 1406.55 2233. 160 3.70174 - - - - - - 192 3.70155 - - - - - - Nucp 64 1.117808 2.24512 4.52035 8.82954 15.7462 19.0123 22.6798 30. 96 1.117799 2.24494 4.51143 8.81161 16.5496 19.3194 22.7602 30. 128 1.117795 2.24491 4.55260 8.82423 16.5379 19.7897 22.8603 30. 160 1.117793 - - - - - - 192 1.117792 - - - - - - v u = x y E(t) = max((t) (t 10dt)) При достаточно малом шаге по времени метод моделирует эволюцию по времени. Но, так как в этой работе исследуется только стационарное решение, сходимость по шагу не исследо валась. Левый график на рисунке 2 типичный, наблюдался при всех расчетах при достаточно мелких сетках. От числа Ре лея зависит лишь время выхода на стационар, что отображено на таблице 2. Правый график на рисунке 2 можно наблюдать при больших числах Релея и при недостаточно мелких сетках, но при увеличении числа узлов в сетке картина меняется на уже описанную выше. Оба графика соответствуют расчету при Сборник статей молодых ученых факультета ВМК МГУ, № 10, Задача о свободной конвекции 3 10 Ra=1e+8 Ra=1e+ N=128 N= 1 10 10-1 10- -3 - 10 E(t) E(t) 10-5 10- -7 - 10 -9 - 10 0 0.1 0.2 0.3 0.4 0 0.1 0.2 0.3 0. t t Рис. 2. Зависимость Е(t) от времени. Слева сетка 128128, справа сетка 9696. Число Релея Ra = 108.

числе Релея равном 108. Первый на сетке в 128128 точек, вто рой на сетке 9696, недостаточно мелкой для верного описания процесса. В момент времени t=0 скорости равны 0. При чис лах Релея до 108 получено стационарное решение, как и в ра боте [2]. Расчетов при больших числах Релея не проводилось.

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

Таблица 2. Зависимость времени выхода на стационар tc от числа Релея Ra. Примерные значения.

103 104 105 106 107 2107 4107 Ra tc 0.95 0.7 0.55 0.4 0.25 0.21 0.16 0. С увеличением числа Релея время выхода на стационар уменьшается.

Анализ решения и графики Все решения получены на сетке 128128 точек. Во всех расчетах коэффициент растяжения выбран таким образом, что Сборник статей молодых ученых факультета ВМК МГУ, № 10, 176 В. О. Пиманов и И. С. Калачинская отношение линейного размера центральной ячейки к крайней равно 10. Считается, что такой сетки достаточно для верного описания процесса.

0. 10^ 10^ 0. T 10^6 10^ 0. 0. 0 0.2 0.4 0.6 0.8 X Рис. 3. Зависимость температуры T от x координаты вдоль гори зонтальной прямой, проходящей через центр каверны (y = 0.5), при различных числах Релея Ra = 103, 104, 105, 106.

1 0.8 0. 10^ 10^ 0.6 0. T T 10^ 10^ 0.4 0. 0.2 0. 0 0 0.2 0.4 0.6 0.8 1 0 0.05 0.1 0.15 0. X X Рис. 4. Зависимость температуры T от x координаты вдоль гори зонтальной прямой, проходящей через центр каверны (y = 0.5), при числах Релея Ra = 107, 2107, 4107, 108. Графики отличает различ ный масштаб по горизонтальной оси.

Сборник статей молодых ученых факультета ВМК МГУ, № 10, Задача о свободной конвекции 10^ 10^ 10^ V 10^ - - - - 0 0.2 0.4 0.6 0.8 X Рис. 5. Изменения вертикальной проекции скорости v вдоль гори зонтальной линии, проходящей через центр каверны (y = 0.5), при числах Релея Ra = 103, 104, 105, 106.

10^8 10^ 2000 1500 1000 500 V V 10^7 10^ -500 - -1000 - -1500 - -2000 - 0 0.2 0.4 0.6 0.8 1 0 0.05 0.1 0.15 0. X X Рис. 6. Изменение вертикальной проекции скорости v вдоль гори зонтальной линии, проходящей через центр каверны (y = 0.5), при числах Релея Ra = 107, 2107, 4107, 108. Графики отличает различ ный масштаб по горизонтальной оси.

В каждом из приведенных на рисунках 3,4,5,6,7,8,9 графи ков наблюдается симметрия решения относительно центра ка верны. Это легко объяснимо, хотя заранее нельзя быть уве ренным в таком результате. Начальные и граничные условия Сборник статей молодых ученых факультета ВМК МГУ, № 10, 178 В. О. Пиманов и И. С. Калачинская 10^ U U 0 10^ - - 10^ - 10^5 - -50 - 10^ 10^6 - 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 Y Y Рис. 7. Изменение горизонтальной проекции скорости u вдоль вер тикальной линии, проходящей через центр каверны (x = 0.5). На при числах Релея Ra = 103, 104, 105, 106, на пра левом рисунке вом Ra = 107, 2107, 4107, 108.

10^ 10^ Nu Nu 10^ 10^ 10^ 10^ 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 Y Y Рис. 8. Изменения числа Нуссельта Nu вдоль холодной стенки (x = 1.0). На левом рисунке при числах Релея Ra = 103, 104, 105, 106, на правом Ra = 107, 2107, 4107, 108.

имеют ту же центральную симметрию. Но, вероятно, начиная с некоторого числа Релея ситуация изменится, и также, вероят но, при этом же числе Релея перестанет существовать и стаци онарное решение.

Сборник статей молодых ученых факультета ВМК МГУ, № 10, Задача о свободной конвекции 10^ 10^ Nu Nu 10^5 10^ 10^4 10^ 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 Y Y Рис. 9. Изменения числа Нуссельта Nu вдоль горячей стенки (x=1.0). На левом рисунке при числах Релея Ra = 103, 104, 105, 106, на правом Ra = 107, 2107, 4107, 108.

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

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

Plot-построения На рисунке 3 представлены изолинии функции тока, x проекции скорости (u), y-проекции скорости (v), температу ры (T) столбцы. Число Релея принимает значения Ra = 3,104,105,106,107,2107,4107,108 строки. Изолинии вы браны так, чтобы показать особенности решения, и в их по строении нет определенной структуры, которую можно было бы указать здесь.



Pages:     | 1 | 2 || 4 | 5 |   ...   | 6 |
 





 
© 2013 www.libed.ru - «Бесплатная библиотека научно-практических конференций»

Материалы этого сайта размещены для ознакомления, все права принадлежат их авторам.
Если Вы не согласны с тем, что Ваш материал размещён на этом сайте, пожалуйста, напишите нам, мы в течении 1-2 рабочих дней удалим его.