КОМПЛЕКС ПРОГРАММ VOFDDE 1.0 ДЛЯ КОЛИЧЕСТВЕННОГО И КАЧЕСТВЕННОГО АНАЛИЗА ДРОБНОГО ОСЦИЛЛЯТОРА ДУФФИНГА

КОМПЛЕКС ПРОГРАММ VOFDDE 1.0 ДЛЯ КОЛИЧЕСТВЕННОГО И КАЧЕСТВЕННОГО АНАЛИЗА ДРОБНОГО ОСЦИЛЛЯТОРА ДУФФИНГА

Ким Валентин Александрович

младший научный международной интегративной научно-исследовательской лабораторией экстремальных явлений Камчатки, Камчатский государственный университет имени Витуса Беринга,

РФ, г. Петропавловск-Камчатский

Паровик Роман Иванович

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

РФ, Камчатский край, с. Паратунка

 

VOFDDE 1.0 SOFTWARE FOR QUANTITATIVE AND QUALITATIVE ANALYSIS OF FRACTIONAL DUFFING OSCILLATOR

Valentin Kim

Junior Researcher, International Integrative Research Laboratory of Extreme Phenomena of Kamchatka, Vitus Bering Kamchatka State University,

Russia, Petropavlovsk-Kamchatsky

Roman Parovik

Doctor of Physical and Mathematical Sciences, Head of the International Integrative Research Laboratory of Extreme Phenomena of Kamchatka, Vitus Bering Kamchatka State University, Petropavlovsk-Kamchatsky, Leading Researcher of the Laboratory for Simulation of Physical Processes Institute of Cosmophysical Research and Radio Wave Propagation, Far Eastern Branch of the Russian Academy of Sciences,

Russia, Kamchatka Territory, Paratunka

 

АННОТАЦИЯ

В работе подробно рассмотрены процедуры библиотеки VOFDDE 1.0, написанной на языке символьной математики Maple для количественного и качественного анализа дробного осциллятора Дуффинга с производной дробного переменного порядка. Дается описание каждой процедуры на языке Maple и для некоторых из них дана визуализация. Комплекс программ VOFDDE 1.0 был зарегистрирован в Роспатенте.

ABSTRACT

The paper considers in detail the procedures of the VOFDDE 1.0 library written in the language of symbolic mathematics Maple for the quantitative and qualitative analysis of the fractional Duffing oscillator with a fractional variable order derivative. A description of each procedure in the Maple language is given and visualization is given for some of them. The VOFDDE 1.0 software package has been registered with Rospatent.

 

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

Рассмотрим следующую задачу Коши для нелинейного осцилляционного уравнения:

                (1)

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

Оператор дробного переменного порядка  понимается в смысле Римана-Лиувилля:

Для количественного и качественного анализа задачи (1) при  была разработана библиотека VOFDDE 1.0.

Комплекс программ VOFDDE 1.0

В этом разделе представлен программный комплекс для среды символьной математики MAPLE 2021, разработанный для решения задач: математического моделирования, сопоставления результата с реальными данными наблюдений, а также визуализации результатов [1]. Данный программный комплекс включает в себя: исполняемый файл, определяющий логику программы, а также библиотеку <<VOFDDEext>> содержащую необходимые для исследования функции.

Комплекс программ VOFDDE 1.0 (сокращение от Variable Order Fractional Differentional Equation) позволяет получить решение задачи Коши (1), строить спектры максимальных показателей Ляпунова, а также амплитудно-частотных (АЧХ) и фазо-частотных (ФЧХ) характеристик в виде массива данных и записывать данные в текстовый файл, а пользователь может визуализировать эти данные в виде графиков или таблиц. VOFDDE 1.0. содержит 9 основных процедур: Duf, Oscillogram, LapLam, AFC, AnalythicalAFC, APFC3d, Quality, WriteFile_txt, EFDS().  

Основное отличие разработанного нового программного комплекса, является возможность использовать широкий функционал среды Maple 2021, что отразилось в следующем:

• реализация нелокальной явной конечно-разностной схемы (ЯКРС), рассмотренной в [2]-[5];
• реализация нелокальной неявная конечно-разностная схема (НКРС);
• реализация модифицированный метод Ньютона (ММН), для решения схемы НКРС;
• реализация алгоритма Беннетина (Вольфа) для построения спектров максимальных показателей Ляпунова.
• визуализация результатов исследования вынужденных колебаний (построение АЧХ, ФЧХ и добротности)
• реализация численного анализа схем: НКРС, ЯКРС.

 

Рисунок 1. Процедуры библиотеки VOFDDEext

 

Файл с кодом вызова библиотеки имеет вид.

Файл: Code_VOFDDE_lib,load.mw. Основной интерфейс «Worksheet».

restart;
currentdir(worksheetdir);
mylibdir:= cat(kernelopts(homedir), kernelopts(dirsep), "maple", kernelopts(dirsep), "toolbox", kernelopts(dirsep), "personal", kernelopts(dirsep), "lib");
FileTools:-MakeDirectory(mylibdir, 'recurse');
LibraryTools:-Create(cat(mylibdir, kernelopts(dirsep), "packages.mla"));
libname := mylibdir, libname;
libname := libname[1], libname[2];
libname[1];
LibraryTools:-Save('VOFDDEext', libname[1]);
restart;
with(VOFDDEext);

Процедура Duf()

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

Duf := proc(T, N, lambda, omega, eta, delta, w) – создание процедуры.

local x, y, q, h, A, c, i, j, k, f; – объявление локальных переменных.

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

h := evalf(T/N); A := [];(пустой массив данных) B := []; x[0] := 0; y[0] := 0; x[1] := h*y[0] + x[0]; – присвоение значений локальным переменным.

for k from 0 to N do

         q[k] := evalf(0.8*cos(0.5*k*h)^2);

end do; – определение значений дробной производной переменного порядка.

for k from 0 to N do

         c[k, 0] := 1;

         for j from 1 to k do

                   c[k, j] := (1-(1+q[k])/j)*c[k, j-1];

         end do;

end do; – определение значений коэффициентов Грюнвальда-Летникова.

for k from 1 to N-1 do

         x[k+1] := (2-lambda*h^(2-q[k])-h^2*omega)*x[k]-x[k-1]-h^2*(eta*x[k]^3-delta*cos(w*k*h))-lambda*h^(2-q[k])*add(c[k, i]*x[k-i], i = 1..k);

         y[k] := (x[k+1]-x[k])/h;

         A := [op(A), [x[k+1], y[k]]]; (запись в массив данных).

end do; – создание массива данных численного решения задачи Коши (1).

return A; – возвращение массива данных.

end proc; – завершение создания процедуры.

Вызов и визуализация процедуры в виде фазовой траектории изображена на рис. 2.

 

Рисунок 2. Визуализация результатов процедуры Duf в виде фазовой траектории

Процедура Oscillogram()

Также, как и Duf, процедура Oscillogram, выводит на экран массив численного решения, только с отличием в том, что вместе со значениями численного решения, выводятся и моменты времени, в которые эти решения были получены. Этот массив можно визуализировать в виде осциллограммы. Процедура представлена в виде кода.

Oscillogram := proc(T, N, lambda, omega, eta, delta, w)

local  x, y, q, h, A, c, i, j, k; - локальные переменные

h := evalf(T/N); A := []; x[0] := 0; y[0] := 0; x[1] := h*y[0]+x[0];

for k from 0 to N do

         q[k] := evalf(0.8*cos(0.5*k*h)^2);

end do;

for k from 0 to N do

         c[k, 0] := 1;

         for j from 1 to k do

                   c[k, j] := (1-(1+q[k])/j)*c[k, j-1];

         end do;

end do;

for k from 1 to N-1 do

         x[k+1] := (2-lambda*h^(2-q[k])-h^2*omega)*x[k]-x[k-1]-h^2*(eta*x[k]^3-delta*cos(w*k*h))-lambda*h^(2-q[k])*add(c[k, i]*x[k-i], i = 1..k);

         A := [op(A), [k*h, x[k]]];

end do; - создание массива данных численного решения

return A;– возвращение массива данных

end proc;

Вызов процедуры

 

Рисунок 3. Визуализация результатов процедуры Oscillogram в виде осциллограммы

 

Массив значений максимальных показателей Ляпунова

Процедура LapLam запускается в виде LapLam(T, N, , , , ). Данная процедура возвращает нам массив показателей Ляпунова относительно коэффициента . Код процедуры:

LapLam := proc(T, N, omega, eta, delta, w)

local  lambda, x, y, xx, yy, XX, YY, Vx, Vy, Vx1, Wy, VxNorm, WyNorm, Vy1, q, h, A, c, i, j, k, l, n, summa, Lap1, Lap2; – локальные переменные.

h := evalf(T/N); A := []; B := []; x[0] := 0; y[0] := 0; xx[0] := 1; yy[0] := 0; XX[0] := 0; YY[0] := 1; n := 100;

for k from 0 to N do

         q[k] := evalf(0.8*cos(0.5*k*h)^2);

end do;

for k from 0 to N do

         c[k, 0] := 1;

         for j from 1 to k do

                   c[k, j] := (1-(1+q[k])/j)*c[k, j-1];

         end do;

end do;

for l from 1 to n do

   summa[1] := 0;

   summa[2] := 0;

   lambda[l] := 0.1*l;

   for k from 1 to N-1 do

           x[k] := h*y[k-1]+x[k-1];

           y[k]:=y[k-1]+h*(-omega^2*x[k-1]-eta*x[k-1]^3+delta*cos(w*h*(k-1))-lambda[l]*h^(-q[k])*add(c[k, i]*x[k-i], i = 0..k-1)); – численная схема

           xx[k]:=h*yy[k-1]+xx[k-1];

           yy[k]:=yy[k-1]+h*((-3*eta*x[k-1]^2-omega^2)*xx[k-1]-lambda[l]*h^(-q[k])*add(c[k, i]*xx[k-i], i = 0..k-1)); – уравнения в вариациях полученные по формуле [2]

           XX[k] := h*YY[k-1]+XX[k-1];

           YY[k]:=YY[k-1]+h*((-3*eta*x[k-1]^2-omega^2)*XX[k-1]-lambda[l]*h^(-q[k])*add(c[k, i]*XX[k-i], i = 0..k-1)); – вторая пара уравнений в вариациях.

           Vx:= Vector([xx[k], yy[k]]);

Vy:= Vector([XX[k], YY[k]]); – векторы значений вариационных переменных

           Vx1:= Vx/VectorNorm(Vx, 2);

           Wy := Vy-DotProduct(Vy, Vx1)*Vx1;

           VxNorm := VectorNorm(Vx, 2);

    WyNorm := VectorNorm(Wy, 2);

    Vy1 := Wy/WyNorm; – ортогонализация Грамма-Шмидта

    summa[1] := summa[1]+ln(VxNorm);

    summa[2] := summa[2]+ln(WyNorm);

     xx[k+1] := Vx[1]/VxNorm;

     yy[k+1] := Vx[2]/VxNorm;

     XX[k+1] := Wy[1]/WyNorm;

     YY[k+1] := Wy[2]/WyNorm;

    end do; – алгоритм Беннетина.

          Lap1[l] := summa[1]/N;

          Lap2[l] := summa[2]/N; – вычисление показателей Ляпунова

          A := [op(A), [lambda[l], Lap1[l]]]; – запись в массив

         end do;

         return A;

         end proc; – возвращение массива данных

Пример вызова процедуры изображен на рис. 4.

 

Рисунок 4. Спектр максимальных показателей Ляпунова.

 

Визуализация AFC

Резонансные кривые строятся с помощью процедур AFC() и AnalythicalAFC(). Запуск: AFC(T, N, M, , , , ), AnalythicalAFC(, , , ), где M - количество прогонок. AFC возвращает массив значений АЧХ, рассчитанных численно, а AnalythicalAFC строит кривую АЧХ по аналитической формуле [2]. Процедурный код:

AFC:= proc(T, N, M, lambda, omega, eta, delta)

local  x, xx, xxx, y, yy, yyy, q, qq, qqq, h, h1, A, AA, AAA, A1, AA1, AAA1, AFC1, AFC2, AFC3, c, cc, ccc, w, i, j, k, f, beta, m, q1, q2, q3, p1, c1, U1, W1, p2, c2, U2, W2, p3, c3, U3, W3, rr4, rrr4, rrrr4;

h := evalf(T/N); h1 := evalf(3/M); A1 := []; AA1 := []; AAA1 := []; x[0] := 0; y[0] := 0; xx[0] := x[0]; yy[0] := y[0]; xxx[0] := x[0]; yyy[0] := y[0]; x[1] := h*y[0] + x[0]; xx[1] := x[1]; xxx[1] := x[1];

   for k from 0 to N do

    q[k] := evalf(0.8*cos(0.5*k*h)^2);

    qq[k] := evalf(1.5/5+1/(2*T)*k*h);

    qqq[k] := evalf(4/5-1/(2*T)*k*h);

   end do; – определение функций дробного порядка производной, обладающих свойствами периодичности и монотонности.

         for k from 0 to N do

          c[k, 0] := 1;

          cc[k, 0] := 1;

          ccc[k, 0] := 1;

                   for j to k do

          c[k, j] := (1-(1+q[k])/j)*c[k, j-1];

          cc[k, j] := (1-(1+qq[k])/j)*cc[k, j-1];

          ccc[k, j] := (1-(1+qqq[k])/j)*ccc[k, j-1];

                   end do;

         end do;

     for j from 0 to M do

          w[j] := j*h1;

          for k from 1 to N-1 do

      x[k+1] := (2-lambda*h^(2-q[k])-h^2*omega)*x[k]-x[k-1]-h^2*(eta*x[k]^3-delta*cos(w[j]*k*h))-lambda*h^(2-q[k])*add(c[k, i]*x[k-i], i = 1..k);

     y[k] := (x[k+1]-x[k])/h;

     xx[k+1]:=(2-lambda*h^(2-qq[k])-h^2*omega)*xx[k]-xx[k-1]-h^2*(eta*xx[k]^3-delta*cos(w[j]*k*h))-lambda*h^(2-qq[k])*add(cc[k, i]*xx[k-i], i = 1..k);

     yy[k]:=(xx[k+1]-xx[k])/h;

   xxx[k+1]:=(2-lambda*h^(2-qqq[k])-h^2*omega)*xxx[k]-xxx[k-1]-h^2*(eta*xxx[k]^3-delta*cos(w[j]*k*h))-lambda*h^(2-qqq[k])*add(ccc[k, i]*xxx[k-i], i = 1..k);

   yyy[k] := (xxx[k+1]-xxx[k])/h;

    end do; – численное решение для переменных порядков дробной производной, обладающих свойствами периодичности и монотонности.

          A[j] := max(seq(abs(x[f]), f = 1500..N));

          AA[j] := max(seq(abs(xx[f]), f = 1500..N));

          AAA[j] := max(seq(abs(xxx[f]), f = 1500..N)); – для разных схем выбираются максимальные значения амплитуд колебаний на больших отрезках времени.

    A1 := [op(A1), [w[j], A[j]]];

    AA1 := [op(AA1), [w[j], AA[j]]];

    AAA1 := [op(AAA1), [w[j], AAA[j]]]; – включение в массив значений AFC.

         end do;

         return A1, AA1, AAA1;

         end proc; – завершение процедуры численного построения AFC.

AnalythicalAFC:=proc(lambda, omega, eta, delta) – аналитическое построение AFC.

local  w, i, j, k, f, beta, m, q1, q2, q3, p1, c1, U1, W1, p2, c2, U2, W2, p3, c3, U3, W3, rr4, rrr4, rrrr4;

beta := 2;

m := -s^(beta-2)*cos(beta*Pi/2);

q1(t):= 0.8*cos(omega*t);

q2(t):= 1.5/5+1/(2*T)*t;   

q3(t):= 4/5-1/(2*T)*t;

psi(1-q1(t)):=-1/2*(10^(1/3)-1)+sum(1/kk-1/(kk+1-q1(t)), kk = 1..n);

p1:=s^(beta-1)*sin((1/2)*beta*Pi)+lambda*s^(q1(t)-1)*sin((1/2)*q1(t)*Pi)-lambda*s^(q1(t)-2)*(diff(q1(t), t))*((ln(s)-abs(psi(1-q1(t))))*cos((1/2)*q1(t)*Pi)-(1/2)*Pi*sin((1/2)*q1(t)*Pi));

sigma1:=sqrt(xi^beta+lambda*s^q1(t)*cos((1/2)*q1(t)*Pi)+3*A^2*eta*(1/4)+lambda*s^(q1(t)-1)*(diff(q1(t), t))*((ln(s)-abs(psi(1-q1(t))))*sin((1/2)*q1(t)*Pi)+(1/2)*Pi*cos((1/2)*q1(t)*Pi)));

U1 := -m*s^2+sigma1^2;

W1 := p1*s; – определение параметров из [2].

psi(1-q2(t)):=-1/2*(10^(1/3)-1)+sum(1/kk-1/(kk+1-q2(t)), kk = 1..n);

p2:=s^(beta-1)*sin((1/2)*beta*Pi)+lambda*s^(q2(t)-1)*sin((1/2)*q2(t)*Pi)-lambda*s^(q2(t)-2)*(diff(q2(t), t))*((ln(s)-abs(psi(1-q2(t))))*cos((1/2)*q2(t)*Pi)-(1/2)*Pi*sin((1/2)*q2(t)*Pi));

sigma2:=sqrt(xi^beta+lambda*s^q2(t)*cos((1/2)*q2(t)*Pi)+3*A^2*eta*(1/4)+lambda*s^(q2(t)-1)*(diff(q2(t), t))*((ln(s)-abs(psi(1-q2(t))))*sin((1/2)*q2(t)*Pi)+(1/2)*Pi*cos((1/2)*q2(t)*Pi)));

U2 := -m*s^2+sigma2^2;

W2 := p2*s;

psi(1-q3(t)):=-1/2*(10^(1/3)-1)+sum(1/kk-1/(kk+1-q3(t)), kk = 1..n);

p3:=s^(beta-1)*sin((1/2)*beta*Pi)+lambda*s^(q3(t)-1)*sin((1/2)*q3(t)*Pi)-lambda*s^(q3(t)-2)*(diff(q3(t), t))*((ln(s)-abs(psi(1-q3(t))))*cos((1/2)*q3(t)*Pi)-(1/2)*Pi*sin((1/2)*q3(t)*Pi));

sigma3:=sqrt(xi^beta+lambda*s^q3(t)*cos((1/2)*q3(t)*Pi)+3*A^2*eta*(1/4)+lambda*s^(q3(t)-1)*(diff(q3(t), t))*((ln(s)-abs(psi(1-q3(t))))*sin((1/2)*q3(t)*Pi)+(1/2)*Pi*cos((1/2)*q3(t)*Pi)));

U3 := -m*s^2+sigma3^2;

W3 := p3*s;

rr4 := implicitplot(A11 - delta/sqrt(U1^2 + W1^2) = 0, s = 0 .. 3, A11 = 0 .. 4, style = line, color = blue);

rrr4 := implicitplot(A2 - delta/sqrt(U2^2 + W2^2) = 0, s = 0 .. 3, A2 = 0 .. 4, style = line, color = red);

rrrr4 := implicitplot(A3 - delta/sqrt(U3^2 + W3^2) = 0, s = 0 .. 3, A3 = 0 .. 4, style = line, color = green);

return rr4, rrr4, rrrr4; – визуализация аналитических AFC.

end proc;

Пример работы процедур AFC() Analythical(AFC) (рис. 5).

 

Рисунок 5. Численное и аналитическое АЧХ

 

AFC3d()

Процедура APFC3d(). Запуск: APFC3d(, , , , A). Процедура строит AFC и PFC в трехмерном пространстве. Показывает зависимость амплитуды и фазы колебаний от частоты внешней силы и порядка дробной производной. Код процедуры имеет следующий вид.

APFC3d := proc(lambda, omega, eta, delta, AA)

local A, w, i, j, k, f, beta, m, q, p1, c1, U1, W1, p2, c2, U2, W2, rr4, rrr4;

beta := 1.955;

m := -s^(beta - 2)*cos(beta*Pi/2);

q(t) := 0.8*cos(0.5*t)

psi(1-q(t)):=-1/2*(10^(1/3)-1)+sum(1/kk-1/(kk+1-q(t)), kk = 1..n);

p1:=s^(beta-1)*sin((1/2)*beta*Pi)+lambda*s^(q(t)-1)*sin((1/2)*q(t)*Pi)-lambda*s^(q(t)-2)*(diff(q(t), t))*((ln(s)-abs(psi(1-q(t))))*cos((1/2)*q(t)*Pi)-(1/2)*Pi*sin((1/2)*q(t)*Pi));

sigma1:=sqrt(xi^beta+lambda*s^q(t)*cos((1/2)*q(t)*Pi)+3*A^2*eta*(1/4)+lambda*s^(q(t)-1)*(diff(q(t), t))*((ln(s)-abs(psi(1-q(t))))*sin((1/2)*q(t)*Pi)+(1/2)*Pi*cos((1/2)*q(t)*Pi)));

U1 := -m*s^2+sigma1^2;

W1 := p1*s;

psi(1-q3(t)):=-1/2*(10^(1/3)-1)+sum(1/kk-1/(kk+1-q3(t)), kk = 1..n);

p2:=s^(beta-1)*sin((1/2)*beta*Pi)+lambda*s^(q(t)-1)*sin((1/2)*q(t)*Pi)-lambda*s^(q(t)-2)*(diff(q(t), t))*((ln(s)-abs(psi(1-q(t))))*cos((1/2)*q(t)*Pi)-(1/2)*Pi*sin((1/2)*q(t)*Pi));

sigma2:=sqrt(xi^beta+lambda*s^q(t)*cos((1/2)*q(t)*Pi)+3*A^2*eta*(1/4)+lambda*s^(q(t)-1)*(diff(q(t), t))*((ln(s)-abs(psi(1-q(t))))*sin((1/2)*q(t)*Pi)+(1/2)*Pi*cos((1/2)*q(t)*Pi)));

U2 := -m*s^2+sigma2^2;

W2 := p2*s;

rr4 := implicitplot3d(A - delta/sqrt(U1^2 + W1^2) = 0, q = 0 .. 1, s = 0 .. 3, A = 0 .. 4, style = surface, color = blue);

rrr4 := plot3d(-arccot((W2/U2)^(-1)), q2 = 0 .. 1, s = 0 .. 3, style = surface, color = blue); – построение 3d графиков AFC и PFC.

         return rr4, rrr4;

         end proc;

Результат выполнения кода изображен на рис. 6.

 

Рисунок 6. АЧХ и ФЧХ

 

Поверхности добротности

Quality(, , , , A, q). Процедура строит поверхности добротности, в зависимости от порядка дробной производной, амплитуды колебаний и частоты внешней силы.

Quality := proc(lambda, omega, eta, delta, AA, qq)

local A, w, i, j, k, f, beta, m, q, Q1, Q2, rr4, rrr4;            

beta := 2;

m := -s^(beta - 2)*cos(beta*Pi/2);

Q1:=sqrt((omega^beta)-lambda*(s^q)*cos((1/2)*q*Pi)+3*eta*(AA^2)*(1/4))/((s^(beta-1))*sin((1/2)*beta*Pi)+lambda*(s^(q-1))*sin((1/2)*q*Pi)); – зависимость добротности от показателя дробной производной и частоты внешней силы

Q2:=sqrt(omega^beta+lambda*s^qq*cos((1/2)*qq*Pi)+3*eta*A^2*(1/4))/(s^(beta-1)*sin((1/2)*beta*Pi)+lambda*s^(qq-1)*sin((1/2)*qq*Pi)); – зависимость добротности от собственной амплитуды и частоты внешней силы

rr4 := plot3d(Q1, s = 0 .. 3, q = 0 .. 1, style = surface, color = blue);

rrr4 := plot3d(Q2, s = 0 .. 3, A = 0 .. 5, style = surface, color = blue);

return rr4, rrr4; – визуализация поверхностей добротности (рис. 7)

end proc;

 

Рисунок 7. Добротность

 

Запись в файл

Процедура WriteFile txt. Запуск: WriteFile txt(result solving::list:=[], path::string:=""), где result solving::list:=[] - данные, записываемые в файл, path::string:="" - имя файла и путь к нему. Процедура позволяет представить результаты в виде текстового файла (рис. 8.).

WriteFile_txt:=proc(result_solving::list:=[], path::string:="")

         uses FileTools;

         local fw, i, path_to_write;

         if (nops(result_solving) = 0) then

                   error "reboot";

         end if;

                   if ( length(path) > 0 ) then

                            path_to_write := path;

                   else

                            path_to_write := CreateFileName();

                            path_to_write := cat(path_to_write, ".txt");

                   end if; – создание чистого файла

                   if (nops(result_solving) > 0) then

                            #-- write in file –

fw := fopen(path_to_write, WRITE, TEXT); – запись данных результатов процедур.

                            #print(path_t        o_write);

                            for i to nops(result_solving) do

                                      fprintf(fw, convert(result_solving[i], string));

                                      if i < nops(result_solving) then

                                               fprintf(fw, "\n");

                                      end if;

                            end do;

                            fclose(fw);

                   end if;

         end proc;

 

Рисунок 8. Пример выгружаемого текстового файла

        

Метод Ньютона для НКРС

Процедура EFDS(). Запуск: EFDS(N, T, , , , , , ). Процедура выводит массив данных численного решения задачи Коши (1), полученного по неявной схеме [6] с помощью модифицированного метода Ньютона. Код процедуры.

EFDS := proc(T, N, lambda, omega, b, delta, phi, epsilon)

local  x, y, q, tau, h, r, A, B, c, i, j, k, n, m, f, ff, JJ, XX0, XX1, FF0, JJ0;

tau := evalf(T/N); A := []; B := []; x[0] := 0; y[0] := 0; x[1] := tau*y[0] + x[0]; ff[1] := x[1]-tau*y[0]-x[0];

for k from 0 to N+1 do

         q[k] := evalf(0.8*cos(0.5*k*tau)^2);

end do;

for k from 0 to N+1 do

         c[k, 0] := 1;

         for j from 1 to k do

                   c[k, j] := (1 - (1 + q[k])/j)*c[k, j - 1];

         end do;

end do;

for k from 1 to N do

ff[k+1]:=(x[k+1]-2*x[k]+x[k-1])/tau^2-delta*evalf(cos(phi*tau*(k+1)))+lambda*tau^(-q[k+1])*add(c[k+1,i]*x[k-i+1],i=0..k+1)+x[k+1]^3+x[k+1];

end do; – создание итерационной функции для EFDS [6].

         for n from 2 to N+1 do

                   for m from 2 to N+1 do

                    JJ[n, m] := diff(ff[n], x[m]);

                   end do;

         end do; – определение матрицы Якоби для итерационной функции

         for l from 0 to N+1 do

                   x[l] := evalf(10^(-8));

         end do; – координаты начального вектора приближений

                   h := 0; – номер итерации

                   r := 10*epsilon;

while r > epsilon do

XX0 := matrix(N, 1, [seq(x[l], l = 1 .. N)]); – начальный вектор приближений

FF0 := matrix(N, 1, [seq(ff[k], k = 1 .. N)]); – начальный вектор свободных членов

JJ0 := matrix([seq([seq(JJ[n, m], m = 2 .. N+1)], n = 2 .. N+1)]); – начальная матрица Якоби

is(norm(JJ0) <> 0);

XX1:= evalm(XX0-JJ0^(-1)*FF0); – основная формула метода Ньютона

r := evalf(norm(evalm(XX1-XX0))); – проверка условия сходимости метода Ньютона

         h := h+1;

         print('%-итерация');

         for l from 1 to N do

                   x[l] := XX1[l, 1];

         end do;

end do; – модифицированный метод Ньютона 

for l from 1 to N do

         y[l]:=(x[l]-x[l-1])/tau;

         A := [op(A), [x[l], y[l]]];

         B := [op(B), [l*tau, x[l]]]; 

end do; – создание массивов данных

return A, B;

end proc;

Результат работы процедуры EFDS() изображен на рис. 9.

 

Рисунок 9. Численное решение по неявной схеме [6], полученное с помощью метода Ньютона

 

Заключение

В статье дано подробное описание процедур пользовательской библиотеке VOFDDE 1.0, написанной на языке Maple для количественного и качественного анализа дробного осциллятора Дуффинга (1). Приведены инструменты визуализации и вывода результатов моделирования. Программа VOFDDE 1.0 была зарегистрирована в Роспатенте.

Работа выполнена при финансовой поддержке Гранта Президента РФ МД-758.2022.1.1.

 

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

1. Дьяконов В.П. Maple в математических расчетах. М.: ДМК Пресс, 2011.  800 с.
2. Kim V.A., Parovik R.I. Mathematical model of fractional Duffing oscillator with variable memory // Mathematics. 2020.  vol. 8, no. 11. pp. 20 – 34.
3. Kim V.A., Parovik R.I. Some Aspects of the Numerical Analysis of a Fractional Duffing Oscillator with a Fractional Variable Order Derivative of the Riemann-Liouville Type // AIP Conference Proceedings. 2021. vol. 8, no. 11. pp. 20 – 34.
4. Kim V.A. Duffing oscillator with external harmonic action and variable fractional Riemann-Liouville derivative characterizing viscous friction // Bulletin KRASEC. Physical and Mathematical Sciences. 2016. vol. 13, no. 2. pp. 46–49.
5. KimV.A., Parovik R.I. Application of the Explicit Euler Method for Numerical Analysis of a Nonlinear Fractional Oscillation Equation  // Fractal Fract. 2022. vol. 6,  no. 5. pp. 274 – 293.
6. Ким В.А., Паровик Р.И. Неявная конечно-разностная схема для осциллятора Дуффинга с производной переменного дробного порядка типа Римана-Лиувилля // Вестник КРАУНЦ. Физико-математические науки. 2022. Т. 40. № 3. С. 179-198.