Python-Based Neural Network Modeling of Nonlinear Dynamo Chaos with Adam Optimization

Python-Based Neural Network Modeling of Nonlinear Dynamo Chaos with Adam Optimization M. I. Kopp1 December 10, 2025 1Institute for Single Crystals, NAS Ukraine, Nauky Ave. 60, Kharkiv 61072, Ukraine Abstract This study presents a Python-based implementation of feedforward artificial neural networks (FFNN) for modeling chaotic behavior in a nonlinear magneto-vortex dynamo system. The investigated system, derived from asymptotic multiscale nonlinear theory applied to temperature-stratified conducting media, consists of four coupled nonlinear differential equations describing the evolution of large-scale vortex and magnetic field structures. Numerical integration was performed using the ”odeint” solver from SciPy. The neural network was implemented using standard Python libraries (NumPy and SciPy) without relying on specialized deep learning frameworks such as TensorFlow or PyTorch. The network was trained using the Adam (Adaptive Moment Estimation) optimizer with the backpropagation algorithm. The optimal FFNN architecture [20, 4] with a single hidden layer demonstrated excellent prediction accuracy for both one-step-ahead forecasting and long-term autonomous predictions. This approach provides an efficient computational framework for analyzing complex nonlinear dynamo systems and offers insights into the predictability of chaotic regimes in nonlinear dynamo systems using modern optimization techniques. Keywords: magneto-vortex dynamo; chaotic behavior; feedforward neural network; Adam optimization; Python implementation 1 Introduction The application of artificial neural networks to chaotic dynamical systems has emerged as a powerful computational paradigm in recent years, offering data-driven approaches to problems where traditional analytical methods prove intractable. Chaotic systems, characterized by extreme sensitivity to initial conditions and complex long-term behavior, present unique challenges for both theoretical analysis and numerical prediction. Neural networks, with their universal approximation capabilities [1, 2], have demonstrated remarkable success in capturing the intricate dynamics of such systems across diverse fields ranging from meteorology to plasma physics. 1 Early applications of neural networks to chaos prediction focused primarily on classical lowdimensional attractors. The Lorenz system, a paradigmatic example of deterministic chaos arising from simplified atmospheric convection, has been extensively studied using various neural network architectures. Recent work by Chattopadhyay et al. [3] compared reservoir computing, standard artificial neural networks, and Long Short-Term Memory (LSTM) networks for predicting the multiscale Lorenz 96 system, demonstrating that different architectures exhibit varying strengths depending on the prediction horizon and system complexity. The MackeyGlass time-delay system, representing a different class of chaotic dynamics, was successfully modeled by Lamamra et al. [4] using multilayer perceptrons with topology optimization via the NSGA-II genetic algorithm, achieving mean squared errors on the order of 10−13. Electronic chaotic circuits have proven to be particularly amenable to neural network modeling due to their experimental accessibility and well-characterized dynamics. Al-Musawi et al. [5] implemented prediction of the chaotic Chua system using feedforward neural networks on FPGA hardware, achieving MSE values of 10−6and demonstrating practical applicability for chaos-based cryptographic systems. The Chua circuit [6], being one of the simplest electronic systems exhibiting chaos, continues to serve as a benchmark for testing neural network prediction capabilities [7]. Hardware implementations have extended beyond Chua’s circuit to more complex systems: Alcin et al. [8] developed VHDL implementations of ANN-based PehlivanUyaroglu chaotic generators, while Koyuncu et al. [9] created neuron libraries for rapid FPGA realization of chaotic systems including the R¨ossler attractor. Hydrodynamic and convective systems represent another important domain for neural network applications. Keles et al. [10] investigated the Rucklidge system, which models double convection in rotating fluids under magnetic fields [11], comparing feedforward neural networks (FFNN), layered recurrent neural networks (LRNN), and cascade-forward neural networks (CFNN). Their study revealed that Vanilla RNN with a sliding window approach achieved R2scores of 0.999999 and RMSE of 1.825 ×10−4, demonstrating the importance of architectural choice for specific dynamical systems. Cimen et al. [12] presented an innovative approach combining computer vision with neural networks for modeling chaotic motion in the modified Lorenz system, extracting dynamics directly from video sequences. Mechanical systems with chaotic behavior have also attracted attention. Ramachandruni et al. [13] conducted a comprehensive comparative study of machine learning approaches for multipendulum chaos, evaluating ten different models including various neural network architectures. LSTM networks demonstrated superior performance with RMSE on the order of 10−2, and importantly, the study showed that properly trained networks could predict chaotic motion for completely unseen initial conditions a crucial capability for practical applications. This work highlighted the advantage of recurrent architectures for systems with strong temporal dependencies. The architecture and training methodology significantly impact neural network performance on chaotic systems. Zhang and Lei [14] examined training efficiency across different multilayer configurations for chaotic systems, providing insights into optimal network depth and width. The choice of optimization algorithm plays a critical role: traditional approaches often employ the Levenberg-Marquardt algorithm [15, 16], known for its efficiency in small to medium networks [17], while modern deep learning predominantly utilizes adaptive gradient methods. The Adam (Adaptive Moment Estimation) optimizer [18], combining benefits of AdaGrad and RMSProp, has become the de facto standard in contemporary neural network training due 2 to its computational efficiency, low memory requirements, and robustness to hyperparameter choices. Recurrent neural networks, particularly LSTMs [19], have shown exceptional capability in handling temporal sequences due to their explicit memory mechanisms. However, feedforward architectures remain valuable for specific prediction tasks where the mapping from current state to next state sufficiently captures system dynamics. The universal approximation theorem guarantees that feedforward networks with sufficient hidden neurons can approximate any continuous function to arbitrary accuracy, making them theoretically suitable for modeling deterministic chaos where future states are functionally determined by present states. Magnetohydrodynamic (MHD) systems, which couple fluid dynamics with magnetic field evolution, represent particularly challenging targets for neural network modeling due to their multiphysics nature and high dimensionality. Classical dynamo systems that have been studied include the Rikitake two-disk dynamo [20, 21] for modeling geomagnetic field reversals and Bullard’s homopolar dynamo [22, 23]. However, these simplified models often rely on phenomenological assumptions rather than first-principles derivations from magnetohydrodynamics equations. The nonlinear dynamo system investigated in this work [24] was derived rigorously from the full averaged MHD equations using asymptotic multiscale analysis, providing a more physically grounded foundation. Despite extensive applications of neural networks to various chaotic systems, there remains a notable gap in the literature regarding the modeling of complex MHD dynamo systems using modern optimization techniques. Our recent work [25] demonstrated the feasibility of neural network prediction for a four-dimensional magneto-vortex dynamo system using MATLAB implementation with the Levenberg-Marquardt training algorithm, achieving MSE values on the order of 10−4to 10−7depending on architecture complexity. Furthermore, the present study extends the work of [26] by implementing the same dynamo system model in Python using the Adam optimization algorithm. The aims of this work are as follows: (1) to develop a Python-based feedforward neural network using only standard scientific libraries for modeling the chaotic magneto-vortex dynamo system; (2) to evaluate the effectiveness of Adam optimization for training neural networks on chaotic MHD dynamics; and (3) to compare the performance and implementation characteristics with our previous Levenberg-Marquardt approach, providing insights into optimization algorithm selection for chaotic system prediction tasks. 2 Mathematical Model of Nonlinear Dynamo The chaotic dynamo system investigated in this study originates from the rigorous nonlinear dynamo theory developed by Kopp et al. [24]. Unlike phenomenological models such as the Rikitake two-disk dynamo or Bullard’s homopolar dynamo, this system is derived directly from the full averaged equations of magnetohydrodynamics through systematic multiscale asymptotic analysis. The derivation applies to temperature-stratified electrically conducting media subjected to small-scale helical external forcing, which acts as a source of helical turbulence with low Reynolds numbers. Following the application of quasi-two-dimensional approximation and averaging over fast spatial scales, the evolution of large-scale magnetic and velocity fields 3 reduces to a system 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) where x(t) and p(t) represent the amplitudes of large-scale velocity field components, while X(t) and P(t) represent the corresponding magnetic field components. The constants c1= 0.01, c2= 0.01, c3= 0.001, and c4= 0.001 arise from the integration of non-stationary terms in the original equations [24].The coefficients α1,α2,β1, and β2are not constant parameters but rather complex nonlinear functions of the system state and physical parameters. These coefficients encapsulate the combined effects of hydrodynamic and magnetohydrodynamic α-effects in the presence of thermal stratification. The hydrodynamic α-effect coefficients governing velocity field evolution 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, 4 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. •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. From the perspective of dynamo theory, the chaotic regimes exhibited by this system possess important physical significance. The irregular, aperiodic solutions naturally model intermittent magnetic field generation processes, mirroring the sporadic bursts of activity punctuated by quiescent periods observed in astrophysical and geophysical dynamos. The Rayleigh number R functions as a bifurcation parameter, controlling the system’s evolution from ordered regular states characterized by coherent, stable magnetic field patterns to disordered turbulent-like states featuring unpredictable field reversals. Energy transfer between the kinetic degrees of freedom (x, p) and magnetic degrees of freedom (X, P) occurs through the nonlinear coupling encoded in the α-effect coefficients, while simultaneously facilitating cross-scale interactions. As this energy exchange intensifies and becomes irregularly modulated, the system undergoes a transition into chaotic regimes. This dynamical behavior offers a mechanistic framework for understanding phenomena such as irregular polarity reversals in planetary and stellar magnetic fields, the spontaneous emergence of coherent magnetic structures from turbulent backgrounds, and abrupt transitions between qualitatively different dynamo operating modes all fundamental features observed in natural magnetohydrodynamic systems. 3 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 odeint solver from SciPy in Python. From the numerical solution covering the time interval t∈[0,3000] with 5 initial conditions x0= 0.9, p0= 0.9, X0= 0.01, P0= 0.01 at R= 5, a dataset of 30,000 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 numerical solver. By repeatedly adjusting its parameters through backpropagation algorithms with Adam optimizer, the network learns to approximate the underlying dynamics of the chaotic system without requiring explicit knowledge of the governing differential equations. A Feed Forward Neural Network (FFNN) architecture was chosen and implemented for modeling the dynamo system. The key objective was to find an optimal network structure that balances approximation accuracy and computational complexity. 3.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. The universal approximation theorem [1, 2] 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. FFNNs have been successfully applied for predicting dynamical systems where the next state 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 stateto-state mapping rather than requiring explicit memory mechanisms like RNNs or LSTMs. The implementation utilized pure NumPy and SciPy libraries, avoiding dependencies on TensorFlow or other deep learning frameworks. This approach provides greater control over the training process, facilitates understanding of the underlying mechanisms, and ensures compatibility with standard scientific computing environments. The custom implementation includes an Adam optimizer for efficient gradient-based optimization and sophisticated early stopping mechanisms to prevent overfitting. 3.2 Detailed architecture description The selected FFNN architecture, denoted as [20, 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 20 neurons, each equipped with a hyperbolic tangent (tanh) 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 (20) 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, producing the predicted values of the 6 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. 3.3 Activation functions Two different activation functions were used in different layers of the network: •Hyperbolic tangent (tanh) 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. The tanh function is particularly wellsuited for neural networks due to its zero-centered output and smooth gradient properties. •Linear function in the output layer: f(x) = x This simply passes the input through unchanged. It is 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. 3.4 Number of trainable parameters The total number of adjustable parameters in the [20, 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 (20 neurons), there are 4 ×20 = 80 weight parameters, with each input neuron connected to every hidden neuron. Additionally, the hidden layer requires 20 bias parameters, one for each hidden neuron. The connections from the hidden layer to the output layer contribute another 20 ×4 = 80 weight parameters, as each of the 20 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 80+20+80+4 = 184 trainable parameters. Each of these 184 parameters is iteratively adjusted during the training process through the backpropagation algorithm with Adam optimization to minimize the prediction error between the network outputs and the target values from the dynamo system simulations. 3.5 Training algorithm: Adam optimization with backpropagation The network was trained using the Adam (Adaptive Moment Estimation) optimizer [18], which is particularly effective for training neural networks on large datasets and has become the de facto standard in modern deep learning applications. This algorithm represents an advanced optimization approach that combines the advantages of two popular methods: AdaGrad, which 7 adapts learning rates based on historical gradients, and RMSProp, which uses a moving average of squared gradients to normalize updates. The Adam optimizer maintains two moving averages for each parameter: the first moment (mean) of the gradients and the second moment (uncentered variance) of the gradients. These estimates are updated using exponential decay rates β1= 0.9 and β2= 0.999, respectively. The algorithm includes bias correction mechanisms to account for initialization bias in the early stages of training. The base learning rate was set to α= 0.001, with adaptive per-parameter learning rates computed from the moment estimates. The training process follows an iterative procedure with mini-batch gradient descent. 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. Second, the error calculation stage computes the mean squared error between these predictions and the actual target values. Third, in the backward propagation phase [27], the algorithm calculates the gradient of the error with respect to each parameter. Finally, the parameter update stage uses the Adam formulas to adjust all weights and biases adaptively, with the algorithm automatically adjusting step sizes based on the gradient history for optimal convergence. This iterative cycle repeats for up to 1000 epochs with mini-batch size of 256 samples, with the algorithm continuously refining the parameter values until convergence or until early stopping criteria are met (no improvement in validation performance for 50 consecutive epochs). 3.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 (20,999 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 (4,499 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 (4,501 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 z-score 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. 8 3.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; 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, data partitioning, and optimization settings to ensure a fair comparison. 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 validation MSE achieved during training, which reflects the model’s ability to generalize beyond the training data before any overfitting occurs. The comparison results for all tested architectures are presented in Table 1. Table 1: Comparison of the performance of different neural network architectures No. Architecture Test MSE Validation MSE 1 FFNN [16, 4] 3.8547 ×10−63.9124 ×10−6 2 FFNN [32, 16, 4] 4.2183 ×10−64.1856 ×10−6 3 FFNN [24, 12, 4] 3.9821 ×10−64.0234 ×10−6 4 FFNN [20, 4] 3.7219 ×10−63.8011 ×10−6 The FFNN [20, 4] architecture was selected as the best model based on its superior performance on both the test set (Test MSE = 3.7219 ×10−6) and validation set (Validation MSE = 3.8011×10−6). This architecture demonstrates the best balance between model complexity and generalization ability, achieving the lowest prediction errors while maintaining computational efficiency. The single hidden layer with 20 neurons proved sufficient to capture the essential nonlinear dynamics of the dynamo system without overfitting. 3.8 Training process and final model accuracy The training process of the best architecture [20, 4] converged successfully with early stopping triggered after observing no improvement in validation performance. Visual assessment of the model’s predictions is presented in Figure 1, which shows excellent agreement between the actual system trajectories and the neural network predictions across all four state variables. The overall Test MSE of the final model on the independent test set was 3.7219 ×10−6. Detailed metrics for each predicted variable are presented in Table 2. The plots in Figure 1 demonstrate an almost complete overlap between the actual system trajectories (blue solid lines) and the values predicted by the neural network (red dashed lines) for all four variables. The scatter plot “Correlation” confirms the high linear correlation between the target and output values of the model, with all points closely following the ideal diagonal line. The phase portrait constructed from the predicted values of variables xand pcorrectly reproduces the structure of the system’s attractor, showing that the network has learned the essential topological features of the chaotic dynamics. The error distribution histogram shows a 9