On numerical solution of the Cauchy problem for one type of ordinary differential equations

Full Text

Введение

Разработка специального математического и алгоритмического обеспечения, ориентированного на решение различных классов задач системного анализа, оптимизации принятия решений является важной проблемой и представляет значительный интерес [1]. Математический аппарат системного анализа состоит из достаточно большого набора инструментов, одним из которых являются обыкновенные дифференциальные уравнения, описывающие динамику изучаемых систем [2, 3]. При этом важна оценка эффективности решения задач системного анализа, основанных на использовании адекватных математических моделей. При их построении приходится учитывать многие параметры, что приводит к т.н. явлению жесткости. Его суть определяется тем, что для адекватного описания процессов в любой точке наблюдения приходится использовать функции, которые меняются достаточно быстро, и функции, которые меняются медленно, т.е. рассматриваемая физическая система характеризуется сильно различающимися характерными временами. Если не учитывать «малые» параметры при моделировании реальных процессов, данные обстоятельство может существенным образом повлиять на адекватность построенной модели [6]. Поэтому в системах уравнений, описывающих изучаемый процесс, необходимо учитывать большое число параметров, что приводит к явлению жесткости [5, 6]. Такие системы встречаются, например, при решении задач радиофизики, физики плазмы, биофизики, изучении нейронов [5, 6, 14]. Задачи, связанные с быстро затухающими переходными процессами, в которых жесткость проявляется естественным образом, охватывают самые разные области, включая изучение демпфирующих систем, анализ систем управления, задачи химической кинетики. При этом момент окончания процесса анализа системы (интегрирования системы дифференциальных уравнений) в таких задачах может определяться только в процессе решения системы, так как связан с ее асимптотическим поведением. Кроме того, при анализе реальных процессов на основе построенной модели в виде системы обыкновенных дифференциальных уравнений может происходить перестройка модели в связи с изменением параметров, от которых зависит изучаемый процесс. Это означает, что указанную систему приходиться решать многократно; при этом трудоемкость вычисления правых частей системы может быть достаточно высокой.

Соответственно, для эффективного решения задач системного анализа путем решения задачи Коши для системы обыкновенных дифференциальных уравнений необходимы эффективные методы численного интегрирования жестких систем дифференциальных уравнений. Необходимо отметить, что явление жесткости также характерно из-за большого различия собственных значений матрицы системы дифференциальных уравнений. Отличительной особенностью жестких систем является то, что на одних временных отрезках решение оказывается существенно нелинейным, а на других − близко к константе. Для решения такого рода систем используются неявные методы интегрирования [5, 6]. При этом, в силу существенно различного поведения решения на разных участках, требуется управлять выбором шагов интегрирования так, чтобы, с одной стороны, обеспечить приемлемую точность решения, а с другой стороны, когда трудоемкость вычисления правых частей системы высока, минимизировать число узлов интегрирования, используемой конечно-разностной формулы на всем отрезке интегрирования. Проблеме построения численных методов решения задачи Коши посвящено достаточно много работ [4-7].

В данной статье рассматриваются стратегии численного интегрирования систем обыкновенных дифференциальных уравнений, основанные на использовании конечно-разностной формулы [8]; приводится доказательство того, что задача о построении оптимального алгоритма удовлетворяет условиям многоэтапности, приведенным в [13], и описывается реализуемая версия оптимального алгоритма. При этом под оптимальным алгоритмом численного интегрирования понимается алгоритм, который минимизирует число вычислений правых частей системы дифференциальных уравнений при условии выполнения ограничений, определяемых точностью вычислений. Для реализации оптимального алгоритма требуется решение вспомогательных задач. В работе приведено исследование свойств оптимальной стратегии интегрирования, на основе которого предлагается модификация рассматриваемых стратегий. Она позволяет исследователю упростить реализацию, отказавшись от решения вспомогательных задач. При этом эксперименты показывают, что качество решения существенно не меняется, более того, в ряде случаев наблюдается уменьшение количества вычислений правых частей.

Статья является продолжением работ авторов [8-12].

Постановка задачи

Решается задача Коши для системы обыкновенных дифференциальных уравнений, при этом ее решение предполагается единственным и четырежды непрерывно дифференцируемым на отрезке $\left[x_0,t\right]$. Отметим, что значение параметра $t$ не является известным, а находится в процессе интегрирования системы. Для отыскания значений $Y_i=\left(y_1^i,\dots ,y_n^i\right)$, $i=\mathrm{1,2,}\dots$, решения системы будем использовать конечно-разностную формулу [8].

Пусть $X\left(\mathrm{0,}m\right)=\left(h_1,h_2,\dots ,h_{2m-1},h_{2m}\right)$ стратегия интегрирования по конечно-разностной формуле на отрезке $\left[x_0, x_0+z\right]$, $z$ – заданная положительная константа. При этом шаги интегрирования выбираются таким образом, что $h_{2i-1}={\tau }_i$, $h_{2i}={\tau }_i$, $i=\mathrm{1,}\dots ,m$, где ${\tau }_0$ – заданная константа.

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

$$\begin{gathered} \sum _{i=1}^m{\tau }_i\Rightarrow max, \\ {f}_i\left({\text{τ}}_{i-1},{\text{\hspace{0.17em}τ}}_i\right)\le \mathrm{0,} \\ {\text{τ}}_i\ge \mathrm{0,} i=\mathrm{1,}\dots ,m. \end{gathered}$$ (1)

Здесь ${\text{τ}}_0=\overline{{\text{τ}}_0}$ – положительная вещественная константа,

$f_i\left({\tau }_{i-1},{\tau }_i\right) = 14{\tau }_i^5+12{\tau }_{i-1}{\tau }_i^4+3{\tau }_{i-1}^2{\tau }_i^3-8{\text{Δ}}_i{\tau }_i-3{\text{Δ}}_i{\tau }_{i-1}$, (2)

где ${\text{Δ}}_i=2{\epsilon }_i{K}^{-1}$, ${\epsilon }_i>0$, $K>0$, $i=\mathrm{1,}\dots ,m$. Вещественная константа $K>0$ ограничивает абсолютное значение четвертой производной от компонент вектора решения 𝑌 (𝑥): $\underset{1\le j\le n}{\max }\left|y^{IV}{\left(x\right)}^j\right|\le K$,

${\epsilon }_i$ – заданная точность вычислений на $i$-ом шаге интегрирования.

Последовательность ${\text{τ}}_i$ будем называть оптимальной $\text{τ}-$последовательностью. Стратегию интегрирования $X\left(\mathrm{0,}m\right)$, полученную с помощью оптимальной $\tau -$последовательности, будем называть оптимальной стратегией интегрирования.

Построение оптимальной стратегии, основанной на использовании конечно-разностной формулы

Обозначим через

$u_i\left({\text{τ}}_{i-1},{\text{τ}}_i\right)=70{\text{τ}}_i^4+48{\text{τ}}_{i-1}{\text{τ}}_i^3+9{\text{τ}}_{i-1}^2{\text{τ}}_i^2-8{\text{Δ}}_i$ (3)

производную функции $f_i\left({\text{τ}}_{i-1},{\text{τ}}_i\right)$ по ${\text{τ}}_i$, а через

$v_i\left({\text{τ}}_{i-1},{\text{τ}}_i\right)=12{\text{τ}}_i^4+6{\text{τ}}_{i-1}{\text{τ}}_i^3-3{\text{Δ}}_i$ (4)

производную функции $f_i\left({\text{τ}}_{i-1},{\text{τ}}_i\right)$ по ${\text{τ}}_{i-1}$, $i=\mathrm{1,}\dots ,m$.

Лемма 1.

При выполнении условия

$f_i\left({\text{τ}}_{i-1},{\text{τ}}_i\right)=14{\text{τ}}_i^5+12{\text{τ}}_{i-1}{\text{τ}}_i^4+3{\text{τ}}_{i-1}^2{\text{τ}}_i^3-8{\text{Δ}}_i{\text{τ}}_i-3{\text{Δ}}_i{\text{τ}}_{i-1}=0$, (5)

где ${\text{τ}}_{i-1}\ge 0$, ${\text{τ}}_i>0$, $i=\mathrm{1,}\dots ,m$, имеют место следующие оценки:

$u_i\left({\text{τ}}_{i-1},{\text{τ}}_i\right)>0$, (6)

$v_i\left({\text{τ}}_{i-1},{\text{τ}}_i\right)>0$, (7)

 $0<1-\frac{v_i\left({\text{τ}}_{i-1},{\text{τ}}_i\right)}{u_i\left({\text{τ}}_{i-1},{\text{τ}}_i\right)}<1$. (8)

Доказательство.

Покажем справедливость неравенства (6), где ${\text{τ}}_{i-1}\ge 0$, ${\text{τ}}_i>0$. Учитывая соотношения (3) и (5) получаем, что:

${\text{τ}}_i{u}_i\left({\text{τ}}_{i-1},{\text{τ}}_i\right)=70{\text{τ}}_i^5+48{\text{τ}}_{i-1}{\text{τ}}_i^4+9{\text{τ}}_{i-1}^2{\text{τ}}_i^3-8{\text{Δ}}_i{\text{τ}}_i=f_i\left({\text{τ}}_{i-1},{\text{τ}}_i\right)++58{\text{τ}}_i^5+36{\text{τ}}_{i-1}{\text{τ}}_i^4+ 6{\text{τ}}_{i-1}^2{\text{τ}}_i^3+3{\text{Δ}}_i{\text{τ}}_{i-1}=58{\text{τ}}_i^5+36{\text{τ}}_{i-1}{\text{τ}}_i^4+6{\text{τ}}_{i-1}^2{\text{τ}}_i^3+3{\text{Δ}}_i{\text{τ}}_{i-1}>0$

В силу того, что ${\text{τ}}_i>0$, имеем $u_i\left({\text{τ}}_{i-1},{\text{τ}}_i\right)>0$.

Аналогично доказывается справедливость неравенства (7) путем рассмотрения выражения $\left(\frac{8}{3}{\text{τ}}_i+{\text{τ}}_{i-1}\right)v_i\left({\text{τ}}_{i-1},{\text{τ}}_i\right)$.

Учитывая (4) и (5), получаем, что:

$\left(\frac{8}{3}{\text{τ}}_i+{\text{τ}}_{i-1}\right)v_i\left({\text{τ}}_{i-1},{\text{τ}}_i\right)= 32{\text{τ}}_i^5+28{\text{τ}}_{i-1}{\text{τ}}_i^4++ 6{\text{τ}}_{i-1}^2{\text{τ}}_i^3-8{\text{Δ}}_i{\text{τ}}_i-3{\text{Δ}}_i{\text{τ}}_{i-1}=f_i\left({\text{τ}}_{i-1},{\text{τ}}_i\right)+ 18{\text{τ}}_i^5+16{\text{τ}}_{i-1}{\text{τ}}_i^4+3{\text{τ}}_{i-1}^2{\tau }_i^3=18{\text{τ}}_i^5+16{\text{τ}}_{i-1}{\text{τ}}_i^4++ 3{\text{τ}}_{i-1}^2{\text{τ}}_i^3>0$,

откуда следует, что $v_i\left({\text{τ}}_{i-1},{\text{τ}}_i\right)>0$ при ${\text{τ}}_{i-1}\ge 0$, ${\text{τ}}_i>0$ и $f_i\left({\text{τ}}_{i-1},{\text{τ}}_i\right)=0$, $i=\mathrm{1,}\dots ,m$.

Из (6) и (7) следует, что $1-\frac{v_i\left({\text{τ}}_{i-1},{\text{τ}}_i\right)}{u_i\left({\text{τ}}_{i-1},{\text{τ}}_i\right)}<1$.

Покажем теперь, что:

$1-\frac{v_i\left({\text{τ}}_{i-1},{\text{τ}}_i\right)}{u_i\left({\text{τ}}_{i-1},{\text{τ}}_i\right)}>0$. (9)

Так как $u_i\left({\text{τ}}_{i-1},{\text{τ}}_i\right)>0$, при ${\text{τ}}_{i-1}\ge 0$, ${\text{τ}}_i>0$ и $f_i\left({\text{τ}}_{i-1},{\text{τ}}_i\right)=0$, $i=\mathrm{1,}\dots ,m$, то (9) равносильно неравенству:

$w_i\left({\text{τ}}_{i-1},{\text{τ}}_i\right)=u_i\left({\text{τ}}_{i-1},{\text{τ}}_i\right)-v_i\left({\text{τ}}_{i-1},{\text{τ}}_i\right)>0$. (10)

Подставляя в (10) выражение для $u_i\left({\text{τ}}_{i-1},{\text{τ}}_i\right)$ и $v_i\left({\text{τ}}_{i-1},{\text{τ}}_i\right)$ из (3) и (4), получаем, что нужно убедиться в справедливости неравенства:

$w_i\left({\text{τ}}_{i-1},{\text{τ}}_i\right)=58{\text{τ}}_i^4+42{\text{τ}}_{i-1}{\text{τ}}_i^3+9{\text{τ}}_{i-1}^2{\text{τ}}_i^2-5{\text{Δ}}_i>0$

при выполнении равенства (5), где ${\text{τ}}_{i-1}\ge 0$, ${\text{τ}}_i>0$.

Для этого рассмотрим выражение $\left(\frac{8}{5}{\tau }_i+\frac{3}{5}{\tau }_{i-1}\right)w_i\left({\tau }_{i-1},{\tau }_i\right)$ и получим, что $\left(\frac{8}{5}{\tau }_i++\frac{3}{5}{\tau }_{i-1}\right)w_i\left({\tau }_{i-1},{\tau }_i\right)=f_i\left({\tau }_{i-1},{\tau }_i\right)+\frac{394}{5}{\tau }_i^5+\frac{276}{5}{\tau }_{i-1}{\tau }_i^4+\frac{183}{5}{\tau }_{i-1}^2{\tau }_i^3+\frac{27}{5}{\tau }_{i-1}^3{\tau }_i^2=\frac{394}{5}{\tau }_i^5++\frac{276}{5}{\tau }_{i-1}{\tau }_i^4+\frac{183}{5}{\tau }_{i-1}^2{\tau }_i^3+\frac{27}{5}{\tau }_{i-1}^3{\tau }_i^2>0$.

Следовательно, $w_i\left({\tau }_{i-1},{\tau }_i\right)>0$ и неравенство (9) доказано при ${\text{τ}}_{i-1}\ge 0$, ${\text{τ}}_i>0$ и $f_i\left({\text{τ}}_{i-1},{\text{τ}}_i\right)=0$, $i=\mathrm{1,}\dots ,m$. Таким образом, лемма 1 доказана.

Пусть ${\text{τ}}_0\in \left[\mathrm{0,}{\tilde{\tau }}_i\right]$, где ${\tilde{\tau }}_i={\left(\frac{4}{7}{\text{Δ}}_i\right)}^{\frac{1}{4}}$, $i=\mathrm{1,}\dots ,m$.

Лемма 2.

Уравнение $f_i\left({\tau }_{i-1},{\tau }_i\right)=0$ задает однозначную положительную функцию $g_i\left({\tau }_{i-1}\right)$ такую, что выполнение соотношения $f_i\left({\tau }_{i-1},{\tau }_i\right)\le 0$ равносильно выполнению неравенства ${\tau }_i\le g_i\left({\tau }_{i-1}\right)$, ${\tau }_{i-1}\in \left[\mathrm{0,}{\tilde{\tau }}_{i-1}\right]$, $i=\mathrm{1,}\dots ,m$.

Доказательство.

При ${\text{τ}}_{i-1}=0$ имеем $g_i\left(0\right)={\left(\frac{4}{7}{\text{Δ}}_i\right)}^{\frac{1}{4}}>0$.

Пусть при фиксированном значении ${\tau }_{i-1}\in \left[\mathrm{0,}{\tilde{\tau }}_{i-1}\right]{\tilde{f}}_i\left({\text{τ}}_i\right)=f_i\left({\text{τ}}_{i-1},{\text{τ}}_i\right)$, где функция $f_i\left({\text{τ}}_{i-1},{\text{τ}}_i\right)$ определена в (2). Используя (2), непосредственно находим, что функция ${\tilde{f}}_i\left({\text{τ}}_i\right)$ является непрерывной и выпуклой функцией при ${\text{τ}}_i\ge 0$. При этом ${\tilde{f}}_i\left(0\right)=-3{\text{Δ}}_i{\text{τ}}_{i-1}<0$, а ${\tilde{f}}_i\left({\tilde{\tau }}_i\right)=\frac{27}{7}{\text{Δ}}_i{\text{τ}}_{i-1}+3{\text{τ}}_{i-1}^2{\tilde{\tau }}_i^3>0$, так как ${\text{Δ}}_i>0$, ${\text{τ}}_{i-1}>0$, ${\tilde{\tau }}_i>0$, $i=\mathrm{1,}\dots ,m$. Следовательно, при каждом фиксированном значении ${\tau }_{i-1}\in \left[\mathrm{0,}{\tilde{\tau }}_{i-1}\right]$ существует единственное значение ${\text{τ}}_i=g_i\left({\text{τ}}_{i-1}\right)\in \left(\mathrm{0,}{\tilde{\tau }}_i\right)$ такое, что ${\tilde{f}}_i\left(g_i\left({\text{τ}}_{i-1}\right)\right)=f_i\left({\text{τ}}_{i-1},g_i\left({\text{τ}}_{i-1}\right)\right)$, причем ${\tilde{f}}_i\left({\text{τ}}_i\right)=f_i\left({\text{τ}}_{i-1},{\text{τ}}_i\right)\le 0$, для каждого ${\tau }_{i-1}\in \left[\mathrm{0,}{\tilde{\tau }}_{i-1}\right] ,{\text{ τ}}_i\in \left[\mathrm{0,}{g}_i\left({\text{τ}}_{i-1}\right)\right]$, $i=\mathrm{1,}\dots ,m$, откуда следует равносильность соотношений $f_i\left({\text{τ}}_{i-1},{\text{τ}}_i\right)\le 0$ и ${\text{τ}}_i\le g_i\left({\text{τ}}_{i-1}\right)$.

Лемма 2 доказана.

Теорема 1.

Решение задачи (1) существует.

Доказательство.

Приведем задачу (1) к задаче на минимум:

$$\begin{gathered} -\sum _{i=1}^m{\tau }_i\Rightarrow min \\ f_i\left({\text{τ}}_{i-1},{\text{τ}}_i\right)\le 0, i=\mathrm{1,}\dots ,m, \\ {\tau }_i\ge 0, i=\mathrm{1,}\dots ,m, \end{gathered}$$ (11)

${\text{τ}}_0=\overline{{\text{τ}}_0}>0$ - заданная величина.

Пусть $D_m$ − допустимое множество задачи (11). Множество $D_m$ является ограниченным. Действительно, при заданном значении ${\text{τ}}_0\in \left(\mathrm{0,}\widetilde{{\text{τ}}_0}\right)$, с учетом доказательства леммы 2, имеем, что если выполнено соотношение $f_1\left({\text{τ}}_0,{\text{τ}}_1\right)\le 0$, то $0\le {\text{τ}}_1\le g_1\left({\text{τ}}_0\right)$. В силу того, что функция $g_i\left({\text{τ}}_{i-1}\right)$ строго монотонно убывает при $f_i\left({\text{τ}}_{i-1},{\text{τ}}_i\right)\le 0$ из леммы 2 следует, что $0\le {\text{τ}}_i\le g_i\left(0\right)={\left(\frac{4}{7}{\text{Δ}}_i\right)}^{\frac{1}{4}}$, $i=\mathrm{2,}\dots ,m$. Поскольку функции $f_i\left({\text{τ}}_{i-1},{\text{τ}}_i\right)$, $i=\mathrm{1,}\dots ,m$ непрерывны, то множество $D_m$ замкнуто. Тогда, по теореме Вейерштрасса, задача (1) имеет решение. Теорема 1 доказана.

В силу леммы 1, леммы 2 и теоремы 1 задача (1) является многоэтапной задачей математического программирования [14]. Это означает, что задача (1) может быть решена последовательным решением $m$ локальных задач:

$$\begin{gathered} {\tau }_i\Rightarrow max, \\ {f}_i\left({\tau }_{i-1},{\tau }_i\right)\le 0, i=\mathrm{1,}\dots ,m. \end{gathered}$$ (12)

Таким образом, чтобы построить оптимальную стратегию интегрирования, необходимо решить уравнение (5) для нахождения значения ${\text{τ}}_i$, $i=\mathrm{1,}\dots ,m$.

Исследование оптимальной стратегии

В [16] рассматривалась упрощенная стратегия интегрирования, основанная на том свойстве, что оптимальная $\tau$-последовательности сходится к неподвижной точке при фиксированном значении ${\text{Δ}}_i=\text{Δ }$ на отрезке длины $z$. Приведем доказательство данного факта.

Напомним, что уравнение (5) на отрезке $\left[\mathrm{0,}\tilde{\tau }\right]$, $\tilde{\tau }={\left(\frac{4}{7}\text{Δ}\right)}^{\frac{1}{4}},\tilde{\tau }<z$, при фиксированном ${\tilde{\tau }}_{i-1 }$имеет единственный положительный корень. Пусть $g$ - отображение отрезка $\left[\mathrm{0,}\tilde{\tau }\right]$ в себя, определяемое уравнением (5), тогда $g\left(\overline{\tau }\right)$ будет единственным положительным корнем уравнения (5), где $\overline{\tau }\in \left[\mathrm{0,}\tilde{\tau }\right]$.

Производная $g^{\text{'}}\left(\overline{\tau }\right)$ функции $g\left(\overline{\tau }\right)$, $\overline{\tau }\in \left[\mathrm{0,}\tilde{\tau }\right]$, имеет вид:

$g^{\text{'}}\left(\overline{\tau }\right)=-\frac{v\left(\overline{\tau },g\left(\overline{\tau }\right)\right)}{u\left(\overline{\tau },g\left(\overline{\tau }\right)\right)}$,

где $u\left(\overline{\tau },g\left(\overline{\tau }\right)\right)=70g^4\left(\overline{\tau }\right)+48g^3\left(\overline{\tau }\right)+9\overline{{\text{τ}}^2}{g}^2\left(\overline{\tau }\right)-8$, $v\overline{\tau }\left(,g\left(\overline{\tau }\right)\right)=12g^4\left(\overline{\tau }\right)+6g^2\left(\overline{\tau }\right)-3$. При выполнении равенства (5)

$u\left(\overline{\tau },g\left(\overline{\tau }\right)\right)>0$,

$v\left(\overline{\tau },g\left(\overline{\tau }\right)\right)>0$,

т.е. функция $g\left(\overline{\tau }\right)$ дифференцируема при $\overline{\tau }\in \left[\mathrm{0,}\tilde{\tau }\right]$ и строго монотонно убывает на этом отрезке. Кроме того

$\left|g^{\text{'}}\left(\overline{\tau }\right)\right|\le \frac{3}{16}$

в силу (5), т.е. отображение $g$ на отрезке $\left[\mathrm{0,}\tilde{\tau }\right]$ является сжимающим отображением и имеет единственную неподвижную точку

${\overline{\tau }}^{*}={\left(\frac{11}{29}\text{Δ}\right)}^{\frac{1}{4}}$.

Исследуем зависимость числа итераций $\overline{m}$, необходимых для достижения неподвижной точки ${\overline{\tau }}^{*}$, от значения $∆$ и ${\tau }_0$. Для этого будем разрешать уравнение (5) методом дихотомии с заданной точностью $\delta =\mathrm{0,000001}$. В результате вычислительных экспериментов на тестируемых системах дифференциальных уравнений выявлены характерные интервалы изменения значений для $∆$ и ${\tau }_0$: [10^-15, 10⁵] и [10^-6, 10¹] соответственно. Вычислительный эксперимент показал, число итераций $\overline{m}$ при варьировании значения ${\tau }_0$ для фиксированного значения $∆$ слабо изменяется. Для примера приведем в табл. 1 результаты эксперимента при $\text{Δ}=\mathrm{0,0001,}{\overline{\tau }}^{*}=\mathrm{0,078478.}$

Таблица 1. Сходимость оптимальной τ-последовательности при фиксированном Δ = 0,0001

Table 1. The convergence of optimal τ-sequence for fixed Δ = 0.0001

При этом при уменьшении значения $∆$ число итераций $\overline{m}$ уменьшается. В табл. 2 приведены результаты эксперимента для фиксированного значения шага ${\tau }_0 = \mathrm{0,0001}$. Данный факт позволяет предложить использовать для вычисления шагов интегрирования значения неподвижной точки на подынтервалах с быстрым изменением решения (например, при $\text{Δ}\le {10}^{-7}$), т.е. не использовать вспомогательные методы для решения уравнения (5) для указанных подынтервалов.

Таблица 2. Сходимость оптимальной τ-последовательности при фиксированном τ₀ = 0,0001

Table 2. The convergence of optimal τ-sequence for fixed τ₀ = 0.0001

Реализуемая версия оптимального алгоритма

Общую схему алгоритма интегрирования, основанного на оптимальной стратегии, можно описать следующим образом.

1. Определим начальный шаг интегрирования $h_0={\text{τ}}_0$, где ${\text{τ}}_0$ − заданное положительное число.
2. Определим величину $z > \mathrm{0,} z \in R$, определяющая размер подотрезка $\left[x_0,x_0+z\right]$, на котором строится решение системы.
3. Построим оценку $K$ четвертой производной от компонент решения системы дифференциальных уравнений для подотрезка $\left[x_0,x_0+z\right]$. По трем узлам из истории вычислений, если у нас есть информация о предыдущих узлах, и узлу сделанного с шагом $h=z$ построим полином Лагранжа четвертого порядка. Далее вычислим четвертую производную для данного полинома и в качестве оценки $K$ возьмем ее абсолютное значение. Если информации о предыдущих вычислениях нет (в начале работы метода), необходимо трижды проинтегрировать систему с шагом $h=-{\tau }_0$.
4. Определим шаги интегрирования $h_{2i-1}={\tau }_i$, $h_{2i}={\tau }_i$, а $h_{2i-2}={\tau }_{i-1}$, $i=\mathrm{1,2,} \dots$. Значение ${\tau }_i$ находим, применяя одну из возможных модификаций: разрешая уравнение (5) при фиксированном ${\tau }_{i-1}$, если $∆$ превышает пороговое значении $\overline{∆}$, и используя в качестве величины шага значение ${\overline{\tau }}^{*}$ в противном случае.
5. Если $x_m>x_0+z$, переходим к следующему подотрезку, принимая $x_0=x_m$, и переходим к шагу 2.
6. Процесс продолжается до тех пор, пока не будет выполнено условие останова.

Вычислительный эксперимент

Продемонстрируем работоспособность предложенной модификации для стратегии интегрирования, описанной в [10, 11] и основанной на оптимальной стратегии. Эксперименты проводились для обыкновенных систем уравнений, описывающих модель Ходкина-Хаксли [14]. Данная модель используется в нейродинамике в качестве фундаментальной модели описания динамики нейрона и расширяется, например, добавлением в модель дополнительных ионных каналов и учетом пространственной структуры нейрона.

Приведем здесь результаты экспериментов, где в качестве функции $I\left(t\right)$, определяющей силу тока, подающуюся на мембрану нейрона, выбирались константная $I\left(t\right) = \mathrm{0,5}$ и синусоидальная функции $I\left(t\right) = sin t+1$, отрезок интегрирования выбирался равным [0, 200]. В табл. 3 представлены результаты, отражающие сравнение запусков с применением предложенной модификации и оптимальной стратегии выбора шагов интегрирования. В качестве основных параметров сравнения используется количество шагов $N$, минимальная и максимальная величины шагов $h_{min}$ и $h_{max}$ при варьировании параметра ${\epsilon }_i$, определяющего точность вычисления размера шага. Оптимальная стратегия обозначена через $\alpha$, предложенная модификация – через $\overline{\alpha }$ соответственно. Для $\overline{\alpha }$ пороговое значение $\overline{∆}={10}^{-7}$.

Таблица 3. Сравнение результатов запуска $\mathit{\alpha }$ и $\overline{\mathbf{\alpha }}$

Table 3. The comparison of obtained results for $\mathit{\alpha }$ and $\overline{\mathbf{\alpha }}$

Отметим, что для более низкой точности ${\epsilon }_i$ предложенная модификация не запускается, так как значение параметра $∆$, которая зависит от ${\epsilon }_i$ не попадает в определенный выше диапазон. На рис. 1 и 2 изображены графики компонент решений V(t), m(t), n(t), h(t), полученные для синусоидальной функции для случая ε_i = 0,00001.

Рис. 1. График функции V(t) для случая ε_i = 0,00001

Fig. 1. V(t) function graph for ε_i = 0.00001

Рис. 2. Графики функций m(t), n(t), h(t) для случая ε_i = 0,00001

Fig. 2. m(t), n(t), h(t) functions graphs for ε_i = 0.00001

В табл. 4 представим результаты запусков алгоритма для задачи Коши, приведенной в [8]. Она имеет аналитическое решение, поэтому в качестве основных параметров сравнения будем использовать максимальную разность между аналитическим и численным решением $\sigma$, полученную по всем компонентам решения, и число шагов интегрирования N. Оптимальная стратегия обозначена через $\alpha$, предложенная модификация – через $\overline{\alpha }$ соответственно. Здесь в модификации $\overline{\alpha }$ будем использовать для подотрезков длиной z=0,1 величину шага, равную ${\overline{\tau }}^{*}$.

Таблица 4. Сравнение результатов запуска $\mathit{\alpha }$ и $\overline{\mathbf{\alpha }}$

Table 4. The comparison of obtained results for $\mathit{\alpha }$ and $\overline{\mathbf{\alpha }}$

Заключение

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

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

Исследования проведены с использованием системы тестирования разрабатываемых алгоритмов [15].

Работа выполнена в рамках реализации программы «Передовая инженерная школа» Нижегородского государственного университета им. Н.И. Лобачевского.

About the authors

A. G. Korotchenko

Author for correspondence.
Email: koangr@yandex.ru
ORCID iD: 0000-0002-7150-8114
Russian Federation, Nizhny Novgorod

V. M. Smoryakova

Email: smorykov@mail.ru
ORCID iD: 0000-0003-0357-0976
Russian Federation, Nizhny Novgorod

References

Supplementary files

Supplementary Files

Action

1. JATS XML

Download

2. Fig. 1. V(t) function graph for εi = 0.00001

Download (121KB)

Indexing metadata

3. Fig. 2. m(t), n(t), h(t) functions graphs for εi = 0.00001

Download (133KB)

Indexing metadata