Численная схема для одной интегро-дифференциальной системы, связанной с задачей космического динамо

Полный текст

Введение

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

Феноменологически эредитарная модель космического динамо может быть представлена в интегро-дифференциальном виде[1, 2]:

$\begin{array}{l}\frac{dx}{dt}\mathrm{<mo>=</mo>}-\mu x+\mathrm{<mo\; stretchy="false">(</mo>}\mu -\varphi {\displaystyle {\int }_0^t}K\mathrm{<mo\; stretchy="false">(</mo>}t-\tau \mathrm{<mo\; stretchy="false">)</mo>}Q\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">(</mo>}\tau \mathrm{<mo\; stretchy="false">)</mo>,}y\mathrm{<mo\; stretchy="false">(</mo>}\tau \mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}d\tau \mathrm{.<mo\; stretchy="false">)</mo>}y\mathrm{,}\hfill \\ \frac{dy}{dt}\mathrm{<mo>=</mo><mo\; stretchy="false">(</mo>}D-{\displaystyle {\int }_0^t}K\mathrm{<mo\; stretchy="false">(</mo>}t-\tau \mathrm{<mo\; stretchy="false">)</mo>}Q\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">(</mo>}\tau \mathrm{<mo\; stretchy="false">)</mo>,}y\mathrm{<mo\; stretchy="false">(</mo>}\tau \mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}d\tau \mathrm{.<mo\; stretchy="false">)</mo>}x-y\mathrm{.}\hfill \end{array}$ (1)

В работе будем использовать компактную форму записи системы (1), которая имеет вид:

$\begin{array}{ll}\frac{dx}{dt}\mathrm{<mo>=</mo>}\left(\eta -p{s}^{2}z\right)y-\eta x\mathrm{,}\hfill & \hfill \\ \frac{dy}{dt}\mathrm{<mo>=</mo><mo\; stretchy="false">(</mo>}D-z\mathrm{<mo\; stretchy="false">)</mo>}x-y\mathrm{,}\hfill & \hfill \\ z\mathrm{<mo>=</mo>}{\displaystyle {\int }_0^t}K\mathrm{<mo\; stretchy="false">(</mo>}t-\tau \mathrm{<mo\; stretchy="false">)</mo>}Q\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">(</mo>}\tau \mathrm{<mo\; stretchy="false">)</mo>,}y\mathrm{<mo\; stretchy="false">(</mo>}\tau \mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}d\tau \mathrm{.}\hfill & \hfill \end{array}$ (2)

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

Рассмотрим задачу Коши для дифференциального уравнения

$\frac{dx}{dt}\mathrm{<mo>=</mo>}f\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>,}t\mathrm{<mo\; stretchy="false">)</mo>,}x\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo><mo>=</mo>}a\mathrm{.}$

Любое такое дифференциальное уравнение может быть записано в эквивалентном в виде

$x\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}a+{\displaystyle {\int }_0^t}f\mathrm{<mo\; stretchy="false">(</mo>}\tau \mathrm{,}x\mathrm{<mo\; stretchy="false">(</mo>}\tau \mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}d\tau \mathrm{.}$

Такая конструкция представляет из себя, как говорилось ранее, частный случай интегрального уравнения Вольтерра [4-6]. В общем виде интегрального уравнения Вольтерра предполагается, что подынтегральная функция $f$ зависит не только от $\tau$ но и от $t$. Если функция $f$ зависит только от $\tau$, то можно свести такое интегральное уравнение к дифференциальному уравнению. Если же функция $f$ явно зависит не только от $\tau$, но и от $t$, то вообще говоря, не факт, что можно свести интегральное уравнение к дифференциальному уравнению или системе дифференциальных уравнений.

Рассмотрим на примере нашей систем (2). Перепишем уравнения системы в интегральную форму:

$\left|\begin{array}{c}\frac{\partial x}{\partial t}\\ \\ \frac{\partial y}{\partial t}\\ \\ z\end{array}\right|\mathrm{<mo>=</mo>}\left|\begin{array}{c}x\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo>}\\ \\ y\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo>}\\ \\ 0\\ \end{array}\right|+{\displaystyle {\int }_0^t}\left|\begin{array}{c}\left({\eta }^T-\varphi z\right)y-{\eta }^{T}x\mathrm{,}\\ \\ \mathrm{<mo\; stretchy="false">(</mo>}D-z\mathrm{<mo\; stretchy="false">)</mo>}x-y\\ \\ K\mathrm{<mo\; stretchy="false">(</mo>}t-\tau \mathrm{<mo\; stretchy="false">)</mo>}Q\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">(</mo>}\tau \mathrm{<mo\; stretchy="false">)</mo>}y\mathrm{<mo\; stretchy="false">(</mo>}\tau \mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}\end{array}\right|dt\mathrm{.}$ (3)

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

Разностные схемы

Для начала в системе (2) введем следующие замены:

$x\mathrm{<mo>=</mo>}\left[\begin{array}{c}x\\ y\end{array}\right]\mathrm{,}a\mathrm{<mo>=</mo>}\left[\begin{array}{c}\mu \\ \varphi \\ D\\ \end{array}\right]-\text{векторпараметров}\mathrm{.}$

И перепишем систему в векторном виде:

$\begin{array}{ll}\frac{dx}{dt}\mathrm{<mo>=</mo>}f\left(x\mathrm{,}z\mathrm{,}a\right)\mathrm{,}\hfill & \hfill \\ z\mathrm{<mo>=</mo>}{\displaystyle {\int }_0^t}K\mathrm{<mo\; stretchy="false">(</mo>}t-\tau \mathrm{<mo\; stretchy="false">)</mo>}Q\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">(</mo>}\tau \mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}d\tau \mathrm{.}\hfill & \hfill \end{array}$ (4)

Вводим временную сетку $T$ с шагом $\Delta t$, ведем расчет для $t_k$ равно отстоящих. Через $x\mathrm{<mo\; stretchy="false">[</mo>}k\mathrm{<mo\; stretchy="false">]</mo>}$ и $z\mathrm{<mo\; stretchy="false">[</mo>}k\mathrm{<mo\; stretchy="false">]</mo>}$ обозначаем значения функции в эти самые моменты времени:

$\begin{array}{ll}{t}_k\mathrm{<mo>=</mo>}k\Delta t\mathrm{,}\hfill & \hfill \\ x\mathrm{<mo\; stretchy="false">[</mo>}k\mathrm{<mo\; stretchy="false">]</mo><mo>=</mo>}x\mathrm{<mo\; stretchy="false">(</mo>}{t}_k\mathrm{<mo\; stretchy="false">)</mo>,}\hfill & \hfill \\ z\mathrm{<mo\; stretchy="false">[</mo>}k\mathrm{<mo\; stretchy="false">]</mo><mo>=</mo>}z\mathrm{<mo\; stretchy="false">(</mo>}{t}_k\mathrm{<mo\; stretchy="false">)</mo>.}\hfill & \hfill \end{array}$

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

В качестве разностной схемы для дифференциальной части возьмем метод трапеции [7-9].

$\begin{array}{l}x\mathrm{<mo\; stretchy="false">[</mo>}k+\mathrm{1<mo\; stretchy="false">]</mo><mo>=</mo>}x\mathrm{<mo\; stretchy="false">[</mo>}k\mathrm{<mo\; stretchy="false">]</mo>}+\frac{\Delta t}{2}{m}_1+\frac{\Delta t}{2}{m}_2\mathrm{,}\hfill \\ m_1\mathrm{<mo>=</mo>}f\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">[</mo>}k\mathrm{<mo\; stretchy="false">]</mo>,}z\mathrm{<mo\; stretchy="false">[</mo>}k\mathrm{<mo\; stretchy="false">]</mo>,}a\mathrm{<mo\; stretchy="false">)</mo>,}\hfill \\ m_2\mathrm{<mo>=</mo>}f\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">[</mo>}k\mathrm{<mo\; stretchy="false">]</mo>}+\frac{\Delta t}{2}{m}_1+\frac{\Delta t}{2}{m}_2\mathrm{,}z\mathrm{<mo\; stretchy="false">[</mo>}k+\mathrm{1<mo\; stretchy="false">]</mo>,}a\mathrm{<mo\; stretchy="false">)</mo>.}\hfill \end{array}$ (5)

В расчетных целях его удобнее записать в следующем виде:

$x\mathrm{<mo\; stretchy="false">[</mo>}k+\mathrm{1<mo\; stretchy="false">]</mo><mo>=</mo>}x\mathrm{<mo\; stretchy="false">[</mo>}k\mathrm{<mo\; stretchy="false">]</mo>}+\frac{\Delta t}{2}\left(f\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">[</mo>}k\mathrm{<mo\; stretchy="false">]</mo>,}z\mathrm{<mo\; stretchy="false">[</mo>}k\mathrm{<mo\; stretchy="false">]</mo>,}a\mathrm{<mo\; stretchy="false">)</mo>}+f\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">[</mo>}k+\mathrm{1<mo\; stretchy="false">]</mo>,}z\mathrm{<mo\; stretchy="false">[</mo>}k+\mathrm{1<mo\; stretchy="false">]</mo>,}a\mathrm{<mo\; stretchy="false">)</mo>}\right)\mathrm{.}$ (6)

В качестве квадратурной формы для интегрального члена $–$ квадратурная формула трапеции [3, 11]:

$z\mathrm{<mo\; stretchy="false">[</mo>}k+\mathrm{1<mo\; stretchy="false">]</mo><mo>=</mo>}\frac{h}{2}\widehat{K}\mathrm{<mo\; stretchy="false">[</mo>}k+\mathrm{1<mo\; stretchy="false">]</mo>}Q\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">[</mo>0<mo\; stretchy="false">]</mo><mo\; stretchy="false">)</mo>}+h{\displaystyle \sum _{i\mathrm{<mo>=</mo>1}}^k}\widehat{K}\mathrm{<mo\; stretchy="false">[</mo>}\Delta t\mathrm{<mo\; stretchy="false">(</mo>}k-i\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">]</mo>}Q\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">[</mo>}i\mathrm{<mo\; stretchy="false">]</mo><mo\; stretchy="false">)</mo>}+\frac{h}{2}\widehat{K}\mathrm{<mo\; stretchy="false">[</mo>0<mo\; stretchy="false">]</mo>}Q\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">[</mo>}k+\mathrm{1<mo\; stretchy="false">]</mo><mo\; stretchy="false">)</mo>.}$ (7)

Расчет состояния системы в $k+1$ момент времени сводиться к тому, что из данной системы уравнений мы должны найти чему равняется $x\mathrm{<mo\; stretchy="false">[</mo>}k+\mathrm{1<mo\; stretchy="false">]</mo>}$ и $z\mathrm{<mo\; stretchy="false">[</mo>}k+\mathrm{1<mo\; stretchy="false">]</mo>}$. В тоже самое время хорошо видно, что переменные $x\mathrm{<mo\; stretchy="false">[</mo>}k+\mathrm{1<mo\; stretchy="false">]</mo>}$ и $z\mathrm{<mo\; stretchy="false">[</mo>}k+\mathrm{1<mo\; stretchy="false">]</mo>}$ фигурируют в правых частях равенства. Таким образом имеет дело не с явными формулами рекуррентного типа, которые позволяют легко и просто просчитывать значения, а с неявными уравнениями решение данной системы уравнений, размерность которой определяется размерностью фазового пространства данной системы и позволит рассчитать состояние системы в следующий шаг. Конечно расчет по явным схемам имеет большие преимущества с точки зрения быстродействия, однако хорошо известно, что неявные схемы имеют меньшие ограничения по выбору шага, с точки зрения устойчивости. Поскольку в изучаемых моделях можно ожидать появление хаотических режимов, очень чувствительных к расчетным ошибкам предпочтение было отдано неявным схемам. Поэтому полученную систему нелинейных уравнений относительно $x\mathrm{<mo\; stretchy="false">[</mo>}k+\mathrm{1<mo\; stretchy="false">]</mo>,}y\mathrm{<mo\; stretchy="false">[</mo>}k+\mathrm{1<mo\; stretchy="false">]</mo>}$ будем решать с помощью метода Ньютона [6, 10].

Метод Ньютона для систем:

$X\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo><mo>=</mo>}X\mathrm{<mo\; stretchy="false">[</mo>}n\mathrm{<mo\; stretchy="false">]</mo>}-W^{-1}\mathrm{<mo\; stretchy="false">(</mo>}X\mathrm{<mo\; stretchy="false">)</mo>}F\mathrm{<mo\; stretchy="false">(</mo>}X\mathrm{<mo\; stretchy="false">)</mo>,}$ (8)

где $W\mathrm{<mo\; stretchy="false">(</mo>}X\mathrm{<mo\; stretchy="false">)</mo>}$ $–$ Якобиан матрицы $F\mathrm{<mo\; stretchy="false">(</mo>}X\mathrm{<mo\; stretchy="false">)</mo>}$.

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

$\begin{array}{l}x\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo><mo>=</mo>}x\mathrm{<mo\; stretchy="false">[</mo>}n\mathrm{<mo\; stretchy="false">]</mo>}+\hfill \\ +\frac{h}{2}\mathrm{<mo\; stretchy="false">(</mo>}-\eta x\mathrm{<mo\; stretchy="false">[</mo>}n\mathrm{<mo\; stretchy="false">]</mo>}+\mathrm{<mo\; stretchy="false">(</mo>}\eta -\varphi z\mathrm{<mo\; stretchy="false">[</mo>}n\mathrm{<mo\; stretchy="false">]</mo><mo\; stretchy="false">)</mo>}y\mathrm{<mo\; stretchy="false">[</mo>}n\mathrm{<mo\; stretchy="false">]</mo>}\eta x\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo>}+\mathrm{<mo\; stretchy="false">(</mo>}\eta -\varphi z\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo><mo\; stretchy="false">)</mo>}y\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo><mo\; stretchy="false">)</mo>,}\hfill \\ y\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo><mo>=</mo>}y\mathrm{<mo\; stretchy="false">[</mo>}n\mathrm{<mo\; stretchy="false">]</mo>}+\frac{h}{2}\mathrm{<mo\; stretchy="false">(</mo><mo\; stretchy="false">(</mo>}D-z\mathrm{<mo\; stretchy="false">[</mo>}n\mathrm{<mo\; stretchy="false">]</mo><mo\; stretchy="false">)</mo>}x\mathrm{<mo\; stretchy="false">[</mo>}n\mathrm{<mo\; stretchy="false">]</mo>}-y\mathrm{<mo\; stretchy="false">[</mo>}n\mathrm{<mo\; stretchy="false">]</mo>}+D-z\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo><mo\; stretchy="false">)</mo>}x\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo>}-y\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo><mo\; stretchy="false">)</mo>.}\hfill \end{array}$ (9)

Теперь в каждое уравнение вместо $z\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo>}$ подставим

$z\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo><mo>=</mo>}\frac{h}{2}K\mathrm{<mo\; stretchy="false">[</mo>0<mo\; stretchy="false">]</mo><mo\; stretchy="false">(</mo>}{q}_{11}{x}^2\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo>}+2q_{12}x\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo>}y\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo>}+q_{22}{y}^2\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo><mo\; stretchy="false">)</mo>}+L\mathrm{<mo\; stretchy="false">[</mo>}n\mathrm{<mo\; stretchy="false">]</mo>,}$

где $q_{11}\mathrm{,}{q}_{22}$ $–$ коэффициенты при функции подавлении энергией $q_{12}$ $–$ коэффициент при функции подавлении спиральностью, $K\mathrm{<mo\; stretchy="false">[</mo>0<mo\; stretchy="false">]</mo>}$ $–$ ядро функционала подавления.

$L\mathrm{<mo\; stretchy="false">[</mo>}n\mathrm{<mo\; stretchy="false">]</mo><mo>=</mo>}K\mathrm{<mo\; stretchy="false">[</mo>}n\mathrm{<mo\; stretchy="false">]</mo>}Q\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">[</mo>0<mo\; stretchy="false">]</mo>,}y\mathrm{<mo\; stretchy="false">[</mo>0<mo\; stretchy="false">]</mo><mo\; stretchy="false">)</mo>}+{\displaystyle \sum _{k\mathrm{<mo>=</mo>1}}^n}K\mathrm{<mo\; stretchy="false">[</mo>}n+1-k\mathrm{<mo\; stretchy="false">]</mo>}Q\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">[</mo>}k\mathrm{<mo\; stretchy="false">]</mo>}y\mathrm{<mo\; stretchy="false">[</mo>}k\mathrm{<mo\; stretchy="false">]</mo><mo\; stretchy="false">)</mo>,}$

$Q\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}y\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{q}_{11}{x}^2\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo>}+2q_{12}x\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo>}y\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo>}+q_{22}{y}^2\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo>}$ $–$ квадратичная форма функции подавления(общий вид).

$\begin{array}{l}x\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo><mo>=</mo>}x\mathrm{<mo\; stretchy="false">[</mo>}n\mathrm{<mo\; stretchy="false">]</mo>}+\frac{h}{2}\mathrm{<mo\; stretchy="false">(</mo>}-\eta x\mathrm{<mo\; stretchy="false">[</mo>}n\mathrm{<mo\; stretchy="false">]</mo>}+\mathrm{<mo\; stretchy="false">(</mo>}\eta -\varphi z\mathrm{<mo\; stretchy="false">[</mo>}n\mathrm{<mo\; stretchy="false">]</mo><mo\; stretchy="false">)</mo>}y\mathrm{<mo\; stretchy="false">[</mo>}n\mathrm{<mo\; stretchy="false">]</mo>}-\eta x\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo>}+\hfill \\ +\mathrm{<mo\; stretchy="false">(</mo>}\eta -\varphi \mathrm{<mo\; stretchy="false">)</mo>}\frac{h}{2}K\mathrm{<mo\; stretchy="false">[</mo>0<mo\; stretchy="false">]</mo><mo\; stretchy="false">(</mo>}{q}_{11}{x}^2\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo>}+2q_{12}x\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo>}y\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo>}+q_{22}{y}^2\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo><mo\; stretchy="false">)</mo>}+L\mathrm{<mo\; stretchy="false">[</mo>}n\mathrm{<mo\; stretchy="false">]</mo><mo\; stretchy="false">)</mo>}y\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo><mo\; stretchy="false">)</mo>,}\hfill \\ y\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo><mo>=</mo>}y\mathrm{<mo\; stretchy="false">[</mo>}n\mathrm{<mo\; stretchy="false">]</mo>}+\frac{h}{2}\mathrm{<mo\; stretchy="false">(</mo><mo\; stretchy="false">(</mo>}D-z\mathrm{<mo\; stretchy="false">[</mo>}n\mathrm{<mo\; stretchy="false">]</mo><mo\; stretchy="false">)</mo>}x\mathrm{<mo\; stretchy="false">[</mo>}n\mathrm{<mo\; stretchy="false">]</mo>}-y\mathrm{<mo\; stretchy="false">[</mo>}n\mathrm{<mo\; stretchy="false">]</mo>}+\hfill \\ +\mathrm{<mo\; stretchy="false">(</mo>}D-\frac{h}{2}K\mathrm{<mo\; stretchy="false">[</mo>0<mo\; stretchy="false">]</mo><mo\; stretchy="false">(</mo>}{q}_{11}{x}^2\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo>}+2q_{12}x\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo>}y\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo>}+q_{22}{y}^2\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo><mo\; stretchy="false">)</mo>}+\hfill \\ +L\mathrm{<mo\; stretchy="false">[</mo>}n\mathrm{<mo\; stretchy="false">]</mo><mo\; stretchy="false">)</mo>}x\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo>}-y\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo><mo\; stretchy="false">)</mo>.}\hfill \end{array}$ (10)

Мы получили систему нелинейных алгебраических уравнений неизвестными в которой являются $x\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo>}$ и $y\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo>}$.

Перенесем левую часть уравнения в правую.

$\begin{array}{l}{F}_1\mathrm{<mo>=</mo>}-x\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo>}+x\mathrm{<mo\; stretchy="false">[</mo>}n\mathrm{<mo\; stretchy="false">]</mo>}+\frac{h}{2}\mathrm{<mo\; stretchy="false">(</mo>}-\eta x\mathrm{<mo\; stretchy="false">[</mo>}n\mathrm{<mo\; stretchy="false">]</mo>}+\mathrm{<mo\; stretchy="false">(</mo>}\eta -\varphi z\mathrm{<mo\; stretchy="false">[</mo>}n\mathrm{<mo\; stretchy="false">]</mo><mo\; stretchy="false">)</mo>}y\mathrm{<mo\; stretchy="false">[</mo>}n\mathrm{<mo\; stretchy="false">]</mo>}-\eta x\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo>}+\hfill \\ +\mathrm{<mo\; stretchy="false">(</mo>}\eta -\varphi \mathrm{<mo\; stretchy="false">)</mo>}\frac{h}{2}K\mathrm{<mo\; stretchy="false">[</mo>0<mo\; stretchy="false">]</mo><mo\; stretchy="false">(</mo>}{q}_{11}{x}^2\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo>}+2q^{12}x\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo>}y\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo>}+q_{22}{y}^2\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo><mo\; stretchy="false">)</mo>}+L\mathrm{<mo\; stretchy="false">[</mo>}n\mathrm{<mo\; stretchy="false">]</mo><mo\; stretchy="false">)</mo>}y\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo><mo\; stretchy="false">)</mo>,}\hfill \\ F_2\mathrm{<mo>=</mo>}-y\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo>}+y\mathrm{<mo\; stretchy="false">[</mo>}n\mathrm{<mo\; stretchy="false">]</mo>}+\frac{h}{2}\mathrm{<mo\; stretchy="false">(</mo><mo\; stretchy="false">(</mo>}D-z\mathrm{<mo\; stretchy="false">[</mo>}n\mathrm{<mo\; stretchy="false">]</mo><mo\; stretchy="false">)</mo>}x\mathrm{<mo\; stretchy="false">[</mo>}n\mathrm{<mo\; stretchy="false">]</mo>}-y\mathrm{<mo\; stretchy="false">[</mo>}n\mathrm{<mo\; stretchy="false">]</mo>}+\hfill \\ \mathrm{<mo\; stretchy="false">(</mo>}D-\frac{h}{2}K\mathrm{<mo\; stretchy="false">[</mo>0<mo\; stretchy="false">]</mo><mo\; stretchy="false">(</mo>}{q}_{11}{x}^2\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo>}+2q_{12}x\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo>}y\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo>}+q_{22}{y}^2\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo><mo\; stretchy="false">)</mo>}+\hfill \\ +L\mathrm{<mo\; stretchy="false">[</mo>}n\mathrm{<mo\; stretchy="false">]</mo><mo\; stretchy="false">)</mo>}x\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo>}-y\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo><mo\; stretchy="false">)</mo>.}\hfill \end{array}$ (11)

Теперь вычислим Якобианы системы по следующей формуле:

$W\mathrm{<mo>=</mo>}\left|\begin{array}{cc}\frac{\partial F_1}{\partial x\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo>}}& \frac{\partial F_1}{\partial y\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo>}}\\ & \\ \frac{\partial F_2}{\partial x\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo>}}& \frac{\partial F_2}{\partial y\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo>}}\\ & \end{array}\right|\mathrm{.}$ (12)

Получим:

$W_{11}\mathrm{<mo>=</mo>}-\frac{1}{2}\varphi h^{2}K\mathrm{<mo\; stretchy="false">[</mo>0<mo\; stretchy="false">]</mo>}{q}_{11}x\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo>}y\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo>}-\frac{1}{2}\varphi h^{2}K\mathrm{<mo\; stretchy="false">[</mo>0<mo\; stretchy="false">]</mo>}{q}_{11}{y}^2\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo>}-1-\frac{1}{2}\eta h\mathrm{,}$

$\begin{array}{l}{W}_{12}\mathrm{<mo>=</mo>}-\frac{1}{4}\varphi h^{2}K\mathrm{<mo\; stretchy="false">[</mo>0<mo\; stretchy="false">]</mo>}{q}_{11}x{\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo>}}^2-\varphi h^{2}K\mathrm{<mo\; stretchy="false">[</mo>0<mo\; stretchy="false">]</mo>}q12x\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo>}y\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo>}-\hfill \\ -\frac{3}{4}\varphi h^{2}K\mathrm{<mo\; stretchy="false">[</mo>0<mo\; stretchy="false">]</mo>}{q}_{22}{y}^2\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo>}+\frac{1}{2}h\left(\eta -\varphi L\mathrm{<mo\; stretchy="false">[</mo>}n\mathrm{<mo\; stretchy="false">]</mo>}\right)\mathrm{,}\hfill \end{array}$

$\begin{array}{l}{W}_{12}\mathrm{<mo>=</mo>}{W}_{21}\mathrm{<mo>=</mo>}-\frac{3}{4}{h}^{2}K\mathrm{<mo\; stretchy="false">[</mo>0<mo\; stretchy="false">]</mo>}{q}_{11}{x}^2\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo>}-h^{2}K\mathrm{<mo\; stretchy="false">[</mo>0<mo\; stretchy="false">]</mo>}{q}_{12}x\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo>}y\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo>}-\hfill \\ -\frac{1}{4}{h}^{2}K\mathrm{<mo\; stretchy="false">[</mo>0<mo\; stretchy="false">]</mo>}{q}_{22}{y}^2\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo>}+\frac{1}{2}h\left(D-L\mathrm{<mo\; stretchy="false">[</mo>}n\mathrm{<mo\; stretchy="false">]</mo>}\right)\mathrm{,}\hfill \end{array}$

$W_{22}\mathrm{<mo>=</mo>}-\frac{1}{2}{h}^{2}K\mathrm{<mo\; stretchy="false">[</mo>0<mo\; stretchy="false">]</mo>}{q}_{12}{x}^2\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo>}-\frac{1}{2}{h}^{2}K\mathrm{<mo\; stretchy="false">[</mo>0<mo\; stretchy="false">]</mo>}{q}_{22}x\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo>}y\mathrm{<mo\; stretchy="false">[</mo>}n+\mathrm{1<mo\; stretchy="false">]</mo>}-1-\frac{1}{2}h\mathrm{.}$

Алгоритм решения:

1. Зададим начальные условия $x\mathrm{<mo\; stretchy="false">[</mo>0<mo\; stretchy="false">]</mo><mo>=</mo>}y\mathrm{<mo\; stretchy="false">[</mo>0<mo\; stretchy="false">]</mo>}$.
2. На текущий момент времени $t$ (шаг $i$ ) пусть будут известны значения $x\mathrm{<mo\; stretchy="false">[</mo>}i\mathrm{<mo\; stretchy="false">]</mo>,}y\mathrm{<mo\; stretchy="false">[</mo>}i\mathrm{<mo\; stretchy="false">]</mo>}$.
3. Вычисляем функции $L\mathrm{<mo\; stretchy="false">[</mo>}n\mathrm{<mo\; stretchy="false">]</mo>}$ и $z\mathrm{<mo\; stretchy="false">[</mo>}n\mathrm{<mo\; stretchy="false">]</mo>}$.
4. Вычисляем Якобиан системы $W$.
5. находим обратную матрицу к нашему Якобиану $W^{-1}$.
6. Считаем значение функции $F_1$ и $F_2$ (11) в текущий момент времени.
7. По формуле метода Ньютона (8) считаем $i+1$ значения для $x$ и $y$.
8. Увеличим временной индекс $t$ на значение $h$ (шаг $i$ на $1$ ) и переходим на шаг алгоритма 2.

Сопоставление с динамикой известной системы

Интегральный член является признаком эредитраности модели (2). В то же время его можно исключить для некоторых типов ядер с экспоненциальной асимптотикой за счет расширения размерности фазового пространства модели. Точнее говоря, если ядро является решением линейного дифференциального уравнения с постоянными коэффициентами, то система (2) равносильна некоторой дифференциальной системе, с начальными условиями на дополнительные фазовые переменные. А именно справедлива теорема [12]:

Теорема Если ядро $K\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}$ является решением дифференциального уравнения

$\begin{array}{l}{a}_0{K}^{\mathrm{<mo\; stretchy="false">(</mo>}n\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}+a_1{K}^{\mathrm{<mo\; stretchy="false">(</mo>}n-\mathrm{1<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}+\dots +a_{n-1}{K}^{\prime }\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}+a_{n}K\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>0}\hfill \end{array}$ (13)

с постоянными коэффициентами $a_i$, то интегральное равенство

$z\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\displaystyle {\int }_0^t}K\mathrm{<mo\; stretchy="false">(</mo>}t-\tau \mathrm{<mo\; stretchy="false">)</mo>}Q\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">(</mo>}\tau \mathrm{<mo\; stretchy="false">)</mo>,}y\mathrm{<mo\; stretchy="false">(</mo>}\tau \mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}d\tau \mathrm{.}$

равносильно следующей задаче Коши для функции $z\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}$:

$\begin{array}{l}{a}_0\frac{d^{n}z}{d{t}^n}+a_1\frac{d^{n-1}z}{d{t}^{n-1}}+\dots +a_{n}z\mathrm{<mo>=</mo>}{\displaystyle \sum _{l\mathrm{<mo>=</mo>1}}^n}{a}_{n-l}{\displaystyle \sum _{s\mathrm{<mo>=</mo>0}}^{l-1}}{K}^{\mathrm{<mo\; stretchy="false">(</mo>}s\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo>}\frac{d^{l-1-s}}{d{t}^{l-1-s}}Q\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>,}y\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>,}\hfill \\ z\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo><mo>=</mo>0,}\hfill \\ \frac{d^{l}z\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}}{d{t}^l}{\mathrm{<mo>|</mo>}}_{t\mathrm{<mo>=</mo>0}}\mathrm{<mo>=</mo>}{\displaystyle \sum _{s\mathrm{<mo>=</mo>0}}^{l-1}}{K}^{\mathrm{<mo\; stretchy="false">(</mo>}s\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo>}\frac{d^{l-1-s}}{d{t}^{l-1-s}}Q\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>,}y\mathrm{<mo\; stretchy="false">(</mo>}t{\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo><mo>|</mo>}}_{t\mathrm{<mo>=</mo>0}}\mathrm{,}l\mathrm{<mo>=</mo>1,}\dots \mathrm{,}n-1.\hfill \end{array}$ (14)

Рассмотрим процесс верификации на нашей системе (2), в возьмем в качестве параметров следующие $\varphi \mathrm{<mo>=</mo>0}$, $Q\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">(</mo>}\tau \mathrm{<mo\; stretchy="false">)</mo>}y\mathrm{<mo\; stretchy="false">(</mo>}\tau \mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo><mo>=</mo>}xy$, $K\mathrm{<mo\; stretchy="false">(</mo>}t-\tau \mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{e}^{-bt}$, ${\eta }^T\mathrm{<mo>=</mo>10}$ и $b\mathrm{<mo>=</mo>8<mo>/</mo>3}$. Воспользуемся условием теоремы об исключении интегрально члена из системы (2), получим систему вида:

Получим систему вида:

$\begin{array}{ll}\frac{dx}{dt}\mathrm{<mo>=</mo>}\left(y-x\right){\eta }^T\mathrm{,}\hfill & \hfill \\ \frac{dy}{dt}\mathrm{<mo>=</mo><mo\; stretchy="false">(</mo>}D-z\mathrm{<mo\; stretchy="false">)</mo>}x-y\mathrm{,}\hfill & \hfill \\ \frac{dz}{dt}\mathrm{<mo>=</mo>}xy-bz\mathrm{,}\hfill & \hfill \\ z\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo><mo>=</mo>0.}\hfill & \hfill \end{array}$ (15)

Эта система $\alpha \omega$ -динамо, соответствующая системе Лоренца. Динамика такой системы прекрасно известна[13].

Будем варьировать параметр D и будем получать различные режимы системы.

Для D = 20 получим асимптотически устойчивый режим рис. 1,2.

Рис. 1. Фазовая координаты: (a) X при D = 20; (b) Y при D = 20.

[Figure 1. Phase coordinates: (a) X at D = 20; (b) Y for D = 20.]

Рис. 2. Фазовый портрет системы при D = 20.

[Figure 2. Phase portrait of the system at D = 20.]

D = 220 получим периодический режим рис. 3,4.

Рис. 3. Фазовая координаты: (a) X при D = 220; (b) Y при D = 220.

[Figure 3. Phase coordinates: (a) X at D = 220; (b) Y for D = 220.]

Рис. 4. Фазовый портрет системы при D = 220.

[Figure 4. Phase portrait of the system at D = 220.]

Если же D = 28 получим хаотический режим рис. 5,6.

Рис. 5. Фазовая координаты: (a) X при D = 28; (b) Y при D = 28.

[Figure 5. Phase coordinates: (a) X at D = 28; (b) Y for D = 28.]

Рис. 6. Фазовый портрет системы при D = 28.

[Figure 6. Phase portrait of the system at D = 28.]

Порядок точности численной схемы

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

Погрешность в результате вычислений может возникнуть по следующим причинам:

1. погрешность при моделировании, любое разностное уравнение не является абсолютным эквивалентом дифференциальному уравнению, это основной источник погрешности.
2. округление чисел при вычислении.
3. погрешность в значениях правой части $f\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}y\mathrm{<mo\; stretchy="false">)</mo>}$, погрешность вызвана тем фактам, что рассматривается некоторое приближение функции к правой части дифференциального уравнения. Так же в ходи вычисления на ЭВМ функция $f\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}y\mathrm{<mo\; stretchy="false">)</mo>}$ может быть приближена другими функциями, что в свою очередь вносит дополнительную погрешность в решение уравнения.
4. $y_{i+1}$ определяется и уравнения которое эквивалентно исходному, но не может быть разрешимо в явном виде.

Исследование порядка точности по правилу Рунге заключается в следующем: берется решение на сетки $h$ и сравнивается с решением на сетке с шагом в $2$ раза меньше т.е. $h\mathrm{<mo>/</mo>2}$. И в дальнейшем рассмотрении разностей погрешностей для этих двух вычислений по формуле [14, 15]:

$ϵ\mathrm{<mo>=</mo>}\frac{\mathrm{<mo>|</mo>}{y}_h-y_{h\mathrm{<mo>/</mo>2}}\mathrm{<mo>|</mo>}}{2^p-1}$ (16)

Формула Рунге справедлива для всех вычислительных процессов, для которых выполняется степенной закон. Для определения порядка метода $p$ необходимо проведение априорной оценки погрешности, что не всегда легко осуществить [16]

Английский математик Эйткен предложил способ оценки погрешности для случая, когда порядок $р$ метода неизвестен . Более того, алгоритм Эйткена позволяет опытным путем определить и порядок метода. Для этого необходимо третий раз вычислить значение величины $y$ с шагом $h\mathrm{<mo>/</mo>4}$ [17].

Для удобства введем переменную $k$, которая в данном случаем равна $\frac{1}{2}$. Для вычисление порядка $p$ используем формулу Эйткена:

$k^p\mathrm{<mo>=</mo>}\frac{y_{kh}-y_{k^{2}h}}{y_h-y_{kh}}\mathrm{.}$ (17)

Из формулы (17) получаем следующее соотношение:

$p\mathrm{<mo>=</mo>}\frac{\log W}{\log k}\mathrm{,}$ (18)

где $W\mathrm{<mo>=</mo>}\frac{y_{kh}-y_{k^{2}h}}{y_h-y_{kh}}\mathrm{.}$.

Были проведены численные эксперименты при шагах $h$, $h\mathrm{<mo>/</mo>2}$, $h\mathrm{<mo>/</mo>4}$, $h\mathrm{<mo>/</mo>8}$, $h\mathrm{<mo>/</mo>16}$, $h\mathrm{<mo>/</mo>32}$. Управляющие параметры исследуемой системы (2) были выбраны следующие:

${\eta }^T\mathrm{<mo>=</mo>10}D\mathrm{<mo>=</mo>15}b\mathrm{<mo>=</mo>}\frac{8}{3}K\mathrm{<mo\; stretchy="false">(</mo>}t-\tau \mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{e}^{-bt}Q\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}y\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}xy\mathrm{.}$

Результаты расчетов приведены в таблице 1.

Таблица 1

Расчет порядка точности по формуле Эйткена [Calculation of the order of accuracy using the Aitken formula]

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

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

Заключение

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

Исследования порядка точности предложенной численной схемы проводилось при помощи метода Эйткена. Полученный порядок равен 1.

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

Об авторах

Евгений Анатольевич Казаков

Автор, ответственный за переписку.
Email: Kazakov@ikir.ru
ORCID iD: 0000-0001-7235-4148

младший научный сотрудник лаборатории электромагнитных излучений

Список литературы

Дополнительные файлы

Доп. файлы

Действие

1. JATS XML

Скачать

2. Рис. 1. Фазовая координаты: (a) X при D = 20; (b) Y при D = 20.

Скачать (322KB)

Метаданные

3. Рис. 2. Фазовый портрет системы при D = 20.

Скачать (171KB)

Метаданные

4. Рис. 3. Фазовая координаты: (a) X при D = 220; (b) Y при D = 220.

Скачать (631KB)

Метаданные

5. Рис. 4. Фазовый портрет системы при D = 220.

Скачать (153KB)

Метаданные

6. Рис. 5. Фазовая координаты: (a) X при D = 28; (b) Y при D = 28.

Скачать (540KB)

Метаданные

7. Рис. 6. Фазовый портрет системы при D = 28.

Скачать (344KB)

Метаданные