A Study of the Efficiency of Parallel Computing for Constructing Bifurcation Diagrams of the Fractional Selkov Oscillator with Variable Coefficients and Memory

1. Introduction

Constructing bifurcation diagrams for various dynamic systems is of great practical importance, as it allows us to identify various regimes depending on changes in key parameters [1,2,3]. One such dynamic system, in our opinion, is the Selkov oscillator due to its important application in applied problems, for example, in biology for describing glycolytic reactions [4,5] or in seismology for describing microseisms [6]. From a mathematical perspective, the Selkov oscillator is a dynamic system consisting of two nonlinear first-order ordinary differential equations that describes various oscillatory regimes, including self-oscillations.

A generalization of the Selkov oscillator to the case of memory effects is the fractional Selkov oscillator [7]. Memory effects here indicate that the current state of the system depends on previous states, i.e., on its prehistory [8,9]. These effects can be described using the theory of fractional calculus using fractional derivatives [10,11,12]. In this case, the Selkov fractional oscillator is a nonlinear system of two ordinary differential equations of fractional orders. The introduction of fractional derivatives also entails the introduction of an additional parameter—a characteristic time scale—to match the dimensions of the right- and left-hand sides of the model equations, which is important in the study of dynamic regimes [13].

In the author’s work [7,14], a quantitative and qualitative analysis of the Selkov fractional oscillator was conducted in cases where the orders of the fractional derivatives were constant and where the characteristic time scale was chosen to be unity. Here, the system’s resting points were studied, the Adams-Bashforth-Multon numerical algorithm was implemented to construct oscillograms and phase trajectories, and regular and chaotic regimes were also investigated using the maximum Lyapunov exponents [34].

Further development of the Selkov fractional oscillator involved introducing derivatives of fractional order variables [16,17], which also resulted in non-constant coefficients in the model equations [18,19,20]. Quantitative and qualitative analyses were conducted using an adapted Adams-Bashforth-Multon algorithm, various oscillatory modes were investigated, and 2D and 3D bifurcation diagrams were constructed. These algorithms were implemented in the ABMSelkovFracSim 2.0 software package using the Python programming language and the PyCharm environment [21,22].

Constructing bifurcation diagrams is a rather computationally intensive task, as it requires repeatedly solving a system of fractional differential equations numerically. This leads to significant computational time expenditures [19,20]. Therefore, a parallel version of the bifurcation diagram calculation, depending on the characteristic time scale, was developed using Python tools . However, the parallel computation algorithm has not been studied for its efficiency depending on CPU threads. Therefore, in this paper, we examine its efficiency by addressing this gap and taking into account the available number of CPU threads, similar to the work [23].

The research plan in this article is structured as follows. In Section 2, we introduce some background information and define the key concepts used in the article. Section 3 presents the problem statement. Section 4 presents the Adams-Bashforth-Multon numerical algorithm for solving the problem. Section 5 describes a parallel algorithm implemented in the ABMSelkovFracSim 2.0 software suite for calculating bifurcation diagrams. Section 6 describes the parallel algorithm for constructing bifurcation diagrams. Section 7 analyzes the parallel algorithm’s performance depending on the number of CPU threads. Finally, we present conclusions on the parallel algorithm’s performance analysis.

2. Preliminaries

Definition 1.

The Gerasimov–Caputo fractional derivative of variable order $0<\alpha \left(t\right)<1$ of a function $x\left(t\right)\in C^1\left[0,T\right]$ has the form [17]:

$${\partial }_{0t}^{\alpha \left(t\right)}x\left(t\right)={\displaystyle \frac{1}{\Gamma (1-\alpha \left(t\right))}}\underset{0}{\overset{t}{\int }}{\displaystyle \frac{\dot{x}\left(\tau \right)d\tau }{{(t-\tau )}^{\alpha \left(t\right)}}},$$

(1)

where $\alpha \left(t\right)$ is a function from the class $C\left[0,T\right]$, $\Gamma \left(·\right)$ is the gamma function.

Remark 1.

In the case where the variable order $\alpha \left(t\right)$ is constant, we arrive at the Gerasimov-Caputo fractional derivative, which has been studied quite well [24,25].

Remark 2.

Here we will not dwell on the properties of the fractional derivative (1), but we will point out that they can be studied in the review articles [16,17], and also refer to the relevant literature sources.

Definition 2.

A bifurcation diagram is a graphical representation of changes in the structure of solutions of a dynamic system when the parameters change. It shows how the stable and unstable states of the system change depending on the values of the parameters.

3. Statement of the Problem

Consider the following dynamic system:

$$\left\{\begin{array}{cc}& {\partial }_{0t}^{{\alpha }_1\left(t\right)}x\left(t\right)=-v_1\left(t\right)x\left(t\right)+w_1\left(t\right)y\left(t\right)+h_1\left(t\right)x^2\left(t\right)y\left(t\right),x\left(0\right)=x_0,\hfill \\ & {\partial }_{0t}^{{\alpha }_2\left(t\right)}y\left(t\right)=v_2\left(t\right)-w_2\left(t\right)y\left(t\right)-h_2\left(t\right)x^2\left(t\right)y\left(t\right),y\left(0\right)=y_0.\hfill \end{array}\right.$$

(2)

where $x\left(t\right),y\left(t\right)\in C^1\left[0,T\right]$—solution functions; $v_1\left(t\right)={\theta }^{1-{\alpha }_1\left(t\right)}$, $v_2\left(t\right)=v_0{\theta }^{1-{\alpha }_2\left(t\right)}$, $w_1\left(t\right)=w_0{\theta }^{1-{\alpha }_1\left(t\right)},w_2\left(t\right)=w_0{\theta }^{1-{\alpha }_2\left(t\right)},h_1\left(t\right)=h_0{\theta }^{1-{\alpha }_1\left(t\right)},h_2\left(t\right)=h_0{\theta }^{1-{\alpha }_2\left(t\right)}$—functions from class $C\left[0,T\right]$; $\theta$—parameter with time dimension; $v_0,w_0,h_0$—given constants; $t\in \left[0,T\right]$—current process time; $T>0$—simulation time; $x_0,y_0$—positive constants responsible for the initial conditions.

Operators of fractional variables of orders $0<{\alpha }_1\left(t\right),{\alpha }_2\left(t\right)<1$ in the dynamic system (2) are understood in the sense of Gerasimov-Caputo (1). Following the work [19], we choose the following four types of functions ${\alpha }_1\left(t\right)$ and ${\alpha }_2\left(t\right)$.

1.

Constant fractional orders

$${\alpha }_1\left(t\right)=a_1=\mathrm{const},{\alpha }_2\left(t\right)=a_2=\mathrm{const},$$

where $a_1,a_2\in (0,1]$ are constant values of the fractional orders. Note that when $a_1=a_2=1$, we arrive at the classical Selkov oscillator proposed in the work [4].

2.

Cosine dependence

$${\alpha }_1\left(t\right)=a_1-k_1\left|cos({\varphi }_{1}t+f_1)\right|,{\alpha }_2\left(t\right)=a_2-k_2\left|cos({\varphi }_{2}t+f_2)\right|,$$

where $a_1,a_2\in (0,1]$ are the initial order values; $k_1,k_2\in {\mathbb{R}}^{+}$ are the oscillation amplitudes; ${\varphi }_1,{\varphi }_2\in {\mathbb{R}}^{+}$ are the oscillation frequencies; $f_1,f_2\in \mathbb{R}$ are the initial phases.

3.

Exponential dependence

$${\alpha }_1\left(t\right)=a_1-k_{1}exp(-{\varphi }_{1}t),{\alpha }_2\left(t\right)=a_2-k_{2}exp(-{\varphi }_{2}t),$$

where $a_1,a_2\in (0,1]$ are the limiting values as $t\to \infty$, $k_1,k_2\in {\mathbb{R}}^{+}$ are the initial deviations, ${\varphi }_1,{\varphi }_2\in {\mathbb{R}}^{+}$ are the rates of exponential decay.

4.

Linear dependence

$${\alpha }_1\left(t\right)=max(a_1-k_{1}t,\epsilon ),{\alpha }_2\left(t\right)=max(a_2-k_{2}t,\epsilon ),$$

where $a_1,a_2\in (0,1]$ are the initial values at $t=0$, $k_1,k_2\in {\mathbb{R}}^{+}$ are the coefficients of linear decrease, $\epsilon =0.01$ is the lower boundary for computational stability.

Definition 3.

The dynamic system (2) will be called a Selkov fractional oscillator with variable coefficients and memory, or simply a Selkov fractional oscillator (SFO).

Remark 3.

It should be noted that in the works of the author [26,27] a special case of the SFO (2) was investigated for $\theta =1$.

However, the characteristic time scale $\theta$, as further studies of bifurcation diagrams [18,19,20] showed, plays an important role in the dynamics of the SFO. Therefore, as test examples, we will use the calculation of bifurcation diagrams depending on the values of the parameter $\theta$.

4. Adams–Bashforth–Multon Method

To study the SFO (2), we use the Adams–Bashforth–Moulton numerical method (ABM) from the family of predictor–corrector methods. The ABM has been studied and discussed in detail in [28,29,30,31]. We adapt this method to solve the SFO (2). To do this on a uniform grid N with step $\tau =T/N$, we introduce the functions $x_{k+1}^p,y_{k+1}^p,k=0,\dots ,N-1$, which will be determined by the Adams–Bashforth formula (predictor):

$$\left\{\begin{array}{c}{x}_{k+1}^p=x_0+{\displaystyle \frac{{\tau }^{{\alpha }_{1,k}}}{\Gamma \left({\alpha }_{1,k}+1\right)}}{\displaystyle \sum _{j=0}^k}{\theta }_{j,k+1}^1\left(-v_{1,j}{x}_j+w_{1,j}{y}_j+h_{1,j}{x}_j^2{y}_j\right),\hfill \\ y_{k+1}^p=y_0+{\displaystyle \frac{{\tau }^{{\alpha }_{2,k}}}{\Gamma \left({\alpha }_{2,k}+1\right)}}{\displaystyle \sum _{j=0}^k}{\theta }_{j,k+1}^2\left(v_{2,j}-w_{2,j}{y}_j-h_{2,j}{x}_j^2{y}_j\right),\hfill \\ {\theta }_{j,k+1}^i={\left(k-j+1\right)}^{{\alpha }_{i,k}}-{\left(k-j\right)}^{{\alpha }_{i,k}},i=1,2.\hfill \end{array}\right.$$

(3)

For the corrector (Adams–Moulton formula), we obtain

$$\left\{\begin{array}{cc}\hfill & x_{k+1}=x_0+K_{1,k}\left(-v_{1,k+1}{x}_{k+1}^p+w_{1,k+1}{y}_{k+1}^p+h_{1,k+1}{x}_{k+1}^{p2}{y}_{k+1}^p\right)+\hfill \\ \hfill & +K_{1.k}\left({\displaystyle \sum _{j=0}^k}{\rho }_{j,k+1}^1\left(-v_{1,j}{x}_j+w_{1,j}{y}_j+h_{1,j}{x}_j^2{y}_j\right)\right),\hfill \\ \hfill & y_{k+1}=y_0+K_{2,k}\left(v_{2,k+1}-w_{2,k+1}{y}_{k+1}^p-h_{2,k+1}{x}_{k+1}^{p2}{y}_{k+1}^p\right)+\hfill \\ \hfill & +K_{2,k}{\displaystyle \sum _{j=0}^k}{\rho }_{j,k+1}^2\left(v_{2,j}-w_{2,j}{y}_j-h_{2,j}{x}_j^2{y}_j\right).\hfill \end{array}\right.$$

(4)

where $K_{1,k}={\displaystyle \frac{{\tau }^{{\alpha }_{1,k}}}{\Gamma \left({\alpha }_{1,k}+2\right)}},K_{2,k}={\displaystyle \frac{{\tau }^{{\alpha }_{2,k}}}{\Gamma \left({\alpha }_{2,k}+2\right)}}$, and weight coefficients in (3) are determined by the following formula:

$${\rho }_{j,k+1}^i=\left\{\begin{array}{cc}& k^{{\alpha }_{i,k}+1}-\left(k-{\alpha }_{i,k}\right){\left(k+1\right)}^{{\alpha }_{i,k}},j=0,\hfill \\ & {\left(k-j+2\right)}^{{\alpha }_{i,k}+1}+{\left(k-j\right)}^{{\alpha }_{i,k}+1}-2{\left(k-j+1\right)}^{{\alpha }_{i,k}+1},1\le j\le k,\hfill \\ & 1,j=k+1,\hfill \\ & i=1,2.\hfill \end{array}\right.$$

Formulas (3) and (4) form the basis of a modified ABM method adapted for a fractional-order system (SFO) with variable exponents (${\alpha }_{1,k},{\alpha }_{2,k}$) and coefficients. The advantage of the ABM method is that it does not require solving nonlinear systems of equations at each step, significantly reducing computational costs. However, the main computational costs here are associated with storing the entire history (memory) of the vectors $x_j,y_j$, which is typical for fractional integro-differentiation methods.

Remark 4.

The presented ABM method is explicit in its implementation, despite the use of the implicit Adams-Multon formula (4). Therefore, the ABM method is conditionally stable, like most explicit methods. It should also be noted that the order of accuracy of the ABM method is higher than that of the nonlocal explicit finite-difference scheme. The properties of the numerical ABM method (3), (4), convergence and stability issues can be studied in the work [18,19,20].

Based on the above Remark 4, the functions ${\alpha }_1\left(t\right)$ in the test examples for this article were chosen to be monotonically decreasing or not rapidly changing on the interval (0,1), and the sampling step $\tau$ was chosen small enough to prevent the solution from exhibiting unphysical oscillations or error growth.

The numerical algorithm ABM (3), (4) was implemented in the ABMSelkovFracSim 2.0 software package for four types of ${\alpha }_1\left(t\right),{\alpha }_2\left(t\right)$ functions. We will refer to them as methods: ABMSelkovFracConst, ABMSelkovFracCos, ABMSelkovFracExp, and ABMSelkovFracLine.

5. ABMSelkovFracSim 2.0 Software Package

The paper [20] provides a detailed description of the ABM software package. We will focus on the module for calculating bifurcation diagrams (the Bifurcation analysis mode, Figure 1).

Figure 1 shows a screenshot of the main window of the Bifurcation Analysis mode. Here we see that this mode’s interface not only allows you to enter values for key FOS parameters and select one of four calculation methods depending on the function types ${\alpha }_1\left(t\right),{\alpha }_2\left(t\right)$, but also to set the range and step size for the key parameter $\theta$. Furthermore, the required number of CPU threads can be selected, which are automatically determined for efficient parallel calculation of the bifurcation diagram.

Also, on the left side of the ABMSelkovFracSim 2.0 software interface in the Bifurcation Analysis mode, we see visualizations of the calculation results in the form of bifurcation diagrams and graphs of the functions $x\left(\theta \right),y\left(\theta \right)$. In Figure 1, the bifurcation diagram shows two ranges of parameter $\theta$ values, within which different dynamic FOS modes exist [19].

Remark 5.

It should be noted that the following useful information is displayed in the interface for the user of the ABMSelkovFracSim 2.0 software package: the time of calculation of the bifurcation diagram, the number of processed points in dynamics and the progress bar of the algorithm’s operation, and the possible presence of errors in the calculation.

Next, we will discuss in more detail the implementation of the parallel part of the bifurcation diagram calculation algorithm in the ABMSelkovFracSim 2.0 software.

6. Parallel Implementation of the Bifurcation Diagram Calculation Algorithm

Parallel implementation of algorithms is an important and relevant task in reducing their execution time. There are many approaches and technologies for parallelizing algorithms; these can be studied, for example, in the articles [32,33]. In this section, we will describe a parallelization method based on the built-in libraries of the Python programming language [21].

Due to the fact that the numerical algorithm ABM (3), (4) is practically impossible to parallelize due to the existence of true dependencies between the grid functions of the solutions $x_k$ and $y_k$, we will parallelize multiple calculations using the ABM algorithm for each fixed value of the parameter $\theta$ from a given interval. That is, the operating principle of our parallel algorithm consists of the following steps: a separate independent process is created for each value of the parameter $\theta$; processes are executed in parallel on multiple processor cores; the results are collected after all calculations are completed. This is the classic "MapReduce" approach, where Map is the distribution of calculations across processes, and Reduce is the collection and merging of results. This approach is optimal for parametric studies, where each calculation is independent of the others.

To implement this approach, we will use Native Python Multiprocessing (Python Multiprocessing, ProcessPoolExecutor) technology from the Python Standard Library (Algorithm 1).

Algorithm 1 Bifurcation analysis with parallel computing

1: Step 1: The user runs the bifurcation analysis 2: Step 2: Create a list of values $\theta =[{\theta }_1,{\theta }_2,\dots ,{\theta }_n]$ 3: for $i=1$ to n do 4:     Create a task $T_i$ to evaluate for $\theta ={\theta }_i$ 5: end for 6: Step 4: Initialize ProcessPoolExecutor with N processes 7: Step 5: Assign tasks $\{T_1,T_2,\dots ,T_n\}$ to processes 8: Step 6:                          ▹ Parallel execution 9: for $k=1$ to N in parallel do 10:     Independently compute results for assigned tasks 11: end for 12: Step 7: Collect all results from processes 13: Step 8: Combine results into a single data structure 14: Step 9: Build a common bifurcation graph

Remark 6.

Let’s highlight the key features of Algorithm 1:

• Data Parallelism: Each value of ${\theta }_i$ is processed independently.
• Process-based Parallelism: Processes (not threads) are used to bypass the Global Interpreter Lock (GIL).
• Master-Worker Pattern: The master process distributes tasks, and workers perform computations.
• Fault Tolerance: An error in one process does not affect others.
• Dynamic Load Balancing: ProcessPoolExecutor automatically distributes tasks.

We do not use other CPU parallelization technologies here for the following reasons. OpenMP is used for parallelization within a single process [34], which is not possible in our case. MPI is used for distributed computing on clusters, which is also not necessary in our case [35]. Threading: Threads in Python have a GIL and are not suitable for numerical calculations [36].

We now turn to the main goal of our article, to studying the efficiency of a parallel algorithm for calculating bifurcation diagrams.

7. Analysis of Algorithm Efficiency Based on Average Execution Time

All calculations in this article were performed on a computer with the following specifications: CPU — AMD Ryzen 7 7800X3D $16\times 4.2$ GHz cores; L2 cache 12 Mb & L3 cache 128 Mb; RAM — 64 Gb; GPU — NVIDIA GeForce RTX 5070 Ti, 16 Gb, 2300 MHz, ALU 8960.

Let us analyze the efficiency of the parallel computation algorithm of the ABMSelkovFracSim 2.0 software package in the mode of calculating and constructing a bifurcation diagram (Bifurcation analysis).

Since each new run of Bifurcation analysis will yield slightly different results for T — execution time in [sec.], and the number of experiments is finite, T can be considered a discrete random variable with some distribution function and mathematical expectation:

$$T_p\left(N\right)={\displaystyle \frac{1}{L}}{\displaystyle \sum _{i=1}^L}{T}_p^i\left(N\right),$$

(5)

where i is the index of the Bifurcation analysis numerical experiment; L is the sample size; p is the number of CPU threads; N is the number of nodes in the uniform grid of the numerical method, i.e., the size of the input data.

Algorithm efficiency is the optimal ratio between computational speedup and the number of CPU threads, compared to the most efficient sequential version of the algorithm. However, when using a large number of CPU threads, the hardware and software limitations of parallelization’s acceleration efficiency become increasingly apparent. These limitations include the high time costs of parallelizing the task and reassembling all computational results [37], known as Amdahl’s law [23]. To represent various data on average execution time, we will use the following parameters in [sec]:

• $T_1\left(N\right)$ — the execution time of a test example of size N required by the sequential Bifurcation analysis algorithm;
• $T_p\left(N\right)$ — the execution time of a test example of size N required by the parallel Bifurcation analysis algorithm on a machine with $p\ge 2$ CPU threads.

Let us calculate $T_p\left(N\right)$ for different numbers of used CPU threads p in [units] with a step of 2.

To obtain efficiency estimates, data on average execution time are considered in terms of (TAECO) [38], applicable to algorithms: T (execution time), A (acceleration, speedup), E (efficiency), C (cost), O (cost-optimal indicator).

• $A_p\left(N\right)$ — is the speedup of the algorithm, in [units], provided by the parallel version of the algorithm compared to the sequential one, and is calculated as follows:

  $$A_p\left(N\right)={\displaystyle \frac{T_1\left(N\right)}{T_p\left(N\right)}}.$$

  (6)
• $E_p\left(N\right)$ — is the efficiency of the algorithm, in [units/thread], of using a given number p of CPU threads and is defined by the ratio:

  $$E_p\left(N\right)={\displaystyle \frac{T_1\left(N\right)}{p·T_p\left(N\right)}}={\displaystyle \frac{A_p\left(N\right)}{p}},$$

  (7)

  moreover, by definition, the sequential algorithm has the greatest efficiency, $E_1\left(N\right)=1$. Therefore, for a parallel algorithm, it is better if $E_p\left(N\right)\to 1$.
• $C_p\left(N\right)$ — is the cost of the algorithm, in [sec. × thread], which is determined by the product of a given number p of CPU threads and the T execution time of the parallel algorithm. The cost is determined by the ratio:

  $$C_p\left(N\right)=p·T_p\left(N\right).$$

  (8)
• $O_p\left(N\right)$ — is the cost-optimal indicator (COI) of the algorithm, in [units × thread], characterized by a cost proportional to the complexity of the most efficient sequential algorithm:

  $$O_p\left(N\right)={\displaystyle \frac{C_p\left(N\right)}{T_1\left(N\right)}},$$

  (9)

  and the closer the value is to 1, the better, i.e., the cheaper the use of the parallel algorithm in terms of engaging CPU threads.

The calculations of efficiency indicators (5)-(9) and their visualization were performed in the MatlabR2025b computer mathematics system.

Example 1.

For the first example, we will choose the ABMSelkovFracCos method and the following parameter values for FOS (2): $T=200,N=3000,x_0=1,y_0=0.5,v_0=0.6,w_0=0.03,h_0=1.3$, $a_1=0.8,a_0=0.9,k_1=0.04,k_2=0.02,{\varphi }_1=3.5,{\varphi }_1=1.5,f_1=f_2=0$.

$${\alpha }_1\left(t\right)=0.8-0.04cos\left(3.5t\right),{\alpha }_2\left(t\right)=0.9-0.02cos\left(1.5t\right).$$

In the Bifurcation analysis mode with all the same parameters, calculation method ABMSelkovFracCos, but with $\theta \in [0.1,2.1]$ and step $h=0.1$. To obtain the average execution time according to (5), we define $L=10$, which will allow for smoother graphs for both average time and efficiency analysis.

The results are shown in Table 1 and Figure 2 below.

Figure 2 presents the results of the efficiency analysis of the parallel Bifurcation analysis algorithm for the ABMSelkovFracCos method compared to its sequential counterpart. The results are presented as histograms for greater clarity. Parallelization, specific to this computer, allows achieving a maximum speedup of up to 8.5 times.

Figure 2 shows a dependence close to exponential: average execution time vs. number of CPU threads. Increasing the number of CPU threads at $p\le 8$ gives a greater performance gain compared to the sequential algorithm, but increasing the number of CPU threads at $p>8$ yields less and less performance gain. In the COI estimate, we see the opposite trend: the cost of using parallel algorithms does not change until around $p=\pm 7$, and then increases linearly and proportionally to the growth in the number of CPU threads. Consequently, at $p=8$ we have the best efficiency coefficient of the parallel algorithm $E_8\left(3000\right)=0.844$ in relation to the cost of using the parallel algorithm $O_8\left(3000\right)=1.183\to 1$, while obtaining a speedup of 6.75 times. As a result, we can say that the Bifurcation analysis mode is optimally parallelized on 7-8 threads, which, depending on the maximum number of CPUs of a particular computer, will allow running several copies of the ABMSelkovFracSim 2.0 software package itself for simultaneous work in several windows.

Example 2.

To evaluate the efficiency of the parallel version of the calculation algorithm, we will choose the ABMSelkovFracExp method.

$${\alpha }_1\left(t\right)=0.8-0.04exp(-3.5t),{\alpha }_2\left(t\right)=0.9-0.02exp(-1.5t).$$

We will select the parameter values from Example 1. To obtain the average time (5), the sample size is $L=10$.

The results are shown in Table 2 and Figure 3 below.

In the case of Test Example 2 for the ABMSelkovFracExp method, parallelization, for this particular computer, allows achieving a maximum speedup of $\approx 9$ times. The estimates for (Figure 3) have the same trends as for Test Example 1 (Figure 2). However, for $p\le 5$, the cost of using parallel algorithms $O_{p\to 5}\left(3000\right)\approx 1$, and therefore $E_{p\to 5}\left(3000\right)\to 1$, i.e., almost as effective as the sequential algorithm, with a 5-fold performance increase.

Example 3.

Let’s consider a parallel version of the algorithm for the ABMSelkovFracLinear method. For simplicity, we set $k_1=k_2=0.001$.

$${\alpha }_1\left(t\right)=max(0.8-0.001t,0.01),{\alpha }_2\left(t\right)=max(0.9-0.001t,0.01),$$

The remaining parameter values are left unchange.

The results are shown in Table 3 and Figure 4 below.

For Test Example 3 for ABMSelkovFracLinear, parallelization allows for a maximum speedup of $\approx 8.5$ times. The efficiency estimates (Figure 4) show trends similar to those for Test Example 2 (Figure 2).

8. Conclusions

This article examines the effectiveness of a parallel algorithm for calculating bifurcation diagrams for the Selkov fractional oscillator, implemented using the standard Python library. The TAECO approach was chosen to evaluate the effectiveness of the parallel algorithm: T (execution time), A (acceleration), E (efficiency), C (cost), O (cost optimality index). Performance graphs were plotted for three ABM calculation methods: ABMSelkovFracCos, ABMSelkovFracExp, and ABMSelkovFracLinear. For these methods and a specific computer with 16 threads, approximately the same results were obtained. The effectiveness of the parallel algorithm is 8-9 times faster than that of the sequential version, while the optimal number of threads is 8 units. Therefore, the implementation of a parallel algorithm for calculating bifurcation diagrams of the fractional oscillator in the ABM software package is justified. Further development of research on this topic may be associated with the use of heterogeneous parallel programming structures, for example, CPU-GPU architecture, similar to the works [39,40].