ров наборов экспериментов. Для каждого объекта системы при его создании генератором указывается произвольное количество метаданных (параметров), необходимое как для сохранения ключевых событий при работе алгоритма, так и для построения агрегированных результатов. Например, такими параметрами могут быть точность работы схемы, название и параметры выбранной стратегии, название вспомогательных алгоритмов и значения их параметров. Таким образом, на первом этапе происходит формирование экспериментов с учетом варьирования интересующих оператора параметров в выбранном им диапазоне.
Второй этап - параллельное выполнение экспериментов и получение первичной информации о результатах, а также накопление ключевых событий. Отметим тот факт, что управление условием останова описывается вне метода интегрирования, другими словами метод численного интегрирования описан как генератор бесконечной последовательности узлов интегрирования. Данное условие позволяет легко включать одну схему интегрирования как вспомогательную для другой. В контексте каждого объекта системы для данного этапа доступен журнал событий, в который записываются интересующие оператора события, что позволяет впоследствии проанализировать поведение алгоритма над конкретной задачей Коши.
Третий этап заключается в построении агрегированной информации и определении отличающихся параметров, для представления в консолидированном отчёте. Здесь, как и на предыдущих этапах, может выполнятся пользовательский код, определенный оператором для дополнительного анализа данных. На этом этапе также генерируются все связанные артефакты (графики полученных решений, события, накопленные в эксперименте и т.п.) специфичные как для отдельных экспериментов, так и всех экспериментов, проведённых для конкретной задачи Коши или с использованием конкретного метода численного интегрирования.
В систему включены реализации некоторых базовых методов интегрирования, например, метод Эйлера и методы Рунге-Кутты разных порядков, реализованы необходимые вспомогательные алгоритмы.
Ввиду модульного подхода программная система не ограничивает оператора использованием только имеющихся схем интегрирования. В частности, в системе реализован адаптер для выполнения методов интегрирования описанных в системе MATLAB над произвольной задачей Коши, заданной оператором, наравне с методами интегрирования, построенными на основе работ [1-5]. Адаптер функционирует по клиент-серверной модели: для выполнения схемы интегрирования, описанной в среде MATLAB, программная система запускает временный TCP сервер на базе технологии WCF (Windows Communication Framework), загружает среду выполнения MATLAB, в которую загружается пользовательский модуль .Net, выступающий в роли TCP клиента. Данный модуль выступает посредником, позволяющим переадресовать запрос на вычисление значения системы в указанном узле из среды MATLAB в программную систему. Таким образом, у программной системы появляется возможность с
90
одной стороны предоставлять среде MATLAB произвольную выбранную оператором систему уравнений, а с другой стороны вносить в журнал события уникальные вычисления системы.
Литература
1. Коротченко А.Г. О задачах математического программирования, имеющих многоэтапный характер [Текст] / Коротченко А.Г. // Вестник нижегородского государственного университета им. Н.И. Лобачевского. – 2011. – Т. 1.
–С. 183-187.
2.Korotchenko A.G. On a method of construction of numerical integration formulas [Text] / Korotchenko A.G., Smoryakova V.M. // AIP Conference Proceedings. – 2016. – V. 1776 №9780735414389. – P.090012-1 - 090012-4.
3.Korotchenko A.G. On a Comparison of Several Numerical Integration Methods for Ordinary Systems of Differential Equations [Text] / Korotchenko A.G.,
Smoryakova V.M. // Lecture Notes in Computer Science.– 2020. – V. 2. № 11974. –
P.406-412.
4.Коротченко А.Г. Введение в многокритериальную оптимизацию. [Текст] / Коротченко А.Г., Кумагина Е.А., Сморякова В.М. // Н.Новгород: Изда-
тельство ННГУ. – 2017. – 55 c.
5. Коротченко А.Г. О сравнении различных методов решения задачи Коши для системы обыкновенных дифференциальных уравнений [Текст] / Коротченко А.Г., Сморякова В.М. // Интеллектуальные информационные системы Труды Международной научно-практической конференции. В 2-х частях. - 2019. - С. 133-137.
6.Мартин Р. Чистый код: создание, анализ и рефакторинг. Библиотека программиста. – СПб.: Питер, 2013. — 464 с.
7.Себеста Р.У. Основные концепции языков программирования. - 5-е изд.
-М.: Вильямс, 2001. – 670 c.
Национальный исследовательский Нижегородский государственный университет им. Н.И. Лобачевского
УДК 519.63
А. Г. Коротченко, В. М. Сморякова
О СТРАТЕГИЯХ ЧИСЛЕННОГО ИНТЕГРИРОВАНИЯ, ОСНОВАННЫХ НА ИСПОЛЬЗОВАНИИ ОДНОЙ КОНЕЧНО-РАЗНОСТНОЙ ФОРМУЛЫ
Проблема построения алгоритмов с заданными свойствами является важной проблемой как с практической, так и с теоретической точек зрения. Указанной проблеме посвящены многие работы такие, например, как [1-6]. При
91
этом соответствующий алгоритм может конструироваться, исходя из некоторой формальной модели вычислений в том числе с учётом его оптимальности по тому или иному критерию на рассматриваемом классе задач. Кроме того, при реализации построенного таким образом алгоритма он может дополняться процедурами, которые улучшают его свойства.
В статье рассматривается алгоритм численного интегрирования систем обыкновенных дифференциальных уравнений построенный в соответствии с описанной выше схемой, см. [1, 3, 4]. Данная работа основана на статье [7].
Пусть требуется найти решение Y =Y (x) (предполагается, что решение Y(x) - единственно и четырежды непрерывно дифференцируемо на отрезке [x0,t]) задачи Коши для системы обыкновенных дифференциальных уравнений:
Y'= F(x,Y), Y(x0) =Y0,
Y=Y(x) = (y1(x),..., yn(x)),x0 ≤ x < t,
где значение параметра t неизвестно априори, а определяется в процессе интегрирования системы. В случае, когда система линейна F(xY,)
= B(x)Y +b(x), где B(x) - n(×n)- матрица, b(x) - n -мерный вектор.
Для определения численных значений Yi = (y1i,..., yni ), i =1,2..., решения сис-
темы (1) будем использовать метод, основанный на применении конечноразностной формулы четвертого порядка [3].
Для того, чтобы реализовать алгоритм, описанный в [3] необходимо определить стратегию перехода от одного отрезка длины z к другому, вычисляя для указанного отрезка оценку четвёртой производной K от компонент решения системы дифференциальных уравнений (1).
Будем определять значение K , используя в качестве аппроксимации решения интерполяционный полином Лагранжа четвёртой степени. При определении решения в дополнительных точках можно использовать какой-либо явный метод. Здесь мы используем метод Рунге-Кутты второго порядка.
Для отрезка [x0,x0 + z] определим методом Рунге-Кутты второго порядка с
шагом h = −τ0 решение Y−3, Y−2, Y−1 в узлах x−3 , x−2 x−1 и с шагом h = z решение Y(x0 + z) в узле x0 + z , где z и τ0 - заданные положительные константы.
Используя узлы x−3 , x−2 x−1, x0 и x0 + z построим полином Лагранжа четвёртого порядка для каждой компоненты решения. Здесь Y0 - начальное условие в узле x0 . В качестве оценки K для отрезка [x0,x0 + z] будем использовать мак-
симальное значение четвертой производной по всем компонентам решения от построенного полинома Лагранжа четвёртого порядка.
Следующий узел интегрирования xi = xi−1 +hi , i =1,2,..., при этом hi опре-
деляется с помощью алгоритма описанного в [3]. Определим несколько процедур перехода к новому подотрезку.
92
Пусть следующий предполагаемый узел xi > x0 + z, тогда выберем в качестве hi = z− xi−1 и найдём решение системы Y(x0 + z) с использованием конечно-
разностной формулы, описанной в [3]. Далее переходим к новому отрезку, где x0 = x0 + z и перезапускаем описанную стратегию. Обозначим данную стратегию
через β1, см. [1]. Здесь выбор величины шага ограничен сверху значением па-
раметра z, что не позволяет использовать оптимальную стратегию выбора шага интегрирования [3] в случае, когда величина шага полученного в ходе интегрирования больше z.
Рассмотрим следующую модификацию указанной стратегии и обозначим её через β2 . Пусть опять узел xi > x0 + z, тогда найдём решение системы Yi(xi) с
использованием конечно-разностной формулы, описанной в [3]. Теперь осуществим переход к новому отрезку пусть x0 = xi и z = max(z,hi).
Обе описанные стратегии β1 и β2 обладают тем недостатком, что прихо-
дится на каждом подотрезке дополнительно вычислять решение системы в вспомогательных узлах, для получения оценки K .
Если рассматриваемая система обладает свойством затухания, для неё можно предложить процедуру перехода β3, описанную в [4], в которой число
вспомогательных вычислений сокращено в сравнении с процедурами β1 и β2 . Пусть K = K0 , значениеK0 определим как описывалось выше для отрезка
[x0,x0 + z]. Для узла xi > x0 + z найдём значение оценки четвёртой производной
решения Ki в узлах (xi−4,Yi−4), (xi−3,Yi−3) , (xi−2,Yi−2), (xi−1,Yi−1) и (xi,Yi) с помощью интерполяционного полинома Лагранжа четвертой степени. Если K ≥ Ki , то
примем K = Ki 2+ K . В противном случае, если K < Ki , переходим к новому от-
резку, где в качестве x0 выбирается значение xi , а z = max(z,hi) и так далее, пока
не выполнится условие окончания процесса интегрирования.
В вычислительном эксперименте рассматривается сравнение описанных выше стратегий. При сравнении использовались три критерия: первый - число
узлов интегрирования N , второй – максимальная разность ∆ между точным и приближённым решениями, определяемая по всем узлам интегрирования и всем компонентам решения, третий M - число вычислений правых частей системы.
В качестве тестовых (модельных) задач использовались системы, начальные условия которых были подобраны таким образом, что отрезок интегрирования разбивается на два интервала, в первом из которых каждая компонента вектора решения существенно не линейна, а на втором она почти не изменяется. Приведём здесь результаты эксперимента (табл. 1-3 ) для системы линейных дифференциальных уравнений, описанной в [1], которая рассматривалась на отрезке [0,5]. Для определения численных значений решения используется ме-
93
тод Гаусса. Здесь начальная величина z = 0.1, начальный шаг h0 = 0.000001, значение параметра εi =ε .
Таблица 1
ε |
N(β1) |
N(β2) |
N(β3) |
0.1 |
53 |
19 |
26 |
0.01 |
61 |
32 |
37 |
0.001 |
77 |
55 |
57 |
0.0001 |
109 |
100 |
91 |
0.00001 |
191 |
187 |
152 |
Таблица 2
ε |
∆(β1) |
∆(β2) |
∆(β3) |
0.1 |
0.05 |
0.01 |
0.007 |
0.01 |
0.003 |
0.003 |
0.002 |
0.001 |
0.0005 |
0.0006 |
0.0004 |
0.0001 |
0.0002 |
0.0001 |
0.0001 |
0.00001 |
0.00002 |
0.00001 |
0.00003 |
Таблица 3
ε |
M(β1) |
M(β2) |
M(β3) |
0.1 |
110 |
34 |
33 |
0.01 |
118 |
50 |
44 |
0.001 |
134 |
78 |
64 |
0.0001 |
166 |
129 |
98 |
0.00001 |
248 |
225 |
159 |
На рисунке представлены результаты запусков стратегий β1, обозначенные маркерами круглой формы, β2 , обозначенные маркерами треугольной формы и β3, обозначенные маркерами ромбовидной формы соответственно.
Координаты маркеров соответствуют значениям критериев ∆ и M . Для наглядности данных значения максимальной разности ∆ между точным и приближённым решениями на рисунке отображены на логарифмической шкале.
94