Neural network modeling of chaos in a nonlinear magnetic-vortex dynamo system

Neural network modeling of chaos in a nonlinear magnetic-vortex dynamo system M. I. Kopp1 December 8, 2025 1Institute for Single Crystals, NAS Ukraine, Nauky Ave. 60, Kharkiv 61072, Ukraine Abstract This study presents the modeling of chaotic behavior in a nonlinear dynamo system using a feedforward artificial neural network (FFNN). The investigated system was obtained based on asymptotic multiscale nonlinear theory applied to temperature-stratified conducting media with small-scale helical external forcing. This small-scale forcing acts as a source of small-scale helical turbulence with a low Reynolds number. The governing equations form a self-consistent nonlinear magneto-vortex dynamo system consisting of four coupled nonlinear differential equations that describe the evolution of large-scale vortex and magnetic field structures. The existence of both regular and chaotic large-scale fields in the stationary regime has been previously established. In this work, the 4th-5th order Runge-Kutta method was employed for the numerical solution of the system, and the generated data were used to train the neural network. The implementation was performed in MATLAB R2017a without using the Neural Network Toolbox. The obtained results demonstrate high accuracy in predicting the system dynamics, achieving a mean squared error (MSE) on the order of 10−6. This approach provides an efficient computational framework for analyzing complex nonlinear dynamo systems and offers insights into the predictability of chaotic regimes in magnetohydrodynamic phenomena. Keywords: magneto-vortex dynamo; chaotic behavior; Runge-Kutta method; artificial neural network 1 Introduction Chaotic dynamo systems are of particular interest in geophysics and astrophysics due to their ability to model magnetic field generation processes in natural objects. Among the classical dynamo systems in which chaos has been investigated, one can distinguish the Rikitake system for modeling geomagnetic dynamo [1] and the homopolar (single-disk) dynamo of Bullard [2]. The Rikitake system, consisting of two disks, demonstrates chaotic behavior and has been successfully applied to describe geomagnetic field reversals of the Earth [3, 4]. Bullard’s homopolar dynamo, despite its simplicity, is also capable of demonstrating complex dynamics 1 under certain parameters [5]. In contrast to these models, the dynamo system [6] investigated in this work was derived from the full averaged equations of magnetohydrodynamics using the method of asymptotic expansions without additional closure hypotheses. This provides a more rigorous physical foundation for the model and allows describing the coupled evolution of largescale vortex and magnetic fields in a temperature-stratified conducting medium. The system is characterized by four nonlinear differential equations with variables, describing the dynamics of vortex and magnetic structures. In recent years, artificial neural networks (ANNs) have been actively applied for modeling and predicting the behavior of chaotic systems. The universal approximation theorem [7, 8] provides the theoretical foundation for using neural networks to approximate complex nonlinear functions. Lamamra et al. [9] used a multilayer perceptron (MLP) with structure optimization using the NSGA-II genetic algorithm for modeling the chaotic Mackey-Glass system, achieving a mean squared error (MSE) of the order of 10−13. Their work demonstrated that optimal neural network topology can be determined through multi-objective optimization, balancing the number of neurons against prediction accuracy. Cimen et al. [10] demonstrated the possibility of modeling chaotic motion of a modified Lorenz system using NAR-type neural networks trained on data obtained through image processing methods. This innovative approach combined computer vision techniques with neural network prediction, opening new possibilities for modeling physically observable chaotic systems. Al-Musawi et al. [11] implemented prediction of the chaotic Chua system using a feedforward neural network on FPGA, utilizing the 32-bit IEEE-754 floating-point format, and applied the developed model for information encryption. Their work achieved an MSE on the order of 10−6and demonstrated the practical applicability of hardware-implemented neural networks for chaos-based cryptographic applications. The Chua circuit, being one of the simplest electronic circuits exhibiting chaotic behavior [12], has found numerous applications in secure communications [13]. Keles et al. [14] investigated the modeling of the chaotic Rucklidge system using various ANN architectures (FFNN, LRNN, CFNN) and different activation functions, determining the most effective configurations for predicting a system with three variables. The Rucklidge system, which models double convection in a rotating fluid layer with an applied vertical magnetic field [15], demonstrated that Vanilla RNN with a sliding window approach achieved an R2score of 0.999999 and RMSE of 1.825 ×10−4. Ramachandruni et al. [16] conducted a comprehensive study on the application of machine learning and neural networks for analyzing and predicting chaos in multi-pendulum systems, comparing 10 different models and neural network architectures. The best results were shown by the Long Short Term Memory (LSTM) network with RMSE of the order of 10−2. Their novel time-step based approach demonstrated the ability to predict chaotic motion for completely unseen initial conditions, which is crucial for real-world applications. Neural networks have also been successfully applied to other chaotic systems. Alcin et al. [17] implemented an ANN-based Pehlivan-Uyaroglu chaotic system in VHDL using the IEEE754 32-bit floating-point standard. Koyuncu et al. [18] developed a multi-layer feedforward ANN on FPGA for a novel four-dimensional hyper-chaotic system, demonstrating the effectiveness of hardware implementations for real-time chaos-based applications. Zhang and Lei [19] examined the training efficiency of multilayer ANN architectures for chaotic systems implementation, providing insights into optimal network depth and width for different types of 2 chaotic dynamics. The application of recurrent neural networks to chaotic time series prediction has been extensively studied. LSTM networks, introduced by Hochreiter and Schmidhuber [20], have shown particular promise due to their ability to capture long-term dependencies in sequential data. Panahi et al. [21] used chaotic artificial neural networks for modeling epilepsy, demonstrating the medical applications of chaos theory. Chattopadhyay et al. [22] employed data-driven predictions using reservoir computing, ANNs, and LSTM networks for the multiscale Lorenz 96 chaotic system, comparing the effectiveness of different architectures. For dynamo systems specifically, adaptive control and synchronization methods have been developed. Vaidyanathan [3] presented adaptive synchronization techniques for the Rikitake two-disk dynamo system, while chaos control methods have been implemented for Tokamak systems with magnetically confined plasma [23]. However, despite successful applications of ANNs to various chaotic systems including the classical attractors (Lorenz, Rossler, Chua) and mechanical systems (pendulums, oscillators), there is a notable absence in the literature of works on modeling complex dynamo systems derived from magnetohydrodynamics equations using artificial neural networks. The present study aims to fill this gap and demonstrate the effectiveness of applying neural network approaches to predicting the behavior of a physical dynamo system [6] with four coupled variables. In our recent work [24], we reported on the prediction of chaotic behavior of largescale vortex and magnetic fields by training a feedforward neural network (FFNN) on time series data obtained by numerically integrating dynamo differential equations [6] using the odeint solver from SciPy and NumPy in Python. The objective of this work is to develop and train a feedforward artificial neural network (FFNN) in MATLAB R2017a for modeling the chaotic dynamics of a nonlinear dynamo system describing the generation of large-scale magnetic and vortex fields. The 4th-5th order Runge-Kutta method is used for numerical solution of the system, and the obtained data are applied for training the neural network, implemented in MATLAB R2017a without using specialized toolboxes. This approach ensures reproducibility and demonstrates that effective chaotic system modeling can be achieved without reliance on proprietary software packages. 2 Mathematical model of nonlinear dynamo The chaotic dynamo system considered in this study is derived from the nonlinear dynamo theory developed by Kopp et al. [6]. After applying the multiscale asymptotic analysis and quasi-two-dimensional approximation, the system describing the evolution of large-scale magnetic and velocity fields in a stratified conducting medium takes the form of four coupled nonlinear ordinary differential equations:                      dx dt = +α2(p, P)·p+c1, dp dt =−α1(x, X)·x+c2, dX dt =−β2(p, P)·P+c3, dP dt = +β1(x, X)·X+c4, (1) 3 where x(t) and p(t) represent the large-scale velocity field components, X(t) and P(t) represent the large-scale magnetic field components, and c1= 0.01, c2= 0.01, c3= 0.001, c4= 0.001 are constant terms arising from the integration of the non-stationary equations [6]. In the system of equations (1), the coefficients α1,α2,β1, and β2are not constants but complex nonlinear functions of the system state and physical parameters. These coefficients represent the combined effects of hydrodynamic and magnetohydrodynamic α-effects in the presence of thermal stratification. The hydrodynamic α-effect coefficients are given by: α1= R(1 + P2 mx2) [(1 + Pr)(1 + P2 mx2) + QX2(Pr−Pm)] 1−QX2Pm 1+P2 mx2 2 [(1 −Pmx2+QX2)2+x2(1 + Pm)2]·d1 , α2= R(1 + P2 mp2) [(1 + Pr)(1 + P2 mp2) + QP2(Pr−Pm)] 1−QP 2Pm 1+P2 mp2 2 [(1 −Pmp2+QP2)2+p2(1 + Pm)2]·d2 ,(2) where the denominators d1and d2are: d1=(1 −Pmx2+QX2)2+x2(1 + Pm)2(1 + P2 rx2)+ +2R(1 −Prx2)(1 + P2 mx2) + QX2(1 + Pmx2)+R2(1 + P2 mx2), d2=(1 −Pmp2+QP2)2+p2(1 + Pm)2(1 + P2 rp2)+ +2R(1 −Prp2)(1 + P2 mp2) + QP2(1 + Pmp2)+R2(1 + P2 mp2). The magnetohydrodynamic α-effect coefficients are: β1=P2 m (1 −Pmx2+QX2)2+x2(1 + Pm)21−R·N1 M1, β2=P2 m (1 −Pmp2+QP2)2+p2(1 + Pm)21−R·N2 M2,(3) with N1= 1−Prx2+QX2(1 + PrPmx2) 1 + P2 mx2+R, M1=[(1 −Pmx2+QX2)2+x2(1 + Pm)2] (1 + P2 rx2) 1 + P2 mx2+ +2R1−Prx2+QX2(1 + PrPmx2) 1 + P2 mx2+R2, N2= 1−Prp2+QP2(1 + PrPmp2) 1 + P2 mp2+R, M2=[(1 −Pmp2+QP2)2+p2(1 + Pm)2] (1 + P2 rp2) 1 + P2 mp2+ +2R1−Prp2+QP2(1 + PrPmp2) 1 + P2 mp2+R2. The system contains several dimensionless parameters that characterize the physical properties of the medium: •R: The Rayleigh number, representing the strength of thermal convection relative to viscous dissipation. This is the key parameter controlling the dynamo instability. 4 •Pm=ν/νm: The magnetic Prandtl number, ratio of kinematic viscosity to magnetic diffusivity. In our simulations, Pm= 1. •Pr=ν/χ: The Prandtl number, ratio of kinematic viscosity to thermal diffusivity. We use Pr= 1 for simplicity. •Q=σB2 0λ2 0/(c2ρ00ν): The Chandrasekhar number, measuring the relative importance of magnetic forces to viscous forces. We set Q= 1. It is easy to see from system (1) that, after grouping the equations into pairs, one pair of oscillators parametrically modulates the effective frequency (or inertial term) of the other. This mechanism generates a two-way interaction between the subsystems, where each pair acts simultaneously as a driver and a responder, shaping the overall dynamics. Such a structure creates genuine bidirectional coupling between the velocity and magnetic components: changes in one subsystem instantaneously alter the restoring forces and effective inertia in the other, resulting in a dynamically adaptive feedback loop. Moreover, the effective frequencies ω1= √α1α2, ω2=√β1β2, are not fixed parameters but depend explicitly on the instantaneous state of the system. This position-dependent frequency modulation produces time-varying resonance conditions and enables parametric excitation. As a result, even initially regular oscillations can undergo strong modulation, giving rise to complex and potentially chaotic behaviour. The dynamical system evolves in a four-dimensional phase space (x, p, X, P) and possesses several geometric features typical of conservative models: it admits fixed points (including the origin), lacks attractors, and is invariant under the inversion (x, p, X, P)→(−x, −p, −X, −P). In such settings, the qualitative behaviour is strongly influenced by the persistence or destruction of invariant structures. According to the Kolmogorov–Arnold–Moser (KAM) theory, weak deviations from integrability preserve most invariant tori, supporting quasi-periodic trajectories, while resonant tori break up and form stochastic layers. In our case the nonlinear coupling terms and the constants Ciact as the effective perturbation, enabling resonances between the two oscillator pairs and opening pathways to chaotic dynamics. The complexity of the coefficient functions further allows for the formation of homoclinic and heteroclinic tangles associated with hyperbolic equilibria. Transverse intersections of their stable and unstable manifolds generate Smale horseshoes, producing sensitive dependence on initial conditions and global chaos in the four-dimensional phase space. From a dynamo-theoretic perspective these chaotic regimes have clear physical meaning. They correspond to intermittent magnetic-field generation, consistent with the irregular bursts and quiet intervals observed in astrophysical dynamos. The control parameter Rgoverns the transition from a laminar regime with coherent magnetic structures to a turbulent-like regime characterised by chaotic reversals. The nonlinear coupling terms mediate energy exchange between the kinetic subsystem (x, p) and the magnetic subsystem (X, P), as well as between different spatial scales through an α-effect mechanism. When this exchange becomes strongly modulated and irregular, the system naturally enters a chaotic state. Such behaviour provides a plausible explanation for irregular field reversals in stars and planets, the spontaneous formation of magnetic structures, and transitions between distinct dynamo regimes. 5 3 Geometric structure and Hamiltonian formalism The nonlinear dynamo system described by equations (1) exhibits a rich geometric structure characteristic of conservative dynamical systems. Despite the presence of constant shift terms c1, c2, c3, c4, the system preserves phase space volume and possesses properties suggestive of underlying Hamiltonian or generalized symplectic structure. 3.1 Conservation of phase space volume A fundamental property of the system (1) is the vanishing divergence of the vector field: ∇·Φ=∂˙x ∂x +∂˙p ∂p +∂˙ X ∂X +∂˙ P ∂P =∂ ∂x[α2(p, P)p+c1] + ∂ ∂p[−α1(x, X)x+c2]+ +∂ ∂X [−β2(p, P)P+c3] + ∂ ∂P [β1(x, X)X+c4] = 0 + 0 + 0 + 0 = 0.(4) This establishes that the system is conservative in the sense of Liouville’s theorem: phase space volumes are preserved under the flow. The constant terms ciact as uniform translations in phase space and do not affect volume preservation. 3.2 Structure matrix and quasi-Hamiltonian form The equations (1) can be written in a structure-matrix formulation: d dt     x p X P     =J(x, p, X, P)·    p x P X     +    c1 c2 c3 c4     ,(5) where the structure matrix is: J(x, p, X, P) =     0α2(p, P) 0 0 −α1(x, X) 0 0 0 0 0 0 −β2(p, P) 0 0 β1(x, X) 0     .(6) The matrix Jis not skew-symmetric in general, since: J21 6=−J12 ⇔α1(x, X)6=α2(p, P).(7) This asymmetry indicates that the system cannot be cast in canonical Hamiltonian form with standard Poisson brackets. However, the block-diagonal structure and volume preservation suggest a generalized geometric framework. 6 3.3 Attempt at Hamiltonian construction For a classical Hamiltonian function H(x, p, X, P) to generate the dynamics via: ˙x=∂H ∂p ,˙p=−∂H ∂x ,˙ X=∂H ∂P ,˙ P=−∂H ∂X ,(8) we would need (ignoring constants momentarily): ∂H ∂p =α2(p, P)p, ∂H ∂x =α1(x, X)x, (9) ∂H ∂P =−β2(p, P)P, ∂H ∂X =−β1(x, X)X. (10) Integrability conditions. For Hto exist, mixed partial derivatives must be equal: ∂2H ∂x∂p =∂ ∂x[α2(p, P)p] = 0,∂2H ∂p∂x =∂ ∂p[α1(x, X)x] = 0.(11) These are trivially satisfied due to variable separation (α2independent of x,α1independent of p). Similarly for the magnetic sector. Integration attempt. From ∂H ∂p =α2(p, P)p, integrating with respect to p: H(x, p, X, P) = Zα2(p, P)p dp +f1(x, X, P),(12) where f1is an integration ”constant”, i.e., function of other variables. From ∂H ∂x =α1(x, X)x: ∂ ∂x Zα2(p, P)p dp +f1(x, X, P)=α1(x, X)x. (13) This simplifies to ∂f1 ∂x =α1(x, X)x . (14) After integrating, we get f1(x, X, P) = Rα1(x, X)x dx +f2(X, P). Continuing this process for the magnetic components and ensuring consistency is nontrivial. The resulting Hamiltonian, if it exists, would have the general form: H(x, p, X, P) = Zα2(p, P)p dp +Zα1(x, X)x dx −Zβ2(p, P)P dP −Zβ1(x, X)X dX+ +cross-terms.(15) However, the cross-terms needed to ensure all consistency conditions are satisfied become exceedingly complex for the specific functional forms in equations (2). 7 3.4 Validation of the generalized Poisson structure for dynamo system The most general geometric framework is that of a Poisson manifold with a state-dependent Poisson tensor Πij(x, p, X, P). Any observable f(x, p, X, P) evolves according to: df dt ={f, g}Π= 4 X i,j=1 Πij(x, p, X, P)∂f ∂zi ∂g ∂zj,(16) where z= (z1, z2, z3, z4) = (x, p, X, P) and the Poisson tensor is: Πij =    0α2(p, P) 0 0 −α1(x, X) 0 0 0 0 0 0 −β2(p, P) 0 0 β1(x, X) 0     .(17) The Poisson tensor must satisfy: •Skew-symmetry: Πij =−Πji •Jacobi identity: 4 X `=1 Πi` ∂Πjk ∂z`+ Πj` ∂Πki ∂z`+ Πk` ∂Πij ∂z`= 0 ∀i, j, k. (18) Verification of Jacobi identity for key components. Consider i= 1 (x), j= 2 (p), k= 3 (X): Π1`∂Π23 ∂z`+ Π2`∂Π31 ∂z`+ Π3`∂Π12 ∂z`= Π12 ∂Π23 ∂p + Π21 ∂Π31 ∂x + Π34 ∂Π12 ∂P = =α2(p, P)·0+(−α1(x, X)) ·0+(−β2(p, P))∂α2(p, P) ∂P =−β2(p, P)∂α2(p, P) ∂P .(19) For the Jacobi identity to hold, we need: β2(p, P)∂α2(p, P) ∂P = 0 or specific relations among αi, βi.(20) This is not automatically satisfied for arbitrary coefficient functions. The Jacobi identity imposes strong constraints on the functional forms of αi(x, X) and βi(p, P). Clearly, when the Jacobi identity fails, the generalized bracket (16) does not define a true Poisson manifold. Instead, the system may exhibit a almost-Poisson or pre-symplectic structure, which are weaker geometric models that still allow some conservation laws and Hamiltonian dynamics [25, 26]. 8 3.5 Nambu mechanics formulation For four-dimensional conservative systems, Nambu mechanics [27] offers a more promising framework that generalizes Hamiltonian dynamics to systems with multiple conserved quantities. For a four-dimensional phase space, the Nambu bracket formulation employs a 4-bracket structure: df dt ={f, C1, C2, C3}N=∂(f, C1, C2, C3) ∂(x, p, X, P)= det      ∂f ∂x ∂f ∂p ∂f ∂X ∂f ∂P ∂C1 ∂x ∂C1 ∂p ∂C1 ∂X ∂C1 ∂P ∂C2 ∂x ∂C2 ∂p ∂C2 ∂X ∂C2 ∂P ∂C3 ∂x ∂C3 ∂p ∂C3 ∂X ∂C3 ∂P      ,(21) where εijkl is the Levi-Civita symbol, C1, C2, C3are three independent conserved quantities (Casimir functions). The equations of motion become: dx dt ={x, C1, C2, C3}N,(22) dp dt ={p, C1, C2, C3}N,(23) dX dt ={X, C1, C2, C3}N,(24) dP dt ={P, C1, C2, C3}N.(25) The fundamental property of the Nambu bracket is that C1, C2, C3are automatically conserved: {Ci, C1, C2, C3}N= 0 for i= 1,2,3. Candidate conserved quantities might include: 1. Generalized energies. C1=Evortex =Zx 0 α1(s, X)s ds +Zp 0 α2(s, P)s ds, (26) C2=Emagnetic =ZX 0 β1(x, s)s ds +ZP 0 β2(p, s)s ds. (27) 2. Enstrophy/magnetic helicity analogs. In fluid and plasma systems, higher-order moments (enstrophy, helicity) are often conserved. For our system: C3=Z[α1(x, X)x2+β1(x, X)X2]w1(x, X)dV +Z[α2(p, P)p2+β2(p, P)P2]w2(p, P)dV, (28) where w1, w2are weight functions ensuring conservation. Given the algebraic complexity of αi, βiin equations (2), analytical derivation of conserved quantities is difficult. The inability to find an explicit Hamiltonian or Nambu formulation analytically doesn’t prevent neural networks from learning the dynamics empirically, demonstrating the power of data-driven approaches for complex physical systems. 4 Numerical solutions of nonlinear dynamo equations Having established the geometric and mathematical foundations of the dynamo system, we now turn to its computational implementation. This section describes the numerical integration 9 Figure 8: Time series of the dynamo system variables for the chaotic regime with initial conditions (x0, p0, X0, P0) = (0.8,0.8,0.01,0.01) at R= 5. Figure 9: Phase space portraits for the chaotic regime at R= 5: (top left) (x, p) projection; (top right) (X, P) projection; (bottom) three-dimensional projections (x, p, X) and (X, P, x). 16 Figure 10: Poincar´e sections for the weakly chaotic regime with initial conditions (0.8,0.8,0.01,0.01) at R= 5. Left: (x, p) section. Right: (X, P) section. Case 2: Initial conditions (0.9,0.9,0.01,0.01) A further slight increase in the initial conditions for the velocity components to x0= 0.9 and p0= 0.9 while maintaining R= 5 results in a transition to a strongly chaotic regime with a sharp increase in the irregularity and amplitude of the magnetic field oscillations, demonstrating the extreme sensitivity to initial conditions characteristic of chaotic systems. The time evolution in Figure 11 shows that, while the velocity components x(t) and p(t) perform weakly chaotic oscillations, the magnetic field components X(t) and P(t) now exhibit strongly irregular, aperiodic behavior with prominent bursts of large amplitude, a clear indication of deterministic chaos. This behavior is more pronounced than in the R= 2 chaotic case. The phase space portraits in Figure 12 confirm the strongly chaotic nature. The (x, p) phase portrait is noticeably broader and more tangled compared to the case: R= 5 and ICs (0.8,0.8,0.01,0.01). The (X, P) portrait shows a much denser, more space-filling trajectory, which is a strong signature of a high-dimensional strange attractor. The three-dimensional projections clearly demonstrate the highly complex and tangled trajectory structure, indicating a loss of any simple or residual periodic structure. The Poincar´e sections in Figure 13 unambiguously confirm the deterministic chaos. Both sections display a dense, space-filling cloud of points, filling extended regions of the phase space more uniformly than in the previous case. The (x, p) section shows a thickened, annular structure with highly scattered points, and the (X, P) section displays a complex, scattered pattern lacking any distinguishable organized structure. The space-filling character of these sections is the definitive signature of a fully developed chaotic behavior 5 Neural network architecture and training Artificial Neural Networks (ANNs) are computational models inspired by biological neural systems, designed to learn complex patterns and relationships from data. Through an iterative training process, neural networks adjust their internal parameters to minimize prediction errors, making them particularly effective for modeling nonlinear dynamical systems where analytical solutions are difficult or impossible to obtain. In this study, the neural network was trained on time series data obtained by numerically integrating the differential equations of the dynamo system (1) using the ODE45 solver in 17 Figure 11: Time series of the dynamo system variables for the strongly chaotic regime with initial conditions (x0, p0, X0, P0) = (0.9,0.9,0.01,0.01) at R= 5. Figure 12: Phase space portraits for the strongly chaotic regime at R= 5: (top left) (x, p) projection showing a highly broadened, multi-layered structure; (top right) (X, P) projection revealing a dense strange attractor; (bottom) three-dimensional projections (x, p, X) and (X, P, x). 18 Figure 13: Poincare sections for the strongly chaotic regime with initial conditions (0.9,0.9,0.01,0.01) at R= 5. Left: (x, p) section showing a dense annular structure. Right: (X, P) section displaying a highly scattered, space-filling point distribution. MATLAB. From the numerical solution covering the time interval t∈[0,3000] with initial conditions x0= 0.9, p0= 0.9, X0= 0.01, P0= 0.01, a dataset of 17272 consecutive state vectors was generated. The training process involves learning the mapping from the current state of the system (x(t), p(t), X(t), P(t)) to its next state (x(t+∆t), p(t+∆t), X(t+∆t), P(t+ ∆t)), where ∆tis the time step determined by the adaptive solver. By repeatedly adjusting its parameters through backpropagation algorithms, the network learns to approximate the underlying dynamics of the chaotic system without requiring explicit knowledge of the governing differential equations. A classical Feed Forward Neural Network (FFNN) architecture was chosen and implemented for modeling the dynamo system (1). The key objective was to find an optimal network structure that balances approximation accuracy and computational complexity. 5.1 Choice of architecture and its justification The FFNN architecture was selected for this research based on several key advantages that make it well-suited for modeling dynamic systems like the dynamo. It is known that the universal approximation theorem [7, 8] proves that an FFNN with just one hidden layer containing a finite number of neurons and a nonlinear activation function can approximate any continuous function on a compact subset of Rnwith any desired accuracy. This fundamental property ensures that the network has sufficient theoretical capacity to model the complex nonlinear dynamics of the dynamo system. In addition, FFNNs have been successfully applied for predicting dynamical systems where the next state of the system is predicted based on the current state, demonstrating their effectiveness for time series modeling in chaotic systems. The feedforward architecture is particularly suitable for our one-step-ahead prediction task, where temporal dependencies are captured through the state-to-state mapping rather than requiring explicit memory mechanisms. The FFNN can be implemented without specialized toolboxes, using only basic MATLAB functions, which provides greater control over the training process and facilitates understanding of the underlying mechanisms. Finally, FFNNs train faster compared to more complex recurrent architectures like RNNs or LSTMs, making them computationally efficient for the large dataset 19 Figure 14: Schematic of the neural network architecture [16,4]. generated from the dynamo system simulations while still maintaining the capability to capture the essential dynamics of the chaotic attractor. 5.2 Detailed architecture description The selected FFNN architecture, denoted as [16,4], consists of three layers with a specific organization of neurons and connections between them. The input layer contains 4 neurons, each corresponding to one of the current state variables of the dynamo system: x(t), p(t), X(t), and P(t). These input neurons receive the instantaneous values of the system state and pass them forward through weighted connections to the next layer. The hidden layer comprises 16 neurons, each equipped with a hyperbolic tangent (tansig) activation function. This layer performs the crucial nonlinear transformation of the input data, extracting complex patterns and relationships between the state variables. The number of neurons in this layer (16) was chosen through empirical testing to provide sufficient representational capacity without introducing unnecessary computational complexity or risking overfitting. The output layer consists of 4 neurons with linear activation functions (purelin), producing the predicted values of the system state at the next time step: x(t+ ∆t), p(t+ ∆t), X(t+ ∆t), and P(t+ ∆t). The linear activation in the output layer allows the network to generate predictions across the full range of real values that the dynamo variables can assume, without artificial constraints imposed by bounded activation functions. The complete architecture is illustrated schematically in Figure 14. 20 5.3 Activation functions Two different activation functions were used in different layers of the network: •Hyperbolic tangent (tansig) in the hidden layer: f(x) = tanh(x) = ex−e−x ex+e−x This function outputs values in the range [−1,1]. It introduces the necessary nonlinearity that allows the network to learn complex patterns. Think of it as a ”squashing” function that converts any input to a value between -1 and 1. •Linear function (purelin) in the output layer: f(x) = x This simply passes the input through unchanged. It’s used in the output layer because our system variables (x, p, X, P) can take any real value, not just values between -1 and 1. 5.4 Number of trainable parameters The total number of adjustable parameters in the [16,4] architecture can be calculated by considering the weights and biases at each layer connection. Between the input layer (4 neurons) and the hidden layer (16 neurons), there are 4 ×16 = 64 weight parameters, with each input neuron connected to every hidden neuron. Additionally, the hidden layer requires 16 bias parameters, one for each hidden neuron. The connections from the hidden layer to the output layer contribute another 16 ×4 = 64 weight parameters, as each of the 16 hidden neurons connects to all 4 output neurons. Finally, the output layer adds 4 bias parameters, one per output neuron. Summing all these components yields a total of 64+16+64+4 = 148 trainable parameters. Each of these 148 parameters is iteratively adjusted during the training process through the backpropagation algorithm to minimize the prediction error between the network outputs and the target values from the dynamo system simulations. 5.5 Training algorithm: Levenberg-Marquardt backpropagation The network was trained using the Levenberg-Marquardt backpropagation algorithm [28, 29], which is particularly efficient for small to medium-sized networks and is widely recognized as one of the fastest training methods for this class of problems [30]. This algorithm represents a sophisticated optimization approach that combines the advantages of two classical methods: the gradient descent method, which takes small, cautious steps down the error surface, and the Gauss-Newton method, which attempts larger, more direct steps toward the minimum by utilizing second-order derivative information [31]. The training process follows an iterative procedure consisting of four main stages. First, during forward propagation, the network processes the input data by passing it through the layers, applying the activation functions, and generating predictions for the system’s next state. 21 Second, the error calculation stage computes the difference between these predictions and the actual target values obtained from the numerical solution of the dynamo equations, quantifying the network’s current performance. Third, in the backward propagation phase [32], the algorithm calculates the gradient of the error with respect to each parameter, determining how much each individual weight and bias contributed to the overall prediction error. Finally, the parameter update stage uses the Levenberg-Marquardt formula to adjust all weights and biases in a direction that reduces the error, with the algorithm adaptively switching between gradient descent behavior (when far from the minimum) and Gauss-Newton behavior (when close to the minimum) for optimal convergence. This iterative cycle repeats for many epochs, typically hundreds or thousands of iterations, with the algorithm continuously refining the parameter values until the network achieves the desired level of accuracy or until other stopping criteria are met, such as reaching a maximum number of epochs or observing no improvement in validation performance for a specified number of consecutive iterations. 5.6 Data preparation The dataset generated from the numerical integration of the dynamo equations was partitioned into three distinct subsets to ensure proper training, validation, and unbiased evaluation of the neural network model. The training set, comprising 70% of the total data (12090 samples), was used to iteratively adjust the network’s weights and biases through the backpropagation algorithm, allowing the model to learn the underlying dynamics of the system. The validation set, containing 15% of the data (2590 samples), served a dual purpose: monitoring the network’s performance during training to prevent overfitting and providing feedback for early stopping when no improvement was observed over consecutive epochs. The remaining 15% of the data (2592 samples) constituted the test set, which was kept completely separate from the training process and used only for final evaluation of the model’s predictive accuracy on previously unseen data. Prior to training, all input and output data were normalized using a standardization procedure to improve training efficiency and numerical stability. The normalization transformation is given by: xnorm =x−µ σ where µrepresents the mean value and σdenotes the standard deviation of each variable computed over the training set. This scaling procedure centers the data around zero with unit variance, ensuring that all four state variables (x, p, X, P) contribute equally to the learning process regardless of their original magnitudes or ranges. The same normalization parameters derived from the training set were consistently applied to both the validation and test sets to maintain data integrity and prevent information leakage between the subsets. 5.7 Architecture testing and selection of the best model To select the optimal configuration, a series of experiments was conducted, testing four different architectures: Architecture 1 (FFNN [16 4]) with one hidden layer containing 16 neurons; 22 Figure 15: MATLAB Neural Network Training interface showing Architecture 4 [20, 4] during the training process. The interface displays the network structure with 4 input neurons, one hidden layer with 20 neurons, and 4 output neurons. Training was performed using the Levenberg-Marquardt algorithm (trainlm) with mean squared error (mse) as the performance metric, completing after 500 epochs with a final performance of 3.61 ×10−7. 23 Table 1: Comparison of the performance of different neural network architectures No. Architecture Test MSE Best Performance (Val) 1 FFNN [16 4] 3.4247 ×10−58.2003 ×10−7 2 FFNN [32 16 4] 2.4290 ×10−41.5240 ×10−7 3 FFNN [24 12 4] 6.1564 ×10−52.8009 ×10−7 4 FFNN [20 4] 3.3631 ×10−53.6099 ×10−7 Architecture 2 (FFNN [32 16 4]) with two hidden layers containing 32 and 16 neurons respectively; Architecture 3 (FFNN [24 12 4]) with two hidden layers containing 24 and 12 neurons; and Architecture 4 (FFNN [20 4]) with one hidden layer containing 20 neurons. Each architecture was trained using identical training parameters and data partitioning to ensure a fair comparison. Figure 15 illustrates the MATLAB Neural Network Training interface during the training process of Architecture 4, showing the network structure, algorithm parameters, and real-time training progress monitoring. The evaluation was based on two complementary criteria. The first criterion was the Mean Squared Error (MSE) on the test set, which measures the model’s accuracy on completely unseen data and represents the ultimate indicator of practical predictive capability. The second criterion was the best validation performance, denoted as Best Performance (Val), which represents the minimum MSE achieved on the validation set during the entire training process. This metric is particularly important because it reflects the model’s ability to generalize beyond the training data before any overfitting occurs. During training, the network’s performance is continuously monitored on the validation set after each epoch, and the Best Performance value captures the lowest error achieved at the optimal point of training typically before the model begins to memorize noise or specific patterns from the training data. The comparison results for all tested architectures are presented in Table 1. Although the [16 4] architecture showed the smallest error on the test set (3.4247×10−5), the FFNN [32 16 4] architecture was selected as the best. This decision is based on its outstanding performance on the validation set (Best Performance = 1.5240 ×10−7), indicating the model’s superior generalization ability and minimal level of overfitting. 5.8 Training process and final model accuracy The training process of the best architecture [32 16 4] is visualized in Figure 16. The plot shows the error dynamics on the training, validation, and test sets. The absence of divergence between the curves and their stable convergence to a low value confirm the correctness of the training procedure and the lack of overfitting. The overall Mean Squared Error (Overall Test MSE) of the final model on the independent test set was 2.4290 ×10−4. Detailed metrics for each predicted variable are presented in Table 2. A visual assessment of the model’s accuracy on test data is presented in Figure 16. The plots demonstrate an almost complete overlap between the actual system trajectories (Real) and the values predicted by the neural network (Predicted). The scatter plot “Real vs Predicted Correlation” confirms the high linear correlation between the target and output values of the model. The phase portrait constructed from the predicted values of variables xand pcorrectly reproduces the structure of the system’s attractor (Figure 16, graph “Predicted Phase Portrait (x-p)”). 24 Figure 16: Training progress and accuracy assessment of the neural network. Upper left graph: error dynamics on the training, validation, and test sets. Other graphs: comparison of real and predicted values of variables on the test set, scatter plot, and predicted phase portrait. Table 2: Detailed prediction quality metrics for each system variable (architecture [32 16 4]) Variable MSE RMSE NMSE x2.7272 ×10−41.6514 ×10−23.7982 ×10−4 p2.7479 ×10−41.6577 ×10−23.9698 ×10−4 X2.4385 ×10−41.5616 ×10−21.1971 ×10−3 P1.8024 ×10−41.3425 ×10−28.9335 ×10−4 25