Рис. 2 .3. Пример объектного представления параметров и функций модели (Simulink).

2.2. Объектно-ориентированная методология моделирования

2.2.1. Построение моделей из элементов

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

Элемент должен представлять собой набор семантически связанных параметров, к которым привязаны определённые алгоритмы вычислительной математики, а точнее – те части алгоритмов, которые рассчитывают именно эти параметры. При формировании объектно-ориентированной библиотеки численных методов математику достаточно запрограммировать некоторое количество классов элементов, которые соответствуют либо типичным математическим объектам (к примеру, системе обыкновенных дифференциальных уравнений), либо типичным объектам рассматриваемой предметной области. Например, в гидромеханике или физиологии человека такими типичными объектами являются трубки (сосуды), резервуары, обменники и т.д. После того, как созданы классы элементов, из них, как из конструктора, эксперт по данной предметной области может строить разнообразные модели безо всякого программирования (см. раздел 2.3.2).

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

2.2.2. Создание типов элементов на основе наследования

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

Например, пусть исходно был запрограммирован класс проводника, то есть элемента, который характеризуется сопротивлением и (по)током через него и предназначен для расчёта электрической цепи или, что то же самое, гидродинамической системы жёстких трубок. Если требуется учесть растяжимость трубок, достаточно унаследовать свойства и методы класса проводника, добавив к ним свойство объёма и процедуру расчёта сопротивления трубки через этот объём. (см. рис. 2.4) Далее класс растяжимой трубки (то есть сосуда, если пользоваться терминологией предметной области – физиологии) можно использовать в качестве родительского класса для эластичных, пластичных и других сосудов, различающихся моделью растяжения.

А теперь предположим, что требуется решать не только эллиптическую задачу расчёта потоков и давлений в системе трубок, но и параболическую задачу переноса веществ или тепла в этой системе. Понятно, что методы расчёта конвективного, диффузионного и конвективно-диффузионного переноса должны существовать в специальных элементах-переносчиках независимо от элементов, связанных с моделированием механики сосудов. (см. рис. 2.4) Чтобы использовать эти методы для задачи переноса в сосудах, достаточно создать элемент, наследующийся от сосуда и от переносчика одновременно, причём никакого нового программного кода создавать не требуется. В объектно-ориентированных языках, не поддерживающих множественное наследование, то же самое достигается путём одиночного наследования от одного класса (например, от сосуда) и агрегации экземпляра другого класса (например, переносчика).

Рис. 2 .4. Пример иерархии наследования классов модели (в объектной нотации Г. Буча).

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

2.3. Пример использования: моделирование организма человека

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

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

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

2.3.1. Описание некоторых блоков модели и проблем их расчёта

Основой комплексной модели организма человека является модель кровеносной системы: во-первых, она служит посредником между всеми другими системами, осуществляя транспорт веществ в организме, во-вторых, параметры кровообращения являются основными показателями состояния человека. Поскольку кровеносное русло включает сосуды существенно различных размеров, их совместное моделирование требует комбинации различных математических подходов. Достаточно изотропную сеть мелких сосудов корректнее описывать моделью сплошной среды, то есть с помощью дифференциальных уравнений в частных производных. Это позволяет не вводить для таких сосудов лишних анатомических характеристик типа их размеров, сводя их свойства к одному параметру непрерывной задачи – проницаемости ткани по отношению к крови. С другой стороны, модель макроскопических сосудов сводится к системе с сосредоточенными параметрами, гидродинамическая часть которой рассчитывается по алгебраическим уравнениям (Кирхгоффа), а транспортная часть – по обыкновенным дифференциальным уравнениям. Проблема в том, что при попытке учесть в рамках такой дискретной модели большое число сосудистых поколений размерность сосредоточенной системы превысит размерность системы сеточных уравнений, аппроксимирующей непрерывную задачу в мелких сосудах.

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

Такой же вывод можно сделать при рассмотрении задачи совместного расчёта параметров сердца и сосудистой системы. Конечно, работа сердца описывается гораздо более сложными уравнениями (см. [11], а также табл. 3 приложения), чем закон пропорциональности между потоком и градиентом давления, характерный для сосудов, поэтому, казалось бы, нельзя рассчитывать их единообразно. Однако если решать задачу отдельно для обоих желудочков сердца при фиксированных свойствах сосудов (интегральном сопротивлении кругов кровообращения) и отдельно – для сосудов при фиксированных параметрах желудочков (артериальном давлении и потоке), то возникает, по крайней мере, два вложенных итерационных процесса (один из них – например, процесс решения задачи в сосудах – должен сходиться на каждой итерации другого процесса – например, процесса расчёта сердечных параметров). Очевидно, эффективность таких расчётов очень низка, и её можно повысить, если объединить итерации схем сердечных и сосудистых уравнений в один процесс. Это возможно только за счёт объектно-ориентированной трактовки элементов сосудистой системы и сердца и их наследования от общего суперкласса (см. раздел 3.4.2). Одновременно такой подход увеличивает гибкость модели: в процедурно-ориентированный алгоритм из двух вложенных процессов весьма затруднительно вставить расчёт дополнительных активных элементов кровеносной системы (например, желудочков сердца, мышечного, дыхательного или других «насосов»).

Ещё более очевидные проблемы, пригодные к решению с помощью объектно-ориентированного подхода, возникают при моделировании транспорта веществ в организме. Вещества могут попадать в организм через дыхательную или пищеварительную систему (в первом случае – за счет конвекции и диффузии, во втором – с участием мышц органов желудочно-кишечного тракта), проникать в кровь (диффузионным путём или за счёт всасывания), транспортироваться кровеносной системой, испытывая по дороге разнообразные превращения и абсорбируясь выделительной системой; наконец, вещества могут разными способами потребляться в органах и тканях. Уравнения для всех этих процессов довольно похожи, однако разные пространственные размерности соответствующих моделей (для пищеварения и выделения – точечные модели, для кровообращения и дыхания – ветвящиеся одномерные и квазитрёхмерные) не позволяют использовать для их расчёта одни и те же процедуры. Особенно сложно при процедурном подходе формализовать условия на границах раздела систем, число типов которых намного больше, чем число самих систем. Однако если представлять организм человека в виде совокупности связанных между собой элементов, то при разработке их методов не нужно заботиться о том, с какой элемент (какой системе принадлежит, является ли «граничным») связан с данным элемент. Кроме того, размерность системы фактически не влияет на объектные методы моделирования (подробнее это описано в разделе 3.2.3), а за счёт наследования алгоритмы расчётов сходных физических процессов в разных функциональных системах организма можно реализовывать только один раз.

Таким образом, даже поверхностный обзор некоторых блоков модели организма человека показывает, что для её реализации необходим объектно-ориентированный подход, – прежде всего, ввиду комплексного характера модели, неоднородности её элементов и сложности связей между ними.

2.3.2. Объектно-ориентированное представление структуры модели

В данном разделе в виде схем представлены результаты объектно-ориентированного анализа некоторых функциональных систем организма человека (более подробно они описаны в [12]). Их модели представляют собой совокупность элементов – объектов различных классов, – соединённых между собой направленными связями, смысл которых описан в разделе 2.2.1. Структура самих элементов подробно приведена в разделе 3.4.1 при рассмотрении объектно-ориентированной библиотеки численных методов, на которой основана данная модель. Классы элементов обозначены на схемах символическими изображениями, соответствие с которыми определяется таблицей 2.2. Каждая схема соответствует некоторой подсистеме организма. Если изображение какого-либо элемента уменьшено по сравнению с остальными и помещено в ромбическую рамку, то он относится к другой подсистеме, являясь связующим звеном между ней и подсистемой рассматриваемой схемы.

Следует заметить, что рисунки 2.5-2.7 не являются теоретическими схемами структуры модели. Они созданы с помощью специального графического редактора расчётной программы, то есть изображённые на рисунках объекты (элементы) создаются при инициализации модели и в ходе моделирования выполняют соответствующие их классам численные алгоритмы. При создании элементов используется, помимо приведённой структурной информации и декларативных связей, также количественная информация и функциональные связи. Функциональные связи (зависимости между параметрами элементов) также имеют прямое отношение к рассматриваемой здесь структуре модели и вводятся в неё с помощью графических редакторов, однако за недостатком места они не приводятся.

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

Рис. 2 .5. Общая структура модели кровеносной системы.

Рис. 2 .6. Структура модели малого круга кровообращения.

Рис. 2 .7. Структура модели большого круга кровообращения.

Таблица 2 .2

Условные обозначения некоторых классов элементов модели организма человека на рис. 2.5-2.7. Номера соответствуют таблицам 1-4 приложения.

2.4. Резюме

В данной главе рассмотрены проблемы математического и имитационного моделирования, которые, как показано, целесообразно решать на основе объектно-ориентированного представления структуры моделей. Путём анализа подходов, принятых в существующих объектных средствах моделирования, выбрана оптимальная трактовка понятия объекта – элемента вычислительных моделей. Эта трактовка позволяет быстро создавать, легко развивать и наглядно представлять не только сами модели, но и методы расчёта, на которых эти модели основаны. Предложенный общий подход проиллюстрирован на примере объектно-ориентированной структуры конкретной модели организма человека. Данная глава является подготовительной по отношению к главе 3 и не претендует на большую научную новизну, лишь уточняя понятия и методики, имеющиеся в литературных источниках и существующих средствах имитационного моделирования.

3. Объектно-ориентированные численные методы

3.1. Постановка задачи и обзор существующих работ

Конечно, математические свойства численных методов не зависят от способа их реализации, и поэтому предметом данной главы отнюдь не является улучшение этих свойств путём их программирования на объектно-ориентированном (а не на процедурном) языке. Объектно-ориентированный подход в вычислительной математике нужен не столько сам по себе, сколько ради применения сложных численных методов к имитационным моделям, необходимость чего детально описана в разделе 2.1.1. Вычислительные алгоритмы на объектно-ориентированной основе являются более гибкими, пригодными к визуальной разработке и повторному использованию в рамках некоторой библиотеки. Эти достоинства, которые свойственны имитационному подходу к моделированию, ориентированному на наглядность и простоту развития моделей, были обсуждены в главе 2. Ниже обсуждается другое, менее очевидное, преимущество ООП – удобство моделирования систем, которые сложны не в имитационном, а в математическом смысле.

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

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

Безусловно, за рубежом существует большое количество работ по анализу применимости объектно-ориентированных языков программирования к решению вычислительных задач [2,4,5]. Основной упор в этих работах делается на преодоление проблемы понижения скорости арифметических операций над большими массивами данных (векторами, матрицами) при использовании таких языков. Другими словами, речь идёт об универсальных рекомендациях вычислителям по эффективному программированию стандартных для процедурных языков конструкций на объектных языках. В частности, в [4] рассматриваются проблемы объектного представления матрицы (вектора), сплошной адресации элементов матриц, разделения данных между матрицей и её подматрицами, уменьшения числа вызовов конструкторов этих объектов. Во многих работах, например, в [2] приводятся результаты сравнения производительности объектных и процедурных языков, и обосновывается выбор наиболее производительных компиляторов. От объектно-ориентированного подхода в таких работах остаётся очень мало; наоборот, целью авторов является приблизить свойства объектно-ориентированных языков к процедурным (процедурные сложно использовать в задачах, связанных, например, с пользовательским интерфейсом). В частности, не делается попыток отказаться от матричного представления данных многих численных методов. Поэтому сами алгоритмы численных методов остаются неизменными, и научных работ по их объектной трактовке не возникает.

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

Естественным приложением научных работ данной тематики (в том числе и настоящей работы) является создание «объектно-ориентированных вычислительных библиотек». Несмотря на то, что их создано очень много, большинство предназначено для решения задач линейной алгебры, вычисления сложных функций и решения недифференциальных уравнений. Более сложные задачи, как обосновывается в [2], должны решаться самостоятельно пользователями этих базовых алгебраических библиотек. В данной работе предлагается иная концепция объектно-ориентированной библиотеки (см. 3.4).

3.2. Объектная интерпретация понятий вычислительной математики

3.2.1. Системы уравнений

В разделе 2.2.1 было обосновано, что элемент вычислительной модели должен представлять собой набор семантически связанных параметров, к которым привязаны расчётные алгоритмы. Другими словами, элемент трактовался как объект (в смысле ООП) и соответствовал либо типичным математическим объектам, либо типичным объектам предметной области. Ниже как раз обсуждается соответствие элемента универсальным математическим понятиям, без привязки к конкретной предметной области. Примером чисто математического элемента является система уравнений, а «семантически связанными параметрами», объединяемыми элементом, в этом случае являются искомые переменные и коэффициенты системы. Если коэффициенты системы – параметры элемента – не являются постоянными, а зависят от решения системы, то эту зависимость также можно представлять как объект или даже совокупность объектов (функций); однако такое представление описывалось выше и не имеет прямого отношения к рассматриваемому здесь математическому объекту.

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

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

Ещё одно математическое достоинство объектов заключается в их инвариантности относительно типа решаемой задачи. Является ли задача статической, динамической или вообще обратной задачей идентификации параметров модели – во всех этих случаях используются одни и те же расчётные классы элементов. Различие между этими задачами заключается только в алгоритме обновления состояния элементов, и этот алгоритм занимает всего несколько строчек кода. В случае динамической задачи обновление происходит вплоть до заданного времени, в случае статической – до получения заданной точности. В случае обратной задачи получение заданной точности происходит много раз в цикле, итерации которого соответствуют различным значениям искомых параметров модели, а изменение этих параметров происходит согласно какому-либо методу решения нелинейных уравнений.

Таким образом, основной элемент объектно-ориентированной модели – система уравнений – рассчитывает значения неизвестных системы (возможно, динамику их изменения во времени) по значениям своих коэффициентов (и правой части), взаимодействуя при этом с другими объектами, если коэффициенты или правая часть зависят от решения. Ниже принципиальные алгоритмы таких расчётов рассматриваются более подробно.

3.2.2. Схемы решения уравнений

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

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

Система уравнений, к какому бы типу она ни относилась, может быть аппроксимирована несколькими схемами (семействами схем), и с каждой схемой, в свою очередь, может быть ассоциировано несколько вычислительных алгоритмов. Как следствие, у соответствующего системе объекта (элемента) может быть несколько методов, которые по-разному решают одну и ту же задачу. Однако такой подход неудобен тем, что разные схемы могут иметь различный набор вспомогательных параметров, которые к тому же не должны иметь ничего общего с параметрами исходной системы уравнений. В связи с этим предлагается выделить вспомогательные параметры из элемента в отдельный объект – схему. Одновременно решается проблема слишком большого количества методов у элемента: различные варианты реализации каждой схемы программируются в соответствующих им объектах, а элементу остаётся лишь вызывать методы схемы, передавая им свои параметры. Таким образом, элемент, содержащий данные об исходной задаче, не хранит и не использует информацию о решающем эту задачу численном методе.

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

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

Проблема оптимизации решения больших систем уравнений путем их расщепления на несколько лучше обусловленных систем меньшей размерности пока не рассматривалась в вычислительной математике. Это связано с тем, что реализовать совместное решение этих систем с помощью процедурного подхода крайне затруднительно. Дело в том, что системы связаны между собой через общие данные, в то время как процедурным образом представленные алгоритмы их решения должны работать независимо. Использование объектного подхода позволяет избежать сложных приёмов разделения одних и тех же данных между процедурами решения разных систем: ничто не мешает одному и тому же параметру входить в несколько различных элементов – систем уравнений (а также в элементы других типов). Самый удобный среди процедурных приём разделения данных существует в языке FORTRAN (оператор EQUIVALENCE), однако он бессилен, если разные системы уравнений решает одна и та же процедура (проблема в том, что процедура, в отличие от класса, не может иметь нескольких экземпляров).

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

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

3.2.3. Пространственные сетки и сеточные шаблоны

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

Ниже для краткости сам узел пространственной сетки (а не система уравнений, соответствующая ему) фигурирует в качестве расчётного элемента. Связям между такими элементами соответствует отношение упорядоченности сеточных узлов по координатам (в случае использования неструктурированных сеток это есть отношение соседства сеточных узлов) – см. рис 3.1.

Упомянутая выше локальность связей означает, что значение какого-либо параметра в рассчитываемом сеточном узле определяется только значениями в узлах шаблона (например, ). Поэтому любой численный метод оказывается привязанным к узлу, то есть к элементу, и для его реализации данный узел не обязан иметь информацию о пространственно-удалённых узлах, с которыми у него нет связей. Это полностью соответствует объектной модели. Если же формулировать аппроксимацию уравнений в частных производных на процедурно-ориентированном языке системы линейных уравнений, то получается следующее: сначала искусственно создаётся проблема в виде огромной разреженной матрицы этой системы, а затем эта проблема преодолевается с помощью каких-нибудь специальных способов хранения и обработки этой матрицы в памяти компьютера.

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

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

Таким образом, при решении уравнений в частных производных ещё более эффективно, чем в других разделах вычислительной математики, используется объектно-ориентированный подход. Это проявляется как в более активном взаимодействии основных расчётных элементов (сеточных узлов), так и в появлении в дополнение к элементам и их параметрам ещё одного, более высокого, уровня объектов – уровня подсистем (областей интегрирования).

3.2.4. Ресурсоёмкость

Важным свойством, характеризующим любой численный метод, является его требовательность к машинным ресурсам. Эти ресурсы делятся на два основных типа – ресурсы процессорного времени и ресурсы оперативной памяти. Считается, что использование объектно-ориентированных языков программирования ведёт к большим затратам и того, и другого (по сравнению с процедурными языками) [4]. Следовательно, приемлемость объектов в качестве вычислительных единиц подлежит обоснованию.

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

Объектно-ориентированный подход может экономить память вычислительных программ за счёт совершенно других механизмов. Дело в том, что в процедурно-ориентированных программах часто встречаются массивы большой длины, заведомо достаточной для хранения всех параметров задачи, хотя часто выделенная под них память используется лишь на несколько процентов. Даже если язык поддерживает задание длины массива во время выполнения программы (в FORTRAN такая возможность появилась только в стандарте 1990-го г.), не всегда этим можно пользоваться, поскольку длина массива во время выполнения может неоднократно меняться. И уж совсем непреодолимые затраты памяти возникают при использовании многомерных массивов, когда для разных значений их медленно меняющегося индекса быстро меняющийся индекс должен изменяться в различных пределах. ООП решает эти проблемы автоматически, поскольку массив является объектом. В частности, элементы многомерного массива – тоже массивы – могут иметь различную длину. Помимо самих массивов, в объектно-ориентированные языки обычно включаются классы – оболочки массивов, которые могут динамически (оптимальным образом) менять их длину в зависимости от реального количества элементов в них. (Правда, динамическое изменение длины плохо отражается на быстродействии.)

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

Использование объектно-ориентированных языков увеличивает расходование ресурсов процессора, но в них существуют механизмы, максимально снижающие этот эффект. Действительно, наследование и полиморфизм уменьшают скорость вычислений, поскольку при вызове какого-либо метода у элемента происходит выбор среди одноимённых методов его класса и родительских классов (этот выбор называется динамическим, или поздним, связыванием). Однако быстродействие можно повысить, например, путём ещё одного наследования от данного класса и объявления полученного класса конечным (не способным иметь подклассов). В языке Java это позволяет компилятору не использовать динамическое связывание методов полученного класса, а также осуществлять автоматическую замену вызова методов их кодом, если он не слишком велик. В некоторых языках (С++) такая замена (inline-подстановка) может быть явно затребована программистом.

Следует также заметить, что с точки зрения быстродействия объектно-ориентированные численные методы имеют и свои преимущества. Дело в том, что в стандартной процедуре решения задачи математической физики значения параметров представляются в виде массивов большой длины (по числу сеточных узлов), доступ к которым осуществляется в одном или нескольких циклах по узлам. Поэтому время расходуется не только на расчёт, но и на сложение инкрементируемых переменных циклов с числами 1, -1, 2, -2 и т. д. (см. рис. 3.1) и на доступ к значениям параметров в массиве (вычисление адреса памяти по индексу массива). В случае объектной реализации численных методов в элементах хранятся явные ссылки на соседние элементы, у которых можно явно запросить значения их параметров, безо всякой арифметики адресов. При этом требуется больше памяти, но зато исключаются значительные затраты на целочисленную арифметику.

3.2.5. Некоторые специфические понятия

Наконец, необходимо рассмотреть то, как согласуются с объектно-ориентированным подходом менее универсальные понятия вычислительной математики, чем описанные выше понятия системы уравнений, схемы и сетки. Прежде всего, существуют многообразные свойства расчётных схем: явность, устойчивость, монотонность, консервативность, гибридность и т.д. Свойства устойчивости и монотонности являются чисто математическими и не имеют отношения к реализации схем; поэтому объектно-ориентированный подход не привносит в эти понятия ничего нового. Почти то же можно сказать о консервативности (свойство, имеющее смысл для гиперболических уравнений), однако практика показывает, что схемы, являющиеся удобными для своей объектной реализации, обычно одновременно являются консервативными. Это связано с тем, что объектные схемы всегда ближе к физике задачи, и реализующие их элементы соответствуют реальным объектам.

Гибридные схемы решения уравнений в частных производных (прежде всего гиперболических, поскольку именно они имеют требующие гибридизации разрывные решения) несколько проще реализуются с помощью объектно-ориентированного подхода, чем с помощью процедурного. Благодаря тому, что схема как объект в разделе 3.2.2 была отделена от самого расчётного элемента (в данном случае – от сеточного узла), для совместного использования в каждом узле двух схем, совершенно не нужно изменять ни код элемента, ни код схем. Достаточно создать класс гибридной схемы (один класс для всех возможных пар элементарных схем!), которая будет в соответствии с градиентом решения интерполировать результаты, получаемые по каждой из двух схем.

Несколько менее впечатляющим выглядит применение объектно-ориентированного подхода к реализации неявных схем. С одной стороны, самый простейший объектно-ориентированный численный метод для уравнений практически любого типа (который не следит за тем, на каком временном слое или на какой итерации берутся значения параметров) фактически является полуявным (треугольным), в то время как простейший процедурный численный метод – явный. С другой стороны, решение какой-либо вспомогательной системы уравнений на верхнем слое объектно-ориентированным способом достаточно сложно. В случае обыкновенных дифференциальных уравнений вспомогательная система уравнений нелинейная, и единственный разумный способ её решения состоит в привязке к основному расчётному элементу – к дифференциальной системе – ещё одного элемента (алгебраической системы). На каждом шаге своего решения дифференциальная система должна изменять некоторые коэффициенты алгебраической, после чего итерировать её схему до получения заданной точности.

Реализация неявных схем в случае пространственно-распределённых задач возможна без использования вспомогательных элементов, но лишь ценой существенного усложнения структуры и функций основных элементов – сеточных узлов. Большинство алгоритмов решения уравнений на верхнем временном слое не являются локальными, в отличие от самих схем для уравнений в частных производных, поэтому теряются некоторые преимущества объектного подхода по отношению к таким уравнениям. К примеру, алгоритм прогонки (даже монотонной прогонки) требует двукратного обращения к методам одного и того же сеточного узла: первый раз для вычисления прогоночных коэффициентов, а второй – для расчёта искомых параметров в этих узлах. Код метода прогонки получается более компактным и развиваемым, если выделить прогоночные коэффициенты в отдельный элемент, связанный с каждым сеточным узлом. Прямой ход прогонки состоит в цепочке вызовов методов расчёта именно этих вспомогательных элементов, и только на обратном ходу прогонки сами сеточные узлы последовательно вызывают друг у друга методы обновления своего состояния. Следует заметить, что распространение простейшего алгоритма прогонки на более сложные задачи (многомерные, с циклическими граничными условиями и т.д.) приводит к серьёзным трудностям, а объектно-ориентированная немонотонная прогонка вообще, по-видимому, невозможна. Если же решать уравнения на верхнем слое итерационными методами, то возникает необходимость в описанных выше вспомогательных элементах алгебраической системы уравнений, привязанных, в данном случае, к каждому сеточному узлу. Тем не менее в разделе 3.3.2.2 показывается, что и в случае неявных схем ООП может стимулировать развитие численных методов.

Выше речь шла о применимости объектов к некоторым типам схем; осталось рассмотреть их эффективность с точки зрения некоторых видов пространственных сеток. Отсутствие регулярной структуры сетки вполне соответствует объектному подходу, поскольку он (в отличие от процедурного) также не предполагает упорядочение узлов по координатам в массивах. Более того, так как для нерегулярных сеток характерно использование в разных узлах шаблонов разного размера и схем различного типа, объектный подход выглядит намного эффективнее процедурного в плане хранения коэффициентов этих схем. Он также облегчает алгоритмы поиска наилучших узлов шаблона.

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

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

Таблица 3 .1

Применимость ООП к некоторым понятиям вычислительной математики

3.3. Объектное представление конкретных численных методов

3.3.1. Обыкновенные дифференциальные уравнения

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

В принципе, при объектном представлении параметров очень удобно хранить в них любое требуемое число их значений на последовательных шагах по времени, что создаёт предпосылки для реализации многошаговых схем. Однако правая часть уравнений обычно состоит не только из параметров, но и из объектов-функций и объектов-операций, для которых хранение истории своего изменения во времени является избыточным. Поэтому целесообразно создавать специальные объекты (класса «параметр»), в каждом из которых хранить последовательность значений одного компонента вектора правой части. Кроме того, сложные задачи, для которых имеет смысл применять объектно-ориентированный подход, обычно приводят к жёстким системам уравнений, в связи с чем необходимо изменять шаг по времени в ходе расчётов, а к этому многошаговые схемы можно приспособить только за счёт существенного усложнения алгоритма. Такие задачи также часто содержат не только сосредоточенные, но и распределённые подсистемы, для расчёта которых многослойные схемы применяются редко. Поэтому использование многошаговых схем для ОДУ нецелесообразно и с точки зрения единообразия расчётов всех частей рассматриваемой системы.

С другой стороны, в пользу объектного представления таких схем можно привести тот факт, что по сравнению с процедурным представлением достаточно просто можно реализовать переход от одношаговых схем к многошаговым (это необходимо для «разгона» многошаговых схем на первых шагах) и обратный переход в процессе расчётов. Для этих переходов достаточно только поменять в объекте-системе ссылку на объект-схему. Следует также обратить внимание, что значения правой части системы ОДУ в случае использования многошаговых схем вычисляются при тех значениях аргументов, которые соответствуют решению системы: .

В случае одношаговых схем порядок их аппроксимации можно сделать выше первого только за счёт вычисления правой части при значениях аргументов, рассчитываемых по сложным формулам: , . При условии объектного представления аргументов (параметров) это приводит к необходимости присваивать им рассчитанные значения непосредственно перед вычислением правой части, после чего нужно «стирать» эти значения из истории изменения аргументов. К счастью, это можно делать чисто техническими способами (не пересматривая саму объектную модель): так как аргументы правой части являются одновременно параметрами элемента, представляющего систему уравнений, он может их изменять без особых проблем. В целом, использование объектно-ориентированных одношаговых методов решения ОДУ (схем типа Рунге-Кутты) предпочтительнее многошаговых (схем Адамса или Гира): не нужно «разгонять» метод на первых шагах и создавать специальные объекты для хранения динамики изменения правой части уравнения, легче изменять шаг по времени и связывать систему с другими расчётными элементами.

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

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

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

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

3.3.2. Уравнения в частных производных

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

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

3.3.2.1. Методы для гиперболических уравнений

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

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

Объектно-ориентированный подход позволяет существенно увеличить быстродействие за счет того, что поиск для каждого узла можно вести направленно, начиная с ближайших к нему соседей. Список соседних узлов, упорядоченный по расстоянию, строится только один раз, а в процессе расчётов у каждого элемента списка вызывается метод, вычисляющий расстояние от него до передаваемой ему в качестве аргумента базовой точки, а также вызывающий такой же метод у своих соседей. Если в этом методе проверять условия приближения узлов к базовой точке, то ( n+1) требуемых узлов находятся по кратчайшему пути. Конечно, и в процедурной программе для каждого узла можно хранить номера некоторого числа ближайших соседей (это число, кстати, фиксировано, в отличие от узлов-объектов). Однако подобный алгоритм в этом случае можно реализовать только через вызов сложной рекурсивной процедуры, требующей доступа к большому количеству не нужных ей данных и, главное, зависящей от размерности задачи. В данном же случае алгоритм поиска не зависит не только от размерности пространства, но и от численного метода (поскольку этот алгоритм содержится в элементе, а не в его схеме).

Схемы чётного (в частности, второго) порядка аппроксимации гиперболических уравнений обладают большой дисперсионной (и малой диффузионной) ошибкой, в связи с чем их решения сильно осциллируют на разрывах. Поэтому среди немонотонных схем имеет смысл рассматривать только схемы третьего порядка на шаблонах, содержащих на нижнем слое ( n+3) узлов. Ближе всего к монотонным находятся схемы, шаблон которых также лежит на минимальном расстоянии от базовой точки. Следовательно, описанный выше для схем первого порядка алгоритм быстрого поиска оптимального шаблона остаётся практически без изменений, а значит, элемент, который реализует наименее осциллирующую схему третьего порядка, можно унаследовать от элемента, реализующего наиболее точную монотонную схему. Как было показано в разделе 3.2.5, алгоритм гибридизации этих двух схем на основе объектного подхода является ещё более простым.

3.3.2.2. Методы для параболических уравнений

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

Рассмотрим сначала случай регулярной сетки. Если на верхнем слое аппроксимировать только значения пространственных производных искомой функции, то рассчитывать неявную схему можно методом прогонки, а в многомерном случае – также методом переменных направлений или методом дробных шагов с независимыми прогонками по каждой координате. С точки зрения ООП метод переменных направлений предпочтительнее метода многоточечной прогонки, поскольку позволяет использовать те же объекты схемы трёхточечной прогонки, что и в одномерных задачах. Для этого даже необязательно создавать специальный класс метода переменных направлений, поскольку n объектов метода трёхточечной прогонки имеют все необходимые данные (одну связь со своим сеточным узлом, две связи с такими же объектами-схемами для соседних узлов и набор коэффициентов схемы) и не должны быть связаны между собой. Таким образом, на этапе подготовки к расчётам «рядом» со структурой связанных внутренних сеточных узлов строится n аналогичных структур из объектов трёхточечной прогонки. Перед каждым шагом узел изменяет в соответствующем ему объекте метода лишь коэффициенты схемы (если коэффициенты уравнения переменные), а структура связей на протяжении всех расчётов остаётся неизменной. Граничные узлы (в которых заданы краевые условия первого рода) привязаны к объектам метода, соответствующим ближайшему внутреннему узлу; кроме того, на разных границах области интегрирования граничные узлы должны быть двух разных типов. Узлы первого типа инициируют прямой ход прогонки в цепочке объектов метода, а узлы второго типа инициируют обратный ход в цепочке внутренних узлов (после того, как до них дойдет очередь прямого хода прогонки – см. рис. 3.3).

В случае явных схем объект области интегрирования содержит неструктурированное множество сеточных узлов, которые он может в произвольном порядке переводить на следующий шаг по времени. Однако области интегрирования, предназначенная для узлов с неявными методами, должна отличать от остальных граничные узлы первого типа, самостоятельно переводя на следующий шаг только эти узлы. Для реализации метода переменных направлений такая область должна также иметь счётчик кратности (в двумерном случае – чётности) шага, а также несколько массивов граничных узлов первого типа (число массивов равно размерности пространства). Если коэффициенты и/или правая часть уравнения зависят от решения, то метод прогонки и метод переменных направлений можно дополнить дополнительными простыми итерациями на каждом (полном) шаге по времени (значения искомой функции на предыдущем временном слое являются хорошим начальным приближением, поэтому более сложные итерационные методы здесь избыточны). За эти итерации (происходящие вплоть до получения заданной точности) также отвечает объект области интегрирования.

Рис. 3 .3. Связи объектов, реализующих метод переменных направлений в случае 2D. Стрелки обозначают направление вычислений, пунктирные линии – связи по данным.

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

Для аппроксимации параболического уравнения с первым порядком в n -мерном случае требуется шаблон, содержащий (2 n+1) сеточных узлов (включая рассчитываемый). Если фиксировать этот шаблон (как это было сделано в случае регулярной сетки), то построение связей происходит только один раз перед расчётами. Если же коэффициенты уравнения сильно изменяются в ходе расчётов, то построение оптимального шаблона можно проводить на каждом шаге или через несколько шагов (как для нелинейных гиперболических уравнений), однако в силу резкого понижения быстродействия это делать нецелесообразно (см. также начало данного раздела). Поэтому каждый внутренний узел необходимо связать с 2 n ближайшими узлами, как и для гиперболических схем (см. 3.3.2.1), что позволяет использовать для «параболических» и «гиперболических» элементов один и тот же суперкласс. Последующий алгоритм связывания объектов схемы многоточечной прогонки, присоединенных ко внутренним узлам, более сложен. Для того, чтобы реализация метода не слишком зависела от размерности, предлагается создавать специальные объекты связей между объектами схемы, которые на этапе расчёта будут хранить соответствующие прогоночные коэффициенты. Эти связи в объекте схемы делятся на n входных и n выходных связей. Входные связи необходимы при вычислении прогоночных коэффициентов для выходных связей на прямом ходу прогонки, а выходные связи используются самими сеточными узлами для расчёта их значений на обратном ходу. Разделение на входные и выходные связи происходит автоматически в ходе процесса связывания, который инициируется граничным сеточным узлом первого типа (см. выше), беспорядочным образом распространяясь по сетке вплоть до граничных узлов второго типа. Следует заметить, что перед беспорядочным связыванием объектов схем внутри области интегрирования необходимо, чтобы все граничные узлы образовали связи с n приграничными объектами схемы; в противном случае граничный узел (как узел первого типа) может связаться с объектом схемы, который уже связал себя с данным узлом (как с узлом второго типа). Запрещение двунаправленных связей между граничным узлом и объектом схемы позволяет самостоятельно (без затухания) распространяться по сетке как процессу построения связей, так и процессам расчёта, основанным на этих связях. Тем не менее целесообразно проводить связывание объектов схемы более упорядоченно; а именно, нужно строить сначала связи с узлом, направление на который составляет наименьший угол с направлением на предыдущий узел цепочки (см. рис. 3.4). Чем более прямыми являются цепочки сеточных узлов, тем более устойчивым является метод прогонки.

Алгоритм перехода некоторой области интегрирования на один шага по времени, как было отмечено выше, состоит из нескольких этапов. В случае регулярной сетки для последовательного выполнения всех этих этапов область интегрирования на каждом шаге инициирует процесс расчёта только один раз (точнее, делает это в цикле, каждая итерация которого полностью рассчитывает одномерную цепочку объектов). В данном случае сначала граничными узлами первого типа рассчитываются все приграничные прогоночные коэффициенты, затем инициируется процесс расчёта объектами схемы всех остальных прогоночных коэффициентов в их связях, и только после его окончания через какой-либо граничный узел второго типа инициируется расчёт внутренних узлов (последние два процесса самостоятельно распространяются по сетке). Разделение на узлы первого и второго типа происходит автоматически при построении связей (если какой-либо объект схемы связался с некоторым граничным узлом, то он помечается как узел второго типа; в противном случае – остаётся узлом первого типа). Поэтому на первом этапе шага по времени можно не вычислять прогоночные коэффициенты в схемах, связанных с граничными узлами второго типа.