CC BY
6
Вестник КРАУНЦ. Физ.-мат. науки. 2022. Т. 41. №4. C. 47-65. ISSN 2079-6641
УДК 519.642.2, 517.938 Научная статья
Математической моделирование в MATLAB циклов солнечной активности по данным роста-спада числа Вольфа
Д. А. Твёрдый, Р. И. Паровик
Институт космофизических исследований и распространения радиоволн ДВО РАН, 684034, c. Паратунка, ул.Мирная, д. 7, Россия E-mail: dimsolid95@gmail.com
В этой статье проводится математическое моделирование динамики солнечной активности. Исследуются данные наблюдений по средне-ежемесячному числу солнечных пятен, называемых числом Вольфа, в период за 24.5 года с мая 1996 года по октябрь 2022 года. Исходя из результатов подобного исследования данных по этому процессу, с применением уравнения Риккати дробного постоянного порядка, о том, что подъем и падение числа Вольфа со временем происходит по кривой, очень близкой к обобщенной логистической кривой, в данной статье также предлагается математическая модель, основанная на уравнении Риккати. Так как уравнение Риккати хорошо описывает процессы, которые подчиняются логистическому закону. Однако уравнение обобщается до интегро-дифференциального уравнения Риккати, введением дробной производной типа Герасимова-Капуто переменного порядка, а дробная производная с переменным порядком, позволяет получить более уточную математическую модель циклов числа Вольфа с насыщением, и позволяет учесть эффект переменной памяти. Все расчёты моделей, обработка данных и визуализации проводятся в программе FDRE 3.0 разработанной в пакете MATLAB. Параметры моделирования уточняются аппроксимацией известных исследуемых данных, при помощи регрессионного анализа. В результате модельные кривые и графики известных за 24.5 года наблюдаемых данных, показывают между собой хорошее соответствие. С помощью уточнённой математической модели делается прогноз на следующие 9 лет, который визуально хорошо согласуется с известными модельными результатами солнечной активности.
Ключевые слова: солнечная активность, число Вольфа, коэффициент детерминации, коэффициент корреляции, математическое моделирование, динамические процессы, эффект насыщения, эредитарность, уравнение Риккати, производная типа Герасимова-Капуто, переменный порядок дробной производной.
d DOI: 10.26117/2079-6641-2022-41-4-47-65
Поступила в редакцию: 19.11.2022 В окончательном варианте: 22.11.2022
Для цитирования. Твёрдый Д. А., Паровик Р. И. Математической моделирование в MATLAB циклов солнечной активности по данным роста-спада числа Вольфа // Вестник КРАУНЦ. Физ.-мат. науки. 2022. Т. 41. № 4. C. 47-65. DOI: 10.26117/2079-66412022-41-4-47-65
Контент публикуется на условиях лицензии Creative Commons Attribution 4.0 International (https://creativecommons.org/licenses/by/4.0/deed.ru)
© Твёрдый Д. А., Паровик Р. И., 2022
Финансирование. Исследования выполнены в рамках государственного задания ИКИР ДВО РАН по теме АААА-А21-121011290003-0
Введение
Исследования солнечно-земных связей — реакций внешних оболочек Земли, в том числе и биосферы, — на изменение солнечной активности (далее БЛ) активно проводились на протяжении всего 20 века Земля получает не только свет и тепло от Солнца, но и подвергается воздействию ультрафиолетового- и рентгеновского-излучений, солнечного ветра и т.д. Колебания мощности этих факторов с изменением уровня БЛ вызывают цепочку явлений в межпланетном пространстве особенно заметны внешние оболочки Земли, изменения в магнитосфере.
Была статистически установлена связь для роста числа аварий в электросетях СШ вблизи авроральной зоны в след и ростом уровня геомагнитной активности. Однако в годы минимальной активности вероятности аварий на опасных и безопасных участках практически уравниваются [1]. Предполагается также, что существует связь между уровнем БЛ и геомагнитной возмущенности и ходом ряда процессов в биосфере Земли, таких как: динамика численности животных, эпидемий, количество сердечно-сосудистых кризов и др.
Изучение солнечно-земных связей является не только фундаментальной научной проблемой, но и имеет большое прогностическое значение. Прогнозы состояния магнитосферы и других оболочек Земли крайне необходимы для решения практических задач в области космонавтики, радиосвязи, транспорта, метеорологии и климатологии, сельского хозяйства, биологии и медицины.
Наиболее сильным и заметным проявлением БЛ является мощная вспышка на Солнце, последствия которой могут сказываться на околоземном пространстве. В частности, ускоренные во вспышке частицы, вторгаясь в ионосферу и стратосферу полярных широт, вызывают длительное, на десятки часов, ухудшение качества коротковолновой радиосвязи, способствуют истощению озонового слоя и т.д. Отсюда и возникает интерес к наблюдению за активностью магнитосферы Солнца, с целью предсказания отписанных и иных возможных последствий для человечества.
Пятна на Солнце — явный признак его активности, образование которых связано с магнитным полем Солнца. После 17 лет наблюдений Генрих Швабе обнаружил, что количество солнечных пятен со временем меняется. Годы минимума солнечных пятен на поверхности Солнца могут быть нулевыми; в годы максимума солнечных пятен их количество измеряется десятками. Максимумы и минимумы чередуются в среднем каждые 11 лет (от 7 до 17 лет); однако могут быть более длительные циклы СА. А их количество называют - числом Вольфа, или относительным числом солнечных пятен, и является важным показателем для предсказания распространения радиоволн с отражением от ионосферы.
Нелокальное нелинейное уравнение с учётом эффекта памяти, и неявная численная схема решения
На процессы насыщения может влиять и наследственность (память или наследственность), что свидетельствует о причинно-следственной связи в динамике про-
цесса. Иначе говоря система хранит информацию о своей предыстории и с точки зрения математики может быть описана с помощью интегро-дифференциальных уравнений. Известный итальянский математик Вито Вольтера часть своих научных трудов [2] посвятил разработке понятия наследственности и применению его в различных областях науки, в частности, в проблемах экологии и физики. Где также отмечает что такие явления, как запаздывающие волны или усталость металлов были известны и ранее, что указано в работе Учайкина В.В. [3]. Исследование того как понятие наследственности можно обосновать математически закономерно продолжается во всём мире, в силу очевидно высокой практической значимости. А в России одним из наиболее значимых учёных по данному направлению является Учайкин Владимир Васильевич, который уже на протяжении более 25 лет, с 1995 года занимается разработкой численных методов решения уравнений в дробных производных.
Рассмотрим следующее эредитарное (наследственное) уравнение:
-г
к(г - а,г)х(а)аа - ь(г)и(г) = адг),г), (1)
0
где и(г) € С2 [0, Т] - неизвестная функция решения, г € [0, Т] - текущий момент времени, Т> 0 - общее время моделирования, ио - заданная константа, Ь(г) >0 -непрерывная функция, ^и(г),г) - нелинейная функция, удовлетворяющая условию Липшица.
Замечание 1. Полная эредитарность, по В. Вольтера [2] согласно принципам эредитарности определяется на интервале (—оо,г). Будем считать наследственность определенной на подинтервале (0,г), т. е. ограниченной.
Важной особенностью является то, что к понятию наследственности приводят разностные ядра интегро-дифференциальных уравнений - называемые функциями памяти. И при выборе в (1) разностного ядра - как степенной функций, позволит естественным образом перейти к хорошо развитому аппарату (точнее операторам) дробного исчисления [4], в частности к производным дробного порядка [5, 6, 7]. Что послужило теоретической основой для выведения численных аппроксимаций таких дробных операторов, и как следствие, схемы численного решения предлагаемой здесь математической модели. Тогда выберем функцию памяти следующим образом:
(г — а)—а(г)
к(г—а,г) = ^ 0<а(г)< 1,
I (1 — а(г))
где а(г) - функция, отвечающая за интенсивность исследуемого процесса.
Замечание 2. Зависимость от г приводит нас к явлению переменной памяти. А особенностью случая степенной функции является то, что процесс постепенно "забывает"о своей предыстории.
Тогда уравнение (1) можно записать в терминах производной типа Герасимова-Капуто (3), в виде задачи Коши для нелинейного дробного уравнения (1) с переменными коэффициентами:
аО(г)и(а)= ь(г)и(г) + ^и(г),г), и(0) = и0, (2)
где оператор дробного переменного порядка 0 < а("Ь) < 1 типа Герасимова-Капуто, действующий на функцию и("Ь) € С[0,Т] имеет вид:
aa(t) f л 1
30t( )u(a) =
Г (1 - a(t))J
Х(а) da, (3)
(t - a)a(t)
где Г(.)-Гамма-функция Эйлера, производная x(t) = -f.
Замечание 3. Операторы вида (3), обобщённые на случай производной дробного переменного порядка иногда называют просто VO "переменного поряд-ка"согласно [8, 9]. Дробное дифференцирование и дробное интегрирование являются обобщениями дифференцирования и интегрирования целого порядка и включают n-ые производные и n-кратные интегралы, где n е N, как частные случаи.
Ввиду нелинейности задачи Коши (2) ее решение будем искать с помощью численного метода конечно-разностных схем [10, 11, 12]. Для чего задачу (9) в дифференциальной постановке, подвергнем дискретизации. Разделим на N равноудалённых узлов сетки отрезок [0,T] с шагом дискретизации h = T/N. Тогда u(t) е C2[0,T] - неизвестная функция решения, на конечной сетке, перейдёт в свой сеточный аналог u(ti) или ui, где i = 1,...,N. Таким же образом и непрерывная гладкая функция 0 < a(t) < 1 перейдёт в a(ti) или ai, что верно и для a(t),b(t),c(t).
Аппроксимацию производной дробного переменного порядка типа Герасимова-Капуто (3) запишем, согласно [13] для t е [ti, ti+1 ] следующим образом:
aa(t) f л 1
Э0.( )u(a) =
'<* ^"Г (1 - a(t))
* x(a)da iv1 i ,
0 й-О^М * ... * A^wi(ui-j-ui-j-
Г (2 - Oi)'
Подставляя (4) в (2), получим дискретный аналог:
i-1
Ai = mT^T, W = (j +1)1-"1 - j1-"1, (4)
Ai Ц wj (ui-j -ui-j-1) = biui + fi, uo = C, (5)
j=o
где С — заданная константа, Ъ = Ъ(") ^ = ^и.,".).
Лемма 1. Аппроксимация ЗЮ^х^а) оператора типа Герасимова-Капуто 30|(")х(а) вида (3) удовлетворяет следующей оценке:
sfW)- 30at(t)u(a)
< Ст2-а, a= max(a(ti)), (6)
где С - независимая от шага дискретизации константа т.
Лемма 2. Дискретная задача Коши (5) аппроксимирует исходную дифференциальную задачу (2) с порядком:
тах|и(^)-и^|= 0(т2-"). (7)
t
o
Запишем (5) в виде нелокальной неявной конечно-разностной схемы (IFDS):
W = (j + 1)1—а - j1-а,
i = 1,...,N, uo = C. (8)
Теорема 1. IFDS (8) безусловно сходится с порядком 2 — сс.
Определение 1. Нелокальная IFDS (8) называется условно устойчивой, если справедлива оценка |Ui — Wi| < C|U0 — W0|, Vi, где C > 0 константа, независящая от шага т, где Ui,Wi - два различных решения схемы (8) в матричной форме: MUi = Fi.
Теорема 2. Нелокальная IFDS (8) безусловно устойчива.
Все доказательства лемм, теорем и оценки полученные в данном исследовании, приведены и опубликованы в статье [14] высоко рейтингового журнала Fractal and Fractional квартиля Q1 с импакт фактором IF = 3.396 за последние 5 лет.
Эредитарная модель SA на основе уравнения Риккати
Следует заметить, что проводились исследования близкие к нашей работе, например в работе Бураева А.В. [15], была исследована динамика SA в период 19982010 гг., и установлена ее связь с селевыми потоками в Кабардино-Балкарской Республике. Так же в этой работе на основе данных о динамике SA показано, что её подъем и падение происходят по кривой, очень близкой к обобщенной логистической кривой [16, 17, 18], и этот процесс нелинеен и фрактален [19]. Однако в этой работе модельное уравнение Риккати с производной дробного постоянного порядка было с постоянными коэффициентами.
В нашем же исследовании, исходя из физических предположений нелинейности и фрактальности [19] процесса, будет строится математическая модель, основанная на задаче Коши для дробного Уравнения Риккати, так как уравнение Риккати хорошо описывает процессы, которые подчиняются логистическому закону [16, 20, 21] и процессы с насыщением, а произвольный порядок дробной производной дает широкий диапазон для уточнения математической модели с насыщением и учитывает эффект переменной памяти динамической системы.
Многочисленные теоретические и практические исследования в мире показывают, что уравнение Риккати представляет заметный интерес, особенно в случае с (дробной) производной. Существует множество приложений уравнения Риккати [22] часто находящих свое применение во многих областях науки. Но для данного исповедования интерес представляет приложения в логистических моделях, целью которых является определение времени насыщения (выхода на плато) [23].
1-1
Ai И W К- -U-j-1) -biUi = fi, j=0
Ai =
T-«t
Г (2 - ai)'
Переход к дробному уравнению Риккати от обобщённого нелинейного случая (2), можно осуществить положив f(u(t),t) = —a(t)u(t)2 + c(t), причем a(t),c(t) G C[0,T] - заданная функция. Тогда задачу для (1) можно поставить следующую задачу Коши:
9j(t)u(ff) + a(t)x2(t) — b(t)u(t) — c(t) = 0, u(0) = U0, (9)
где u0 - заданная константа, 0 < a(t) < a, 0 < b(t) < b, 0 < c(t) < С - непрерывные функции, а b(t) > a(t), a,b,c - константы.
Определение 2. Уравнение в (9) будем называть дробным уравнением Риккати в переменными коэффициентами и дробной производной Герасимова-Капуто (3). А модель (9) будем называть: эредитарной моделью SA.
Дискретная форма (9) в виде IFDS, аналогично случаю обобщенного нелинейного уравнения (8) примет вид :
i—1
Ai И wi (xi—j — xi—j—1) + aix2 — biXi — ci = 0,
j=0
i = 1,...,N, x0 = C — const, (10)
В результате получаем систему нелинейных алгебраических уравнений, представленную схемой (10). Для ее решения воспользуемся методом MNM [12, 14, 24], представив (10) в виде итерационной функции:
i—1
F(ui) = Ai ^ w- (u—j — u—j—1) + aiu2 — biui — ci, j=0
с помощью которой составим итерационный процесс:
F (Um)
Um+1 = Um—JU^, (11)
подробнее об MNM для этой задачи в [14, p.16].
Известно, что обычный метод Ньютона (NM) требует вычисления матрицы Яко-би (то есть решения системы нелинейных алгебраических уравнений) на каждой итерации. Использование MNM позволяет вычислить якобиан только один раз в начальном приближении U0, но это приводит к снижению порядка сходимости метода (по сравнению с NM) до первого порядка.
Замечание 4. Mетод решения (10) - MNM локально устойчивый для IFDS и сходящийся с первым порядком, согласно [14].
Определение 3. IFDS (10) безусловно устойчива согласно теореме 2, а вычислительная точность при использовании MNM стремится к теоретической pAprior = 2 — max(a(t)), что соответствует условию Леммы 2, как показано в компьютерных экспериментах [14, р.17-23].
Определение 4. В качестве начального приближения U0 для MNM (11) можно взять последнее значение u^, полученное по схеме EFDS также предложенной и исследованной в [14, р.14-17], при выполнении условия сходимости EFDS.
Поскольку EFDS условно сходится с первым порядком, согласно Теореме 5 [14, р.14-17], а также [12].
Моделирование динамики 8Л
Известно мало работ в которых бы моделировались динамические процессы с учётом эффекта памяти, и проводилось сопоставление расчётных данных с экспериментальными или наблюдаемыми. Такие работы в основном и имею уклон в физическое (реальное) моделирование процесса [25, 26], но с попыткой предварительно построить некую теоретическую основу для проверки [26, 27].
Ещё меньше аналогичных исследований SA для процессов насыщением и эффектом памяти, что больше бы соответствовало тематике данного исследования. Редко проводится сравнение экспериментальных данных процессов с насыщением и расчетных данных по математической модели.
Рис. 1. Данные среднему значению числа Вольфа с шагом в 1 месяц : а) с Мая 1996 гг по Август 2021 гг; b) с Мая 1996 гг по Октябрь 2022 гг [Figure 1. Data to the average value of the Wolf number in 1-month steps: a) from May 1996 to August 2021; b) May 1996 to October 2022]
Поскольку изменение количества солнечных пятен является наиболее изученным видом SA, мы можем применить описанную (9) модель для аппроксимации данных по замерам числа Вольфа с мая 1996 гг по октябрь 2022 гг, т.е. 23, 24 и середину текущего 25-го циклов SA с шагом 1 месяц.
Для данных наблюдений произведём нормировку количества солнечных пятен на максимум; поэтому данные наблюдений приведены в относительных единицах. Результаты моделирования не нормировались.
Данные по процессам SA, используемые при моделировании получены из базы данных World Data Center для производства, сохранения и распространения числа солнечных пятен, с веб-сайта Королевской бельгийской обсерватории www.sidc.be, проекта Sunspot Index [28] по долгосрочным наблюдениям за Солнцем.
Предложенная ранее модель
Используемые параметры моделирования рассчитывались при различных вариациях значений параметров моделирования, из которых выбиралась наиболее подходящая комбинация по максимальным коэффициентам: корреляции Пирсона (Сог) [29] и детерминации (И2) [30, 31], с нормированными данными наблюдений. Определив таким образом неизвестные функции а("Ь), Ъ("Ь), е("Ь) и сс("Ь) на основе известных данных наблюдений, модель (9) можно скорректировать так:
30t(t)u(a) + a(t)x2(t) — b(t)u(t) — c(t) = 0, u(0) = u0,
a(t) = — 0.5 exp ( —| cos t
T
T
+ 0.25 — ( + 0.25 exp ^tA(X
T
a(t) = b = 0.01,
T
m л с ЛМ AMc c(t) = —0.5 exp I -^r- I cos
T
T
T
T
+ 0.25 — ( ^ + 0.25 exp ^Ac
T
(12)
Тем самым повторим эксперимент с использованием библиотеки «FDREext» для среды Maple 2019 из статьи автора [32], но у же с использованием нового программного комплекса «FDRE 3.0» для пакета MATLAB.
■а 0.8
Й I
о.е
0.4
0.2
, 1
a) c)
1 ' b)
1 лГ 4Vwf. . J
Рис. 2. Результаты моделирования : а) аппроксимация SA при (13), где: R2 = 0.74, Cor = 0.87 ; b) предсказание SA при (14). Нормированные на максимум значения числа Вольфа, с шагом в 1 месяц : с) с Mая 1996 гг по Август 2021 гг [Figure 2. Simulation results: a) SA approximation at (13), where: R2 = 0.74, Cor = 0.87 ; b) SA prediction at (14). Wolf number values normalized to the maximum, in 1-month increments: c) from May 1996 to August 2021]
2
2
Будем аппроксимировать (рис. 2, синяя кривая) наблюдаемые нормированные данные моделью (12), с параметрами:
N = 304, T = 304, u0 = 0.031109,
Да =-1, Ма = 2.2 П, ©а = 0.125, Ла = 2.85,
Ac =-1, Mc = 2.2 п, 0c = 0.125. (13)
А так же сделаем прогноз (рис. 2, красная кривая) средне значения числа Вольфа на 9 лет, с Августа 2021 гг по Август 2031 гг, с помощью (12) при:
N = 424, T = 424, u0 = 0.031109,
Да = -1.3947, Ma = 3.068 п, ©а = 0.1743, Ла = 3.975,
Дс =-1.3947, Mc = 3.068 п, ©с = 0.1743. (14)
Теперь спустя 14 месяцев после эксперимента, результат которого был опубликован в [32] в Феврале 2022 гг, можем частично проверить предсказание на период из (рис. 1, (b) - чёрная кривая) с Августа 2021 гг по Октябрь 2022 гг.
Для чего параметры скорректированной модели (12), возьмём такими:
N = 318, T = 318, u0 = 0.031109,
Да = -1.046, Ма = 2.301 П, ©а = 0.13, Ла = 2.98,
Дc = -1.046, Mc = 2.301 п, ©c = 0.13. (15)
Из (рис. 3, (a) - лазурная кривая) видно что модельная кривая полученная при помощи (12) недостаточно хорошо отражает фрактальность процесса изменения SA, которую мы подразумеваем и закладываем в модельные уравнения. Что выражается в том что для 1-го всплеска числа Вольфа (23-цикл SA), ненормированная модельная кривая визуально хорошо ложится на нормированные данные, но дальше, для 2-го и начала 3-го всплеска, всё хуже. А так же в том, что спад SA проходит дольше чем подъём, а модельная кривая этого совсем не отражает. Что подтверждается падением значений коэффициентов (Cor) и (R2), которые в дальнейшем должны будут ещё упасть.
Результаты моделирования SA с помощью (12) были опубликованы в статье [32], высоко рейтингового журнала Fractal and Fractional квартиля Q1 с им-пакт фактором IF = 3.396 за последние 5 лет.
Рис. 3. Результаты моделирования : а) частичная проверка предсказания SA (аппроксимация) из (рис. 2, (b) - красная кривая) при (15), где: R2 = 0.72, Cor = 0.86 ; b) предсказание SA при (14). Норм. на макс. значения числа Вольфа, с шагом в 1 месяц : с) с Мая 1996 гг по Август 2021 гг; d) с Мая 1996 по Октябрь 2022 гг [Figure 3. Simulation results : a) partial check of SA prediction (approximation) from (Fig. 2, (b) - red curve) at (15), where: R2 = 0.72, Cor =0.86; b) SA prediction at (14). Wolf number values normalized to the maximum, with a step of 1 month: c) from May 1996 to August 2021; d) from May 1996 to October 2022]
Новая модель с корректным учетом эффекта насыщения
Аналогично при помощи регрессионного анализа и его коэффициентов (R2) [30, 31] и (Cor) [29], модель (9) можно скорректировать:
30t(t)u(a) + a(t)x2(t) — b(t)u(t) — c(t) = 0, u(0) = u0,
a(t) = Да cos ^L^ + фа
a(t) = ^ b = 0.01,
2
c(t) = Дс cos( + Фс ) . (16)
Аппроксимируем (рис. 4, синяя кривая) известные данные моделью (16), при:
N = 318, T = 318, u0 = 0.031109,
Да = 0.75, Ма = 2.25 П, Фа = 0, Ла = 2,
П
Дс = 0.25, Мс = 2.25 п, Фс = -. (17)
Рис. 4. Результаты моделирования : а) аппроксимация SA при (17), где: R2 = 0.76, Cor = 0.90 ; b) предсказание SA при (18). Норм. на макс. значения числа Вольфа, с шагом в 1 месяц : с) Май 1996 по Август 2021 гг; d) Май 1996 по Октябрь 2022 гг [Figure 4. Simulation: a) SA approximation at (17), where: R2 = 0.76, Cor = 0.90 ; b) SA prediction at (18). Wolf number values normalized to the maximum, with a step of 1 month: c) from May 1996 to August 2021; d) from May 1996 to October 2022]
Так же с (16) проведём прогноз (рис. 4, красная кривая) за тот же период:
N = 424, Т = 424, и0 = 0.031109,
Аа = 0.75, Ма = 3.099 п, Фа = 0, Ла = 2.7547,
п
Ае = 0.25, Мс = 3.099 п, Фс = -. (18)
6
Но теперь следует обратить внимание на функцию для определения а(1). Так как именно а(1) - функция делает основной вклад в значение Э^и^) дробной производной 3 в каждый момент времени 1 = 0..Т, что и вносит эффект насыщения в модель. Что наглядно видно из (рис. 5), где сопоставляются по точечно модельная кривая и кривая распределения Э^и^) значений в ходе всего численного эксперимента.
Также а(1) - функция задаёт цикличность модельной кривой, позволяя в принципе описывать циклы SA на основе данный по числу Вольфа, что становится ясно из (17) и (18).
В тоже время с(1) - функция "сдвигает"равномерно всё решение вниз или вверх, т.к. не умножается на значение решения в соответствующей ей точке. В отличии от а(1) - функции, которая определяет степень нелинейности для каждой точки решения.
Кроме того из (рис. 2 и 4) видно, что результаты полученные по уточнённой математической модели (16), визуально хорошо согласуются с известными результатами моделирования прогнозов SA. И в тоже время лучше результатов по (12), полученных авторами раннее в [32] исследовании.
Рис. 5. 1) Результаты x(t) моделирования (16) при (17) и (18); 2) распределения значений дробной производной, соответствующих моделей
[Figure 5. 1) Results of x(t) simulations (16) with (17) and (18); 2) the distribution of the values of the fractional derivative corresponding to the models]
Рис. 6. Результаты моделирования иными методами [Figure 6. Simulation results by other methods]
Графики на (рис. 6) получены с сайта spaceweatherlive.com, некоммерческой организации из Бельгии, которая состоит из нескольких веб-сайтов об астрономии, космосе, космической погоде, северном сиянии и связанных с ними темах.
Используемое программное обеспечение
Для решения задач возникающих при моделировании в исследовании [32] использовался разработанный автором программный комплекс, в виде пользовательской библиотеки «FDREext» для среды MAPLE 2021.
В настоящем исследовании, для проведения моделирования динамического процесса: обработки и визуализации данных, численного решения модельной задачи и визуализации результатов, использовался разработанный автором программный комплекс «FDRE 3.0» на языке MATLAB, а все результаты получены с его помощью.
Данный программный комплекс «FDRE 3.0» является обобщением идей «FDREext» реализованных ранее для моделирования процессов дробной динамики, что выражено в следующих качествах программного комплекса:
• большей универсальности - как при работе только с экспериментальными данными, так и только с моделью или совместно;
• заметного удобства - добавления новых вариаций модели и экспериментальных данных, для последующей программной обработки;
• на порядки возросшей скорости работы - так как стало возможно полностью уйти в (10) от символьных вычислений.
Заключение
Показано, что модель (16) при оптимальном (17) выборе соответствующих параметров моделирования a(t) и а(^,Ь(^,сИ, дает хорошее приближение нормированных наблюдаемых данных по числу Вольфа для циклов солнечной активности,
•л
с коэффициентом детерминации R2 = 0.76 и корреляции Cor = 0.90.
Показано, что предложенная в данной статье модель (16) лучше описывает процесс подъёма-спада SA чем изученная ранее (12) модель. Что подтверждается регрессионным анализом результатов моделирования и нормированных данных за один временной промежуток с Мая 1996 года по Октябрь 2022 года. А так же подтверждается лучшим визуальным соответствием модели (16) с прогнозами полученными другими методами в данной области, чем (12).
Кроме того, при иных параметрах (18) модель (16) способна дать некоторый прогноз возможного количества солнечных пятен. Следовательно, можно приблизительно оценить границы текущего 25-го цикла солнечной активности.
Показано, что солнечная активность снижается. И согласно предложенной модели пик солнечной активности ожидается в период с Февраля 2024 года по Август 2024.
Конкурирующие интересы. Авторы заявляют об отсутствии конфликта интересов в отношении авторства и публикации.
Авторский вклад и ответственность. Все авторы внесли свой вклад в эту статью. Авторы несут полную ответственность за предоставление окончательной
версии статьи в печать. Окончательный вариант рукописи был одобрен всеми авторами.
Аббревиатуры
SA Solar Activity
VO variable order
IFDS Implicit Finite-difference Method
EFDS Explicit Finite-difference Method
NM Newton Method
MNM Modified Newton Method
Список литературы
1. Муртазов А. К. Физика земли. Космические воздействия на геосистемы 2-е изд. пер. и доп.. Москва: Юрайт, 2021.268 с. ISBN 978-5-534-11473-7.
2. Volterra V. Functional theory, integral and integro-differential equations. osow: Science, 1982 ISBN 9780598446336..
3. Учайкин В. В. Метод дробных производных. Ульяновск: Артишок, 2008.510 с. ISBN 9785904198015.
4. Нахушев А.М. Дробное исчисление и его применение. Москва: Физматлит, 2003.272 с. ISBN 5-9221-0440-3.
5. Kilbas A.A., Srivastava H.M., Trujillo J.J. Theory and Applications of Fractional Differential Equations. Amsterdam: Elsevier Science Limited, 2006. 523 pp. ISBN 9780444518323.
6. Uchaikin V. V. Fractional Derivatives for Physicists and Engineers. Vol. I. Background and Theory. Berlin: Springer, 2013.373 pp. ISBN 978-3-642-33911-0 DOI: 10.1007/978-3-642-33911-0.
7. Ortigueira M.D., Valerio D., Machado J. T. Variable order fractional systems, Communications in Nonlinear Science and Numerical Simulation, 2019. vol.71, pp. 231-243 DOI: 10.1016/j.cnsns.2018.12.003.
8. Patnaik S., Hollkamp J. P., Semperlotti F. Applications of variable-order fractional operators: a review, Proceedings of the Royal Society A, 2020. vol.476, no. 2234, pp. 20190498 DOI: 10.1098/rspa.2019.0498.
9. Coimbra C. F. M. Mechanics with variable-order differential operators, Annalen der Physik, 2003. vol. 12, no. 11-12, pp. 692-703 DOI: 10.1002/andp.200310032.
10. Parovik R. I. On a finite-difference scheme for an hereditary oscillatory equation, Journal of Mathematical Sciences, 2021. vol.253, no. 4, pp. 547-557 DOI: 10.1007/s10958-021-05252-2.
11. Parovik R. I. Mathematical modeling of linear fractional oscillators, Mathematics, 2020. vol. 8, no. 11, pp. 18-79 DOI: 10.3390/math8111879.
12. Parovik R. I. Tverdyi D. A. Some Aspects of Numerical Analysis for a Model Nonlinear Fractional Variable Order Equation, Mathematical and Computational Applications, 2021. vol. 26, no. 3, pp. 55 DOI: 10.3390/mca26030055.
13. Sun H., et al. Finite difference schemes for variable-order time fractional diffusion equation, International Journal of Bifurcation and Chaos, 2012. vol.22, no. 04, pp. 1250085 DOI: 10.1142/S021812741250085X.
14. Tverdyi D. A., Parovik R. I. Investigation of Finite-Difference Schemes for the Numerical Solution of a Fractional Nonlinear Equation, Fractal and Fractional, 2022. vol. 6(1), no. 23, pp. 1-27 DOI: 10.3390/fractalfract6010023.
15. Бураев А. В. Некоторые аспекты математического моделирования региональных проявлений солнечной активности и их связи с экстремальными геофизическими процессами, Доклады Адыгской (Черкесской) Международной Академии Наук, 2010. Т. 12, №1, С. 88-90.
16. Постан М. Я. Обобщенная логистическая кривая: ее свойства и оценка параметров, Экономика и математические методы, 1993. Т. 29, №2, С. 305-310.
17. Therese A. S. Generalized Logistic Models, Journal of the American Statistical Association, 1988. vol.83, no. 402, pp. 426-431 DOI: 10.1080/01621459.1988.10478613.
18. Rzkadkowski G., Sobczak L. A generalized logistic function and its applications, Foundations of Management, 2020. vol. 12, no. 1, pp. 85-92 DOI: 10.2478/fman-2020-0007.
19. Mandelbrot B. B. The fractal geometry of nature. New York: WH freeman, 1982. 468 pp. ISBN 0716711869.
20. Drozdyuk A. V. Logistic curve. Toronto: Choven, 2019. 270 pp. ISBN 978-0-9866300-2-6.
21. Твёрдый Д. А. Паровик Р.,И. Математическое моделирование некоторых логистических законов с помощью эредитарной динамической системы Риккати / Материалы 11 Всероссийской научной конференции с международным участием (27-30 мая 2019 г.)., Математическое моделирование и краевые задачи. Самара, СамГТУ, 2019, С. 348-352.
22. Taogetusang, Sirendaoerji, Li S. New application to Riccati equation, Chinese Physics B, 2010. vol. 19, pp. 080303 DOI: 10.1088/1674-1056/19/8/080303.
23. Куркин А. А., Куркина О.Е., Пеленовский Е. Н. Логистические модели распространения эпидемий, Труды НГТУ им. Р.Е. Алексеева., 2020. Т. 129, С. 9-18.
24. Torres-Hernandez A., et al. Fractional Newton-Raphson Method Accelerated with Aitken's Method, Axioms, 2021. vol. 10, no. 2, pp. 1-47 DOI: 10.3390/axioms10020047.
25. Жуков С. А. О пьезокерамике и перспективах ее применения, Мир техники и технологий: международный промышленный журнал, 2021. №5, С. 56-60.
26. Bayldon J. M., Daniel I. M. Flow modeling of the VARTM process including progressive saturation effects, Composites Part A: Applied Science and Manufacturing, 2009. vol. 40, no. 8, pp. 1044-1052 DOI: 10.1016/j.compositesa.2009.04.008.
27. Landis C. M. On the strain saturation conditions for polycrystalline ferroelastic materials, Jornal of Applied Mechanics., 2009. vol.70, no. 4, pp. 470-478 DOI: 10.1115/1.1600472.
28. Sunspot Index and Long-term Solar Observations. (ROB) Av. Circulaire, 3 - B-1180 Brussels, : 18.11.2022 https://www.sidc.be/silso/datafiles.
29. Cox D. R. Hinkley D.V. Theoretical Statistics, 1st edition. London: Chapman & Hall/CRC, 1979.528 pp. ISBN 9780412161605.
30. Hughes A. J., Grawoig D.E. Statistics: A Foundation for Analysis. Boston: Addison Wesley, 1971.525 pp. ISBN 978-0201030211.
31. Chicco D., Warrens M. J., Jurman G. The coefficient of determination R-squared is more informative than SMAPE, MAE, MAPE, MSE and RMSE in regression analysis evaluation, PeerJ Computer Scienc, 2021. vol.299, pp. e623 DOI: 10.7717/peerj-cs.623.
32. Tverdyi D. A., Parovik R. I. Application of the Fractional Riccati Equation for Mathematical Modeling of Dynamic Processes with Saturation and Memory Effect, Fractal and Fractional, 2022. vol. 6(3), no. 163, pp. 1-35 DOI: 10.3390/fractalfract6030163.
Твёрдый Дмитрий АлександровичА - кандидат физико-математических наук, ведущий программист лаборатории электромагнитных излучений института космофизических исследований и распространения радиоволн ДВО РАН, Паратунка, Россия, ОЯСГО 0000-0001-6983-5258.
Паровик Роман Иванович А - доктор физико-математических наук, доцент, ведущий научный сотрудник лаборатории моделирования физических процессов института космофизических исследований и распространения радиоволн ДВО РАН, Паратунка, Россия, ОЯСГО 0000-0002-1576-1860.
Vestnik KRAUNC. Fiz.-Mat. Nauki. 2022. vol. 41. no. 4. pp. 47-65. ISSN 2079-6641
MSC 26A33, 85-08 Research Article
Mathematical modeling in MATLAB of solar activity cycles according to the growth-decline of the Wolf number
D.A. Tverdyi, R.I. Parovik
Institute of Cosmophysical Research and Radio Wave Propagation FEB RAS, 684034, p. Paratunka, Mirnaya st., 7, Russia E-mail: dimsolid95@gmail.com
In this article, mathematical modeling of the dynamics of solar activity is carried out. Observational data on the average monthly number of sunspots, called the Wolf number, for the period of 24.5 years from May 1996 to October 2022 are studied. Based on the results of a similar study of data on this process, using the Riccati equation of a fractional constant order, that the rise and fall of the Wolf number over time occurs along a curve very close to the generalized logistic curve, this article also proposes a mathematical model based on the Riccati equation. Since the Riccati equation describes well the processes that obey the logistic law. However, the equation is generalized to the integro-differential Riccati equation by introducing a fractional derivative of the Gerasimov-Caputo type of variable order, and a fractional derivative with a variable order, allows you to get a more precise mathematical model of Wolf number cycles with saturation, and allows you to take into account the effect of variable memory. All model calculations, data processing and visualization are carried out in the FDRE 3.0 program developed in the MATLAB package. Modeling parameters are refined by approximation of known data under study, using regression analysis. As a result, the model curves and graphs of the observed data known for 24.5 years show good agreement with each other. With the help of a refined mathematical model, a forecast is made for the next 9 years, which visually agrees well with the known model results of solar activity.
Key words: solar activity, Wolf number, mathematical modeling, saturation effect, hereditarity, Riccati equation, Gerasimov-Caputo type derivative.
d DOI: 10.26117/2079-6641-2022-41-4-47-65
Original article submitted: 19.11.2022 Revision submitted: 22.11.2022
For citation. Tverdyi D. A., Parovik R. I. Mathematical modeling in MATLAB of solar activity cycles according to the growth-decline of the Wolf number. Vestnik KRAUNC. Fiz.-mat. nauki. 2022, 41: 4,47-65. d DOI: 10.26117/2079-6641-2022-41-4-47-65
Competing interests. The authors declare that there are no conflicts of interest regarding authorship and publication.
Contribution and Responsibility. All authors contributed to this article. Authors are solely responsible for providing the final version of the article in print. The final version of the manuscript was approved by all authors.
The content is published under the terms of the Creative Commons Attribution 4.0 International License (https://creativecommons.org/licenses/by/4.0/deed.ru)
© Tverdyi D. A., Parovik R.I., 2022
Funding. The research was carried out within the framework of the state task of IKIR FEB RAS on the topic AAAA-A21-121011290003-0
ISSN 2079-6641
TBepgbiH A-A., napoBHK P.M.
References
[1] Murtazov A. K. Fizika zemli. Kosmicheskie vozdejstviya na geosistemy 2-e izd. per. i dop. [Physics of the earth. Space impacts on geosystems 2nd ed. additional]. Moscow, YUrajt, 2021, 268 (In Russian).
[2] Volterra V. Functional theory, integral and integro-differential equations. Moscow, Nauka, 1982 (In Russian).
[3] Uchajkin V. V. Metod drobnyh proizvodnyh [Fractional derivative method]. Ulyanovsk, Artichoke, 2008, 510 (In Russian).
[4] Nahushev A. M. Drobnoe ischislenie i ego primenenie [Fractional calculus and its application]. Moscow, Fizmatlit, 2003, 272 (In Russian).
[5] Kilbas A.A., Srivastava H.M., Trujillo J.J. Theory and Applications of Fractional Differential Equations. Amsterdam, Elsevier Science Limited, 2006, 523.
[6] Uchaikin V. V. Fractional Derivatives for Physicists and Engineers. Vol. I. Background and Theory. Berlin, Springer, 2013, 373, DOI: 10.1007/978-3-642-33911-0.
[7] Ortigueira M.D., Valerio D., Machado J.T. Variable order fractional systems, Communications in Nonlinear Science and Numerical Simulation, 2019, 71, 231-243. DOI: 10.1016/j.cnsns.2018.12.003.
[8] Patnaik S., Hollkamp J. P., Semperlotti F. Applications of variable-order fractional operators: a review, Proceedings of the Royal Society A, 2020, 476:2234, 20190498. DOI: 10.1098/rspa.2019.0498.
[9] Coimbra C. F. M. Mechanics with variable-order differential operators, Annalen der Physik, 2003, 12: 11-12, 692-703. DOI: 10.1002/andp.200310032.
[10] Parovik R.I. On a finite-difference scheme for an hereditary oscillatory equation, Journal of Mathematical Sciences, 2021, 253:4, 547-557. DOI: 10.1007/s10958-021-05252-2.
[11] Parovik R.I. Mathematical modeling of linear fractional oscillators, Mathematics, 2020, 8:11, 18-79. DOI: 10.3390/math8111879.
[12] Parovik R. I. Tverdyi D. A. Some Aspects of Numerical Analysis for a Model Nonlinear Fractional Variable Order Equation, Mathematical and Computational Applications, 2021, 26:3, 55. DOI: 10.3390/mca26030055.
[13] Sun H., et al. Finite difference schemes for variable-order time fractional diffusion equation, International Journal of Bifurcation and Chaos, 2012, 22:4, 1250085. DOI: 10.1142/S021812741250085X.
[14] Tverdyi D. A., Parovik R.I. Investigation of Finite-Difference Schemes for the Numerical Solution of a Fractional Nonlinear Equation, Fractal and Fractional, 2022, 6(1):23, 1-27. DOI: 10.3390/fractalfract6010023.
[15] Buraev A. V. Nekotorye aspekty matematicheskogo modelirovaniya regional'nyh proyavlenij solnechnoj aktivnosti i ih svyazi s ekstremal'nymi geofizicheskimi processami [Some aspects of mathematical modeling of regional manifestations of solar activity and their connection with extreme geophysical processes]. Doklady Adygskoj (CHerkesskoj) Mezhdunarodnoj Akademii Nauk [Reports of the Adyghe (Circassian) International Academy of Sciences]. 2010, vol. 12, no. 1, pp. 88-90 (In Russian).
[16] Postan M.Y. Generalized logistic curve: its properties and estimation of parameters, Ekonomika i matematicheskie metody, 1993, 29:2. 305-310 (In Russian).
[17] Therese A. S. Generalized Logistic Models, Journal of the American Statistical Association, 1988, 83:402, 426-431. DOI: 10.1080/01621459.1988.10478613.
[18] Rzkadkowski G., Sobczak L. A generalized logistic function and its applications, Foundations of Management, 2020, 12:1, 85-92. DOI: 10.2478/fman-2020-0007.
[19] Mandelbrot B. B. The fractal geometry of nature. New York, WH freeman, 1982, 468.
[20] Drozdyuk A. V. Logistic curve. Toronto, Choven, 2019, 270 pp., isbn: 978-0-9866300-2-6.
[21] Tverdyi D.A., Parovik R.I. Mathematical modeling of some logistic laws using the hereditary dynamical Riccati system, In: Proceedings of the 11th All-Russian Scientific Conference with International Participation (May 27-30, 2019). Matem. Mod. Kraev. Zadachi. Samara, Samara State Technical Univ., 2019, pp. 348-352 (In Russian).
[22] Taogetusang, Sirendaoerji, Li S. New application to Riccati equation, Chinese Physics B, 2010, 19, 080303. DOI: 10.1088/1674-1056/19/8/080303.
[23] Kurkin A. A., Kurkina O. E., Pelenovskij E. N. Logistic models of the spread of epidemics, Proceedings of NSTU im. R.E. Alekseev., 2020, 129, 9-18.(In Russian).
[24] Torres-Hernandez A., et al. Fractional Newton-Raphson Method Accelerated with Aitken's Method, Axioms, 2021, 10:2, 1-47. DOI: 10.3390/axioms10020047.
[25] Zhukov S.A. O p'ezokeramike i perspektivah ee primeneniya [On piezoceramics and the prospects for its application]. Mir tekhniki i tekhnologij: mezhdunarodnyj promyshlennyj zhurnal [World of Engineering and Technology: International Industrial Journal]. 2021, no. 5, pp. 56-60 (In Russian).
[26] Bayldon J.M., Daniel I. M. Flow modeling of the VARTM process including progressive saturation effects, Composites Part A: Applied Science and Manufacturing, 2009, 40:8, 1044-1052. DOI: 10.1016/j.compositesa.2009.04.008.
[27] Landis C. M. On the strain saturation conditions for polycrystalline ferroelastic materials, Jornal of Applied Mechanics, 2009, 70:4, 470-478. DOI: 10.1115/1.1600472.
[28] Sunspot Index and Long-term Solar Observations, dataset by Royal Observatory of Belgium (ROB) Av. Circulaire, 3 - B-1180 Brussels, acsses: 18.11.2022, https://www.sidc.be/silso/datafiles.
[29] Cox D. R. Hinkley D. V. Theoretical Statistics, 1st edition. London, Chapman & Hall/CRC, 1979, 528.
[30] Hughes A. J., Grawoig D. E. Statistics: A Foundation for Analysis. Boston, Addison Wesley, 1971, 525.
[31] Chicco D., Warrens M. J., Jurman G. The coefficient of determination R-squared is more informative than SMAPE, MAE, MAPE, MSE and RMSE in regression analysis evaluation, PeerJ Computer Sci., 2021, 299, e623. DOI: 10.7717/peerj-cs.623.
[32] Tverdyi D. A. Parovik R. I. Application of the Fractional Riccati Equation for Mathematical Modeling of Dynamic Processes with Saturation and Memory Effect, Fractal and Fractional, 2022, 6(3), 163, 1-35. DOI: 10.3390/fractalfract6030163.
Tverdyi Dmitrii AlexsandrovichA - Ph. D. (Phys. & Math.), Lead programmer laboratory of electromagnetic propogation Institute of Cosmophysical Research and Radio Wave Propagation FEB RAS, Paratunka, Russia, ORCID 0000-0001-6983-5258.
Parovik Roman IvanovichA - D. Sci. (Phys. & Math.), Associate Professor, Leading researcher laboratory of modeling physical processes Institute of Cosmophysical Research and Radio Wave Propagation FEB RAS, Paratunka, Russia, ORCID 0000-0002-1576-1860.
