DynPy—Python Library for Mechanical and Electrical Engineering: An Assessment with Coupled Electro-Mechanical Direct Current Motor Model

Abstract

DynPy is an open-source library implemented in Python (version 3.10.12) programming language which aims to provide a versatile set of functionalities for mechanical and electrical engineers. It enables the user to model, solve, simulate, and report analysis of dynamic systems with the use of a single environment. The DynPy library comes with a predefined collection of ready-to-use mechanical and electrical systems. A proprietary approach to creating new systems by combining independent elements defined as classes, such as masses, springs, dampers, resistors, capacitors, inductors, and more, allows for the quick creation of new, or the modification of existing systems. In the paper examples for obtaining analytical and numerical solutions of the systems described with ordinary differential equations were presented. The assessment of solver accuracy was conducted utilising a coupled electro-mechanical model of a direct current motor, with MATLAB/Simulink (R2022b) used as a reference tool. The model was solved in DynPy with the hybrid analytical–numerical method and fully analytically, while in MATLAB/Simulink strictly numerical simulations were run. The comparison of the results obtained from both tools not only proved the credibility of the developed library but also showed its superiority in specific conditions.

1. Introduction

Simulations are an essential part of modern science, as they provide a multitude of benefits. Starting with the possibility of analysing very complicated processes, impossible or barely impossible to be conducted in real-life conditions, to ‘what if’ scenarios, hypotheses evaluation, understanding of processes, and optimisations, ending with designing and testing new technologies and solutions. The goal of this article is to present the abilities of the independent, Python-programming-language-based DynPy (‘DynamicPython’) [1] library and compare the results obtained from simulations of a dynamic system conducted with well-known, field-tested commercial software. A growing demand for clean energy and transportation introduces new challenges to multiple engineering branches [2,3,4,5]. As a result, a continuous seeking for more reliable, durable, safe, and effective hardware design is seen. This usually means a necessity for analysing yet more complicated phenomena [6], for which nontrivial approaches are needed, most likely to be investigated with the use of some software.

DynPy is an open-source project aimed at providing an engineering tool that allows, among other functionalities: performing symbolic computations, conducting numerical simulations, and creating comprehensive reports of conducted calculations. It supports a few useful and suitable typesetting systems for reporting engineering and scientific activities, e.g., Markdown and LaTeX. Another advantage of the developed library is the significant number of dynamic components it provides. This allows for quickly simulating the behaviour of systems described by mathematical dependencies characterizing their dynamic states. DynPy is conceptualised, developed, and maintained by the authors and other developers. Throughout the years it has been utilized by various researchers in the field of mechanical engineering. Due to the constantly growing capabilities of the developed tool, its use in research could be categorized regarding the formerly available set of functionalities. In [7], the proprietary methodology for creating dynamic systems, as well as numerical and analytical solving algorithms, were utilised. Further, in [8,9], chosen elements of a text editing submodule were additionally used for adding formulas and other minor typesetting activities. Finally, the article in [10] was the first one to be created as a whole strictly with DynPy’s functionalities—from model creation, through simulations, to manuscript generation. The authors are not aware of any other examples but do not exclude such a possibility, as the library is available in a public GitHub (version 2.34.1) repository. Such circumstances are the effect of the current article being the first one explicitly about DynPy, as the state of development allows for presentation to the broader public.

In the case of electrical systems, the description of their dynamics most often takes the form of ordinary differential equations (further referred to as ODE). Due to the complexity of the formulas, solutions can be found in analytical forms or with the use of numerical algorithms. The biggest advantage of the analytical approach lies in the provision of exact solutions for the whole domain, which ensures high repeatability and, thus, reliability of achieved results. Even though it is not always possible to directly obtain closed-form formulas, the stated aspects push scientists to find different mathematical tools, which lead to overwork of hitherto numerically analysed phenomena with analytical methods also in the case of electrical circuits [11,12,13,14]. Although, in the current state, the DynPy library does not provide as specialised functionalities as presented in the cited literature, it is crucial to note its ability in this aspect. Firstly, a broad spectrum of typical electrical elements is predefined with the possibility of further enlargement. Secondly, basic analytical and numerical solvers are implemented. Thanks to the utilisation of an object-oriented paradigm and open-source code, more sophisticated approaches can be introduced on the user’s demand.

As noted earlier, the ability to achieve analytical solutions is case-specific and highly depends on the type of equations describing the system. Numerical simulations became such a powerful tool due to the provision of approximated results to problems, which could not be worked out by traditional means [15,16,17]. Multiple types of numerical solving methods have been introduced throughout the years [18]. In the case of ODE, three main families of algorithms can be recognised: explicit single-step, explicit multistep, and their implicit analogues dedicated to stiff systems of differential equations [19]. Typical single-step methods for ODE may be divided into the Runge–Kutta (RK, fixed step size) and Runge–Kutta–Fehlberg (RKF, variable step size) approaches. Examples of their utilisation in analyses of electrical circuits in both nonstiff and stiff ODEs can be found in [20,21,22,23,24,25,26,27,28]. Moving on to multistep algorithms, again two main groups could be highlighted, known as Adams and predictor–corrector methods. The Adams family methods in analyses of electrical circuits were presented in [29,30,31,32,33,34,35]. Predictor–corrector (PC) algorithms include both explicit and implicit algorithms at once, which results in obtaining a set of properties connecting their benefits. The sequence of operations is as follows: prediction (explicit, P), evaluation (E), correction (implicit, C), evaluation (E) [36,37]. The described rule is outright applicable to predictor–corrector algorithms based on Adams–Bashforth and Adams–Moulton methods, but constant research on the effects of numerical approximations in the development of different approaches including explicit–explicit or implicit–implicit PC algorithms [38,39,40,41,42,43,44,45]. The DynPy environment in its current state provides six types of numerical solving techniques, inherited from the SciPy (version 1.15.0) [46] module: RK45—explicit Runge–Kutta method of order $5\left(4\right)$; RK23—explicit Runge–Kutta method of order $3\left(2\right)$; DOP853—explicit Runge–Kutta method of order 8; Radau—implicit Runge-Kutta method of the Radau IIA family of order 5; BDF—implicit multistep variable-order (1 to 5) method based on a backward differentiation formula for the derivative approximation; and the LSODA —Adams/BDF method with automatic stiffness detection and switching. Simultaneously with other development and maintenance works, actions are taken to create DynPy-embedded, standalone versions of these and other useful algorithms, so a package can become more versatile and less susceptible to third-party changes in the code, possibly due to the previously stated inheritance.

Solving and simulating electrical circuit behaviour can be performed with the use of multiple software tools available on the market. While some of them strictly require the user to formulate code instructions, others have more advanced Graphical User Interfaces which provide the possibility to create a representation of analysed problems by block diagrams or 3D models. Vendors such as MathWorks (Natick, MA, USA) or Wolfram (Champaign, IL, USA) offer platforms of both types [47,48,49,50]. Each of them uses individual language, mutual between different developed environments. The mentioned programs serve a multipurpose function, being used by engineers of imaginably all branches. Throughout the years, many additional packages and add-ins were introduced, so the functionalities of interest for specific sectors of the industry could be grouped and state better fitment to the consumers [51,52]. When such universal tools are not found as a suitable solution for any reason, more task-oriented software is accessible, such as TINA, whose scope is restricted to electrical/electronic circuit modelling and simulations [53]. All of the highlighted examples, despite their undeniable advantages, share one undisputable drawback—to obtain a full experience of their abilities, whether in general or in exact context, a license has to be acquired, which often does not include additional toolboxes or add-ins offered independently at different prices. Although it is possible to find free and open-source alternatives, such as Scilab [54] or CircuitLab [55], they are not flawless—less developed or completely lack predefined electrical elements, one-off types of interfaces, and a smaller group of users effects in worse community support (what becomes a serious problem when diving into analyses of more complex phenomena, which, may even outcome as too complicated for the utilised software), and so on.

2. Materials and Methods

2.1. DynPy Library—General Overview

The decision to use Python programming language for the implementation of the DynPy library was dictated by its numerous advantages. It is used by hundreds of thousands of programmers worldwide in areas such as web scripting, system programming, user interfaces, and product customisation [56,57], and the opinion is becoming widespread that it is an optimal choice for scientists and engineers seeking a high-level language for writing scientific applications [58]. As a high-level language, its syntax is more comprehensible for users, putting more focus on the task that has to be carried out rather than analysing the aspects of computing areas that naturally are not taken into consideration by nonprogrammers, such as memory management. The general purpose character of the language results in a wide variety of up-to-date libraries addressing multiple fields, such as scientific computing, data analysis, data visualisation, and more. Thanks to that, it was possible to use pre-existing modules as a foundation to implement functionalities desired functionalities, instead of writing the code from scratch. NumPy (version 2.2.1), pandas (version 2.2.3), SciPy (version 1.15.0), and SymPy (version 1.13.3) [59] are concurrently used for data processing, computations, simulations, and various mathematical operations. Pandas, in collaboration with Matplotlib (version 3.10.0), is used for presenting data and results in tabular and graphical form, while PyLaTeX (version 1.4.2) provides crucial methods required for text formatting and generating .tex and .pdf files. Therefore, Jupyter Notebook (version 7.3) is perfectly suited as a working environment with the DynPy library [60]. The ability to execute code in separate cells allows for conducting calculations in a specified order. The visualisation of operations facilitates understanding and control over the performed computations. Additionally, the option of combining code, text, and images makes the developed document readable and more accessible to other users.

In the process of selecting a programming paradigm from multiple which Python (version 3.10.12) offers, the decision was made to adopt an object-oriented programming approach, as it enables the creation of structured and modular solutions. Object-oriented code is more transparent and flexible, facilitating project management and development and providing versatility and usefulness. In the context of the project, it enables easier management of multiple modules and provides a structure that facilitates collaboration among multiple software developers.

It is worth mentioning at this point that the work on the DynPy library project follows an Agile approach, which is based on flexibility, adaptability, and continuous adjustment to changing needs. In contrast to the classical technique, which assumes planning and implementation in one cycle, the Agile methodology focuses on short-cycle stages. The goal is divided into smaller parts, and work is carried out step by step with the possibility of adapting to new conditions. This not only provides greater flexibility but also enables faster delivery of useful engineering tools [61]. The project as a whole is not presented further in this work. Instead, individual sections are presented and discussed.

The algorithms implemented in the library allow the solving and simulating of linear and nonlinear systems, both numerically and analytically. The numerical techniques inherited from SciPy provide a possibility for the analysis of linear and weakly and strongly nonlinear problems. The analytical approaches are defined in dynpy.solvers.linear and dynpy.solvers.nonlinear submodules. On top of the typical approaches utilised when working on linear ODEs, a few useful methods are defined which enable dealing with nonlinear problems. Firstly, a functionality of linearisation around an introduced point is provided, where calling the linearized method on the provided equations of motion returns a new set of formulas based on Taylor polynomial expansion of the nonlinear terms. Secondly, perturbational solvers are implemented, such as the Multiple Scales Method. The effectiveness in dealing with such problems in both numerical and analytical ways was tested in the previously referred literature, namely [7,8]. In detail, the library allows for the application of a few approaches for solving ODEs according to the class of the problem under analysis:

• Linear: DynPy provides analytical and numerical solvers for linear systems. The scope of the application of analytical solvers is limited by the order of the equation and its parameters. The system of two second-order equations can be solved analytically in general cases (all parameters given as arbitrary symbols). For selected cases of parameter arrangement, equations of higher order are also possible to solve. All other linear cases can be efficiently solved with numerical solvers as described previously.
• Nonlinear: This class of equations can be mostly solved by numerical solvers provided in DynPy. It results in the TimeDataFrame object (DynPy modification of pandas DataFrame object, providing data structure optimised for handling results in time and frequency domain), a uniform interface for further analysis (e.g., calculations, transformations) and manipulation (embedding the tables and graphs in the report) of results. This approach is consistent with an analytical approach, which the user can handle in a standardised way as presented before. Moreover, the selected cases are possible to solve in an analytical way if the equation is solvable from the theoretical point of view. First-order and Bernoulli equations are implemented currently.
• Perturbational: Selected cases can be solved with the usage of the Multiple Scale Method, which is an analytical method for solving systems with weak nonlinearities or parametrical properties. The equation has to be described by the perturbational parameter (from a practical point of view—a small physical part of the system). DynPy provides several tools that allow to create a perturbational form of the equation and solve it with the same methodology and in the same environment as previously described cases.

The DynPy library is equipped with a series of dynamic models that simulate real systems, the base of which is a ComposedSystem class. An electrical circuit can be constructed by leveraging its functionalities and incorporating components provided by the dynpy.models.electric.elements module (e.g., Resistor) or by utilizing predefined systems, the availability of which continues to expand. An exemplary electrical circuit represents a Thevenin model of the lithium-ion battery, a code snippet of which is provided in Appendix A. Using such a class eliminates the need for independently designing a system while also allowing for easy analysis. Additionally, the models include graphical visualisations, as shown in Figure 1. Expanding the DynPy library with various models is an effective way to enhance its functionality. Thanks to a standardised approach to building models, adding more becomes simpler.

2.2. Solvers

The core of the simulational capabilities of the DynPy library lies in its dynpy.solvers module. Technically, the most important class of all defined in its submodules is ODESystem, which serves as a template for creating objects in the form of differential equations. In the scope of the class body, there is an example demonstrating its use to create a differential equation and subsequently solve it symbolically and numerically using the tools provided by the library, which can be displayed by calling the help(ODESystem) command in the interpreter.

Preparation of the defined model for simulations in the DynPy environment starts with the creation of ODESystem class instance. To satisfy its necessary arguments, equations of motion of the system have to be introduced, as well as dependent variables and the degree of the differential equation. The ODEs can be introduced in two different ways, either by directly defining them as an input argument or by creating a dynamic system class and using its _eoms property defined in the parent class ComposedSystem. The latter method is presented, as it shows DynPy’s functionalities in a wider aspect. A code snippet of such a constructed DCMotor class is presented in Appendix B. Such a class can be defined directly in the interpreter or (preferably) in the proper module. In the case of providing the definition in a module, it is necessary to import it into an interpreter. Next, creating a noncalled instance is suggested, so it gives access to inherited and internally defined objects. A properly constructed class enables the return of equations of motion, which can be directly introduced as an argument to ODESystem. With a similar approach, a matrix of dependent variables can be formulated. The last step is the provision of an ODE order. As a result, a complete representation of the ODESystem type object is obtained. A graphical representation of the described steps and generated output is shown in Figure 2.

The prepared object allows for solving the system of equations analytically using the solution method (Figure 3). By default, a matrix equation is returned, where the left-hand side represents dependent variables and the right-hand side represents a complete solution. An additional as_eq_list method was developed, which converts the original output into a list of equations of motion as Sympy’s Eq class type, so they can be easily reported in the document. For clarity of view, a display command, being a part of the IPython module, was utilised to produce output, including an asterisk to unpack the list elements. The symbolic solution is very extensive, but substituting numerical values provides simplification (Figure 4).

To generate a dataset of numerical values for a chosen time span, three additional methods are called on the solution. Firstly, with_ics, which takes a list of initial conditions (ICs) as an argument. Based on the provided ICs, integration constants are evaluated and substituted. Secondly, by utilizing the subs method, the system’s symbolical parameters are replaced with numbers provided by a user in the form of a dictionary. Lastly, compute_solution is called on the predefined time span. Executed code results in a table of analytical solutions and their derivatives calculated for subsequent time stamps (Figure 5).

An approach to obtain a strictly numerical solution is presented in Figure 6. The substitution of parameter values is provided on equations of motion derived from the DCMotor class instance. Later, the resulting object is treated with a numerised method, which prepares the model for numerical computations. In the same line, with_ics is utilised to input initial conditions, and compute_solution with a defined time span solves the system, returning a dataset of results.

The presented methodology is well suited for analysing a system with a single set of chosen parameter values. The more sophisticated way of conducting simulations, which enables the possibility of analysing the changes in the system’s behaviour for multiple values of a single parameter is provided by NumericalAnalysisDataFrame class. It includes the from_model method, which takes as arguments an object of ODESystem type, analysed parameter (symbol) and a list of its values, reference parameters values (a dictionary, excluding the parameter under investigation), and a list of dependent variables and an index of DataFrame to be generated (where time_span is a set of time stamps, and dcmotor.ivar is a symbolic definition of time as an independent variable derived from DCMotor class, with ivar being independent variable). Due to the necessity of providing multiple sets of initial conditions to properly initiate a DataFrame, the additional method insert_ics was developed so that a simple list of values can be introduced. Such a setup directly outputs a result of numerical simulations (Figure 7) after executing perform_simulations command. To obtain numerical values calculated for the analytical solution, slight modifications have to be performed, as described in Figure 8. For all of the presented approaches, plots are easily generated by adding .plot() after compute_solution or perform_simulations, which can be additionally formatted to show chosen results with pandas library options (as a visible droplevel command, for example).

3. Theory and Calculation

3.1. DC Motor—Tested System

The accuracy and effectiveness of the developed DynPy library in the aspect of solving systems of ODEs and conducting simulations were investigated with a model of a direct current electric motor described by a system of two differential equations.

Based on the equivalent circuit (Figure 9) and Kirchhoff’s second law, the electrical equation of the motor was derived:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill U_z=E+U_{Lw}+U_{Rw}\end{array}\end{array}$$

(1)

where:

• $U_{Lw}$—voltage related to the rotor inductance;
• E—electro-motive force of induction;
• $U_z$—voltage supplying the rotor;
• $U_{Rw}$—voltage across the rotor winding resistance.

The voltage across the resistance of the rotor windings is proportional to the current flowing through them:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill U_{Rw}=R_w{i}_w\end{array}\end{array}$$

(2)

where the following are used:

• $i_w$—rotor winding current;
• $R_w$—equivalent resistance of the rotor windings.

The voltage across the rotor winding inductance is proportional to the rate of change in the current flowing through it (magnetic circuit losses are neglected):

$$\begin{array}{c}\hfill \begin{array}{c}\hfill U_{Lw}=L_w{\dot{i}}_w\end{array}\end{array}$$

(3)

where the following are used:

• $L_w$—equivalent inductance of the rotor windings;
• ${\dot{i}}_w$—rotor winding current rate of change.

When the rotor rotates, a counter-electro-motive force (EMF) is induced in its windings, the value of which is proportional to the angular velocity of the rotor:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill E=k_e{\omega }_s\end{array}\end{array}$$

(4)

where the following are used:

• ${\omega }_s$—angular velocity of the rotor;
• $k_e$—back electro-motive force constant.

By substituting each component into Equation (1), we obtain

$$\begin{array}{c}\hfill \begin{array}{c}\hfill U_z=L_w{\dot{i}}_w+R_w{i}_w+k_e{\omega }_s\end{array}\end{array}$$

(5)

The rotor torque has to overcome torques opposing the movement and is proportional to the current flowing through it. Assuming that the magnetic flux of the stator remains constant, two formulas were derived:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill M_s=M_a+M_l+M_r\end{array}\end{array}$$

(6)

$$\begin{array}{c}\hfill \begin{array}{c}\hfill M_s=k_m{i}_w\end{array}\end{array}$$

(7)

where the following are used:

• $M_s$—rotor torque;
• $M_l$—engine load torque;
• $k_m$—torque constant;
• $M_r$—rotor motion resistance torque;
• $M_a$—rotor angular acceleration torque.

The torque values associated with the angular acceleration and the rotor’s motion can be expressed as

$$\begin{array}{c}\hfill \begin{array}{c}\hfill M_a=J{\dot{\omega }}_s\end{array}\end{array}$$

(8)

$$\begin{array}{c}\hfill \begin{array}{c}\hfill M_r=B{\omega }_s\end{array}\end{array}$$

(9)

where the following are used:

• B—coefficient of viscous friction reduced to the rotor shaft;
• J—moment of inertia reduced to the rotor shaft;
• ${\dot{\omega }}_s$—angular acceleration of the rotor.

Substituting the individual components described in Formulas (8) and (9) into Equation (6) we obtain

$$\begin{array}{c}\hfill \begin{array}{c}\hfill k_m{i}_w=M_l+B{\omega }_s+J{\dot{\omega }}_s\end{array}\end{array}$$

(10)

As a result of the performed steps, equations of motion of the system were obtained:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill -U_z+L_w{\dot{i}}_w+R_w{i}_w+k_e{\omega }_s=0\end{array}\end{array}$$

(11)

$$\begin{array}{c}\hfill \begin{array}{c}\hfill M_l+B{\omega }_s+J{\dot{\omega }}_s-k_m{i}_w=0\end{array}\end{array}$$

(12)

The proposed model was solved with two techniques—with a hybrid approach, which utilises both analytical and numerical solutions at once (HAN), and strictly analytically, with forcing input as the Heaviside function. The HAN solution was completed in two steps. The system was first solved analytically for constant external forcing with the use of DynPy library (the solution was not calculated manually). Later, the obtained solution was updated numerically by alternately introducing dictionaries with parameter values filled in such a manner, so a step character of forcing could be imitated. As input functions in the model, the rotor supply voltage ($U_z$) and the motor load torque ($M_l$) were considered:

$$U_z\left(t\right)=\left\{\begin{array}{cc}10\hfill & \mathrm{for}t=\left\{\right(0,2)\cup (4,6)\cup (8,10\left)\right\}\hfill \\ 0\hfill & \mathrm{for}t=\left\{\right(2,4)\cup (6,8\left)\right\}\hfill \end{array}\right.$$

(13)

$$M_{\mathrm{l}}\left(t\right)=\left\{\begin{array}{cc}0.2\hfill & \mathrm{for}t=\left\{\right(1,3)\cup (5,7)\cup (9,10\left)\right\}\hfill \\ 0\hfill & \mathrm{for}t=\left\{\right(0,1)\cup (3,5)\cup (7,9\left)\right\}\hfill \end{array}\right.$$

(14)

Both are represented by step functions described in Equations (13) and (14), with graphical interpretation shown in the plots (Figure 10 and Figure 11). A detailed view of selected mathematical dependencies of the analysed system derived from the DCMotor class, and the full form of the analytical solution are defined in Appendix C.

The equations of motion for a system with Heaviside step function forcing have the following form:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill L_w{\dot{i}}_w+R_w{i}_w+k_e{\omega }_s-U_z\theta \left(sin\left({\displaystyle \frac{2\pi t}{T}}\right)\right)=0\end{array}\end{array}$$

(15)

$$\begin{array}{c}\hfill \begin{array}{c}\hfill B{\omega }_s+J{\dot{\omega }}_s+M_l\theta \left(-cos\left({\displaystyle \frac{2\pi t}{T}}\right)\right)-k_m{i}_w=0\end{array}\end{array}$$

(16)

where the following are used:

• $\theta$—Heaviside function;
• T—cycle length.

The step function takes a sinewave as an argument, to which additional conditions were provided for proper representation of a demanded character of input voltage and motor load. The analytical solution of Equations (15) and (16) is provided in Appendix C, while the code snippet with the implementation of the HAN method is in Appendix D.

3.2. DC Motor—Conducting Simulations Using the DynPy Library

The parameter values of the model were selected and presented in

Table 1. Parameter values adopted for calculations.

Parameter	Unit	Value
$R_w$	$\Omega$	2.0
$L_w$	H	0.1
$k_e$	${\displaystyle \frac{\mathrm{Vs}}{\mathrm{rad}}}$	0.1
$k_m$	${\displaystyle \frac{\mathrm{Nm}}{\mathrm{A}}}$	0.1
J	${\displaystyle \frac{{\mathrm{kgm}}^2}{{\mathrm{s}}^2}}$	0.1
B	${\displaystyle \frac{\mathrm{Nms}}{\mathrm{rad}}}$	0.5
T	s	4

. Results of simulations for armature current and rotor angular velocity obtained with the HAN method are provided in Figure 12 and Figure 13, while plots presented in Figure 14 and Figure 15 apply to the fully analytical solution.

In Figure 16 and Figure 17, corresponding plots were juxtaposed. The fully analytical solution is visibly smoother in areas where rapid changes in analysed values occur (selected time stamps of 1 $\mathrm{s}$ interval are presented in

Table 2. Results of the simulation conducted using the DynPy library—fully analytical method.

$\mathit{t}\left[\mathbf{s}\right]$	${\mathit{i}}_{\mathit{w}}\left[\mathbf{A}\right]$	${\mathit{\omega }}_{\mathit{s}}\left[\frac{\mathbf{rad}}{\mathbf{s}}\right]$
0.000000	0.000000	−0.000000
1.000000	4.951052	0.979779
2.000000	4.754254	0.596283
3.000000	−0.000552	−0.381824
4.000000	0.384783	−0.001621
5.000000	4.951041	0.972042
6.000000	4.384571	0.594716
7.000000	−0.000541	−0.374180
8.000000	0.739952	0.000653
9.000000	4.951030	0.964474
10.000000	4.950499	0.596505

). The main reason behind this might be identified as a much bigger number of samples utilised to generate a plot (one thousand time stamps), but the aspect of a fully analytical solution being more accurate than HAN, in general, cannot be omitted and will be analysed in future works. Nonetheless, it is safe to assume that even if the HAN approach with reduced step size would result in results indistinguishable from the analytical solution, it will still require more preparation, resulting in higher overall time cost. Since the HAN solution provides the possibility to compare the behaviour of the numerical solution regarding changes in step size and a comparison with commercial software in this aspect, it was selected for further consideration in the study. The integration step is $0.1$ $\mathrm{s}$, and selected time stamps of 1 $\mathrm{s}$ interval are presented in

Table 3. Results of the simulation conducted using the DynPy library—HAN method.

$\mathit{t}\left[\mathbf{s}\right]$	${\mathit{i}}_{\mathit{w}}\left[\mathbf{A}\right]$	${\mathit{\omega }}_{\mathit{s}}\left[\frac{\mathbf{rad}}{\mathbf{s}}\right]$
0.000000	−0.000000	0.000000
1.000000	4.951055	0.981733
2.000000	4.970133	0.596505
3.000000	0.019241	−0.387659
4.000000	0.000164	−0.002446
5.000000	4.951056	0.981718
6.000000	4.970133	0.596505
7.000000	0.019241	−0.387659
8.000000	0.000164	−0.002446
9.000000	4.951056	0.981718
10.000000	4.970133	0.596505

.

3.3. DC Motor—Conducting Simulations Using MATLAB/Simulink

MATLAB/Simulink, as a powerful tool for the simulation and modelling of dynamic systems, has found wide application in various engineering fields. In the context of this paper, it was used as a reference tool for comparison with the results obtained from the DynPy library, since its widespread use in the fields of control systems, drive systems, automation, and many other areas of engineering makes conducted simulations a reliable reference point.

Reference simulations were carried out using both variable and fixed step solvers. This approach enables a detailed analysis of the impact of different integration strategies on simulation results and allows for a direct comparison with the Python implementation, as the DynPy library HAN solution utilises a fixed step approach. The designed model of the DC motor used for conducting simulations in the MATLAB/Simulink environment is depicted in Figure 18.

The results obtained using the variable integration step are presented in

Table 4. Results of the simulation conducted with variable integration step in MATLAB/Simulink.

$\mathit{t}\left[\mathbf{s}\right]$	${\mathit{i}}_{\mathit{w}}\left[\mathbf{A}\right]$	${\mathit{\omega }}_{\mathit{s}}\left[\frac{\mathbf{rad}}{\mathbf{s}}\right]$
0.000000	0.000000	0.000000
1.000000	4.950632	0.981745
2.000000	4.970501	0.596487
3.000000	0.019819	−0.387680
4.000000	−0.000201	−0.002427
5.000000	4.950692	0.981726
6.000000	4.970500	0.596487
7.000000	0.019819	−0.387680
8.000000	−0.000201	−0.002427
9.000000	4.950692	0.981726
10.000000	4.970500	0.596487

(Figure 19 and Figure 20), while those with a fixed step set to be the same as in the simulations using the DynPy library ( $0.1$ $\mathrm{s}$) are in in

Table 5. Results of the simulation conducted with a constant integration step of $0.1$ $\mathrm{s}$ in MATLAB/Simulink.

$\mathit{t}\left[\mathbf{s}\right]$	${\mathit{i}}_{\mathit{w}}\left[\mathbf{A}\right]$	${\mathit{\omega }}_{\mathit{s}}\left[\frac{\mathbf{rad}}{\mathbf{s}}\right]$
0.000000	0.000000	0.000000
1.000000	4.950963	0.982077
2.000000	4.970140	0.596408
3.000000	0.019333	−0.388004
4.000000	0.000157	−0.002349
5.000000	4.950964	0.982063
6.000000	4.970140	0.596408
7.000000	0.019333	−0.388004
8.000000	0.000157	−0.002349
9.000000	4.950964	0.982063
10.000000	4.970140	0.596408

(Figure 21 and Figure 22).

After initial analysis of the simulation results obtained with variable and constant integration steps, significant differences were noticeable both in the achieved values and in the shape of the current intensity in the rotor windings and the rotor angular velocity. Therefore, a decision was made to reduce the constant integration step to a value of $0.01$ $\mathrm{s}$ and conduct further simulations. In Figure 23 and Figure 24 are the results obtained after implementing this modification. Values of selected time stamps are listed in

Table 6. Results of the simulation conducted with a specified constant integration step of $0.01$ $\mathrm{s}$ in MATLAB/Simulink.

$\mathit{t}\left[\mathbf{s}\right]$	${\mathit{i}}_{\mathit{w}}\left[\mathbf{A}\right]$	${\mathit{\omega }}_{\mathit{s}}\left[\frac{\mathbf{rad}}{\mathbf{s}}\right]$
0.000000	0.000000	0.000000
1.000000	4.951055	0.981734
2.000000	4.970133	0.596505
3.000000	0.019241	−0.387659
4.000000	0.000164	−0.002446
5.000000	4.951056	0.981718
6.000000	4.970133	0.596505
7.000000	0.019241	−0.387659
8.000000	0.000164	−0.002446
9.000000	4.951056	0.981718
10.000000	4.970133	0.596505

.

3.4. Application of DynPy for Electrical Circuits

DynPy library enables to modelling of electric circuits in a similar way as it handles electro-mechanical and mechanical systems. Electrical modelling was performed based on a serial DC RLC circuit (Figure 25). A single-degree-of-freedom, second-order system was introduced, in which the Heaviside function provides the variable value of the resistance. Such a definition of equations of motion introduces a more challenging problem, as the nonconstant parameter directly affects the behaviour of the first derivative of the dependent variable, while in the case of the DC motor system, the forcing vector was subjected to periodic changes:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill -U+L\ddot{q}+R_{min}\dot{q}+{\displaystyle \frac{q}{C}}+R\theta \left(sin\left({\displaystyle \frac{2\pi t}{T}}\right)\right)\dot{q}=0\end{array}\end{array}$$

(17)

where the following are used:

• R—amplitude of the variable part of the resistance;
• $R_{min}$—constant value of the resistance;
• L—inductance;
• C—capacitance;
• U—source voltage;
• q—electric charge;
• $\dot{q}$—current;
• $\ddot{q}$—current rate of change.

The simulational model adopts the parameters for modulation control of the current in a power circuit. Values of the system’s parameters are provided in

Table 7. Parameter values adopted for simulations of RLC circuit.

Parameter	Unit	Value
U	V	10.0
R	$\Omega$	19.9
$R_{min}$	$\Omega$	0.1
L	H	1
C	F	100
T	s	4

. Parameters were selected in a way that allows for presenting analysed phenomena in an understandable form. The period was selected arbitrarily, then other parameters were calculated to ensure decay times lower than the period. Under these conditions, the circuit can operate as a current control system that can be utilised to simulate electronic commutation processes or polarisation changes in direct current machines. The resistance signal is shown in Figure 26 and described with the following formula:

$$R\left(t\right)=\left\{\begin{array}{cc}0.1\hfill & \mathrm{for}t=(0,2)\cup (4,6)\cup (8,10)\hfill \\ 20\hfill & \mathrm{for}t=(2,4)\cup (6,8)\hfill \end{array}\right.$$

(18)

To solve the equation, the HAN method was utilised in a similar way as in the case of the DC motor. First, the solution (19) and its first derivative representing current (20) for a system with stationary resistance was returned from an instance of ODEs system class:

$$\begin{array}{cc}\hfill q=& C_1{e}^{t\left(-{\displaystyle \frac{R}{2L}}-{\displaystyle \frac{\sqrt{C{R}^2-4L}}{2\sqrt{C}L}}\right)}+C_2{e}^{t\left(-{\displaystyle \frac{R}{2L}}+{\displaystyle \frac{\sqrt{C{R}^2-4L}}{2\sqrt{C}L}}\right)}-{\displaystyle \frac{2CLU}{\sqrt{C}R\sqrt{C{R}^2-4L}+C{R}^2-4L}}\hfill \\ & +{\displaystyle \frac{2CLU}{\sqrt{C}R\sqrt{C{R}^2-4L}-C{R}^2+4L}}\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill \dot{q}=& C_1\left(-{\displaystyle \frac{R}{2L}}-{\displaystyle \frac{\sqrt{C{R}^2-4L}}{2\sqrt{C}L}}\right)e^{t\left(-{\displaystyle \frac{R}{2L}}-{\displaystyle \frac{\sqrt{C{R}^2-4L}}{2\sqrt{C}L}}\right)}\hfill \\ \hfill & +C_2\left(-{\displaystyle \frac{R}{2L}}+{\displaystyle \frac{\sqrt{C{R}^2-4L}}{2\sqrt{C}L}}\right)e^{t\left(-{\displaystyle \frac{R}{2L}}+{\displaystyle \frac{\sqrt{C{R}^2-4L}}{2\sqrt{C}L}}\right)}+{\displaystyle \frac{CRU}{\sqrt{C}R\sqrt{C{R}^2-4L}+C{R}^2-4L}}\hfill \\ \hfill & +{\displaystyle \frac{\sqrt{C}U\sqrt{C{R}^2-4L}}{\sqrt{C}R\sqrt{C{R}^2-4L}+C{R}^2-4L}}+{\displaystyle \frac{\sqrt{C}U\sqrt{C{R}^2-4L}}{\sqrt{C}R\sqrt{C{R}^2-4L}-C{R}^2+4L}}\hfill \\ \hfill & -{\displaystyle \frac{CRU}{\sqrt{C}R\sqrt{C{R}^2-4L}-C{R}^2+4L}}\hfill \end{array}$$

(20)

Formulas (19) and (20) were later introduced into a slightly modified (inputted system, data dictionaries, and number of time stamps) han_solution_loop function, presented in Appendix D. The HAN method was used twice to obtain results for $0.1$ $\mathrm{s}$ and $0.01$ $\mathrm{s}$ time steps. Simultaneously, the system was created in the MATLAB/Simulink environment (Figure 27) and solved there separately with variable, $0.1$ $\mathrm{s}$, and $0.01$ $\mathrm{s}$ step sizes.

The considered case underwent the entire process of electrical circuit modelling and its further analysis. The MATLAB/Simulink environment was used again as a reference tool for the assessment of the performance of the demonstrated library. A detailed view of the comparison between RLC circuit results generated with DynPy and MATLAB/Simulink is provided in Section 4.2. The results presented there confirm the outcomes obtained previously with the analysis of the DC motor.

4. Results and Discussion

4.1. DC Motor

The results of computer simulations were subjected to comparative analysis to assess the effectiveness of the developed library. The listed plots and tables show a comparison of the armature current and rotor angular velocity of the DC motor obtained using DynPy (with a constant integration step size of $0.1$ $\mathrm{s}$) and the following:

• MATLAB/Simulink with a variable integration step size (Figure 28,

  Table 8. DC motor: comparison of results for chosen time stamps—MATLAB/Simulink with variable step and HAN method with $0.1$ $\mathrm{s}$ step.

  	MATLAB, $\mathbf{\Delta }\mathit{t}\ne \mathit{const}$		DynPy HAN, $\mathbf{\Delta }\mathit{t}=0.1\mathbf{s}$
  $\mathit{t}\left[\mathbf{s}\right]$	${\mathit{i}}_{\mathit{w}}\left[\mathbf{A}\right]$	${\mathit{\omega }}_{\mathit{s}}\left[{\displaystyle \frac{\mathbf{rad}}{\mathbf{s}}}\right]$	${\mathit{i}}_{\mathit{w}}\left[\mathbf{A}\right]$	${\mathit{\omega }}_{\mathit{s}}\left[{\displaystyle \frac{\mathbf{rad}}{\mathbf{s}}}\right]$
  0.100000	4.303243	0.246351	4.317466	0.236242
  1.100000	4.954424	0.889913	4.955557	0.827837
  2.100000	0.662744	0.345930	0.652732	0.359291
  3.100000	0.015932	−0.295513	0.014739	−0.233768
  4.100000	4.242498	0.233902	4.317565	0.234769
  5.100000	4.954449	0.890017	4.955558	0.827828
  6.100000	0.662744	0.345930	0.652732	0.359291
  7.100000	0.015932	−0.295513	0.014739	−0.233768
  8.100000	4.242498	0.233902	4.317565	0.234769
  9.100000	4.954449	0.890017	4.955558	0.827828

  );
• MATLAB/Simulink with a constant step of $0.1$ $\mathrm{s}$ (Figure 29,

  Table 9. DC motor: comparison of results for chosen time stamps—MATLAB/Simulink with $0.1$ $\mathrm{s}$ step and HAN method $0.1$ $\mathrm{s}$ step.

  	MATLAB, $\mathbf{\Delta }\mathit{t}=0.01\mathbf{s}$		DynPy HAN, $\mathbf{\Delta }\mathit{t}=0.1\mathbf{s}$
  $\mathit{t}\left[\mathbf{s}\right]$	${\mathit{i}}_{\mathit{w}}\left[\mathbf{A}\right]$	${\mathit{\omega }}_{\mathit{s}}\left[{\displaystyle \frac{\mathbf{rad}}{\mathbf{s}}}\right]$	${\mathit{i}}_{\mathit{w}}\left[\mathbf{A}\right]$	${\mathit{\omega }}_{\mathit{s}}\left[{\displaystyle \frac{\mathbf{rad}}{\mathbf{s}}}\right]$
  0.100000	6.650000	0.083333	4.317466	0.236242
  1.100000	4.952507	0.827282	4.955557	0.827837
  2.100000	−1.679797	0.512135	0.652732	0.359291
  3.100000	0.017789	−0.233214	0.014739	−0.233768
  4.100000	6.650094	0.081924	4.317565	0.234769
  5.100000	4.952508	0.827273	4.955558	0.827828
  6.100000	−1.679797	0.512135	0.652732	0.359291
  7.100000	0.017789	−0.233214	0.014739	−0.233768
  8.100000	6.650094	0.081924	4.317565	0.234769
  9.100000	4.952508	0.827273	4.955558	0.827828

  );
• MATLAB/Simulink with a constant step of $0.01$ $\mathrm{s}$ (Figure 30,

  Table 10. DC motor: comparison of results for chosen time stamps—MATLAB/Simulink with $0.01$ $\mathrm{s}$ step and HAN method with $0.1$ $\mathrm{s}$ step.

  	MATLAB, $\mathbf{\Delta }\mathit{t}=0.01\mathbf{s}$		DynPy HAN, $\mathbf{\Delta }\mathit{t}=0.1\mathbf{s}$
  $\mathit{t}\left[\mathbf{s}\right]$	${\mathit{i}}_{\mathit{w}}\left[\mathbf{A}\right]$	${\mathit{\omega }}_{\mathit{s}}\left[{\displaystyle \frac{\mathbf{rad}}{\mathbf{s}}}\right]$	${\mathit{i}}_{\mathit{w}}\left[\mathbf{A}\right]$	${\mathit{\omega }}_{\mathit{s}}\left[{\displaystyle \frac{\mathbf{rad}}{\mathbf{s}}}\right]$
  0.100000	4.317996	0.236209	4.317466	0.236242
  1.100000	4.955557	0.827836	4.955557	0.827837
  2.100000	0.652203	0.359324	0.652732	0.359291
  3.100000	0.014740	−0.233768	0.014739	−0.233768
  4.100000	4.318094	0.234735	4.317565	0.234769
  5.100000	4.955557	0.827827	4.955558	0.827828
  6.100000	0.652203	0.359324	0.652732	0.359291
  7.100000	0.014740	−0.233768	0.014739	−0.233768
  8.100000	4.318094	0.234735	4.317565	0.234769
  9.100000	4.955557	0.827827	4.955558	0.827828

  ).

The simulation results obtained using the DynPy library with a constant integration step size of $0.1$ $\mathrm{s}$ and in MATLAB/Simulink using a variable integration step size confirm the agreement between the results outputted from both tools. To present the compatibility of the results more clearly, the characteristics of both simulations are overlaid in (Figure 28). Only minor differences are noticeable in the rotor current plots in areas of dynamic changes in characteristics. The results of the simulation conducted in MATLAB/Simulink are more precise, and the curve representing them is smoother. These differences result from the use of a variable integration step size in the simulation using MATLAB/Simulink. This approach allows for adjusting the integration step size depending on the rate of change in the function, resulting in better accuracy of results in areas of dynamic changes. In the case of rotor angular velocity characteristics, there are no significant differences—the plots are almost identical.

Next, the simulation results were compared, where in both software, a constant integration step of $0.1$ $\mathrm{s}$ was specified (Figure 29). Significant differences in the characteristics of rotor current intensity are visible. In MATLAB/Simulink, the waveform of the rotor current intensity undergoes oscillations in regions of rapid changes in its value, while the corresponding plot generated using the DynPy library remains stable. In the same intervals as the rotor current intensity, there are discrepancies in the characteristics of angular velocity. However, they are less noticeable, and their differences are acceptable. The disagreement in the results may stem from various factors such as the numerical implementation of the algorithm in both environments or the handling of numerical errors.

Due to discrepancies in the simulation results with the specified constant integration step of $0.1$ $\mathrm{s}$, it was decided to reduce its value to $0.01$ $\mathrm{s}$ in the case of MATLAB/Simulink software and conduct further simulations. The juxtaposed characteristics confirm the consistency of the simulation results, and for better presentation, the characteristics are overlaid (Figure 30). The graph is very similar to the plot in Figure 28, with only minor, negligible differences in the values of rotor current intensity and rotor angular velocity in the areas of rapid changes.

4.2. RLC Circuit

Similarly, the comparison was prepared for the RLC circuit:

• With both solvers utilizing constant step of $0.1$ $\mathrm{s}$ (Figure 31,

  Table 11. RLC circuit: comparison of results for chosen time stamps—MATLAB/Simulink with $0.1$ $\mathrm{s}$ step and HAN method with $0.1$ $\mathrm{s}$ step.

  	DynPy, HAN, $\mathbf{\Delta }\mathit{t}=0.1\mathbf{s}$	MATLAB/Simulink, $\mathbf{\Delta }\mathit{t}=0.1\mathbf{s}$
  $\mathit{t}\left[\mathbf{s}\right]$	$\dot{\mathit{q}}$[A]	$\dot{\mathit{q}}$[A]
  0.000000	0.000000	0.000000
  1.000000	0.502298	0.497255
  2.000000	0.502046	0.497013
  3.000000	9.494228	10.368359
  4.000000	17.585402	20.136108
  5.000000	0.491754	0.486170
  6.000000	0.491507	0.485625
  7.000000	9.294929	10.130792
  8.000000	17.216256	19.674735
  9.000000	0.481431	0.475031
  10.000000	0.481189	0.474498

  );
• With both solvers utilizing constant step of $0.01$ $\mathrm{s}$ (Figure 32,

  Table 12. RLC circuit: comparison of results for chosen time stamps—MATLAB/Simulink with $0.01$ $\mathrm{s}$ step and HAN method with $0.01$ $\mathrm{s}$ step.

  	DynPy, HAN, $\mathbf{\Delta }\mathit{t}=0.01\mathbf{s}$	MATLAB/Simulink, $\mathbf{\Delta }\mathit{t}=0.1\mathbf{s}$
  $\mathit{t}\left[\mathbf{s}\right]$	$\dot{\mathit{q}}$[A]	$\dot{\mathit{q}}$[A]
  0.000000	0.000000	0.000000
  1.000000	0.502287	0.497263
  2.000000	0.502034	0.497013
  3.000000	9.898286	10.368359
  4.000000	18.310948	20.136108
  5.000000	0.491532	0.485869
  6.000000	0.491285	0.485625
  7.000000	9.686345	10.130792
  8.000000	17.918876	19.674735
  9.000000	0.481007	0.474736
  10.000000	0.480765	0.474498

  );
• With a variable integration step size in MATLAB/Simulink, $0.1$ $\mathrm{s}$ and $0.01$ $\mathrm{s}$ step size in DynPy (Figure 33,

  Table 13. RLC circuit: comparison of results for chosen time stamps—MATLAB/Simulink with variable step and HAN method with $0.1$ $\mathrm{s}$ and $0.01$ $\mathrm{s}$ step.

  	DynPy, $\mathbf{\Delta }\mathit{t}=0.1\mathbf{s}$	DynPy, $\mathbf{\Delta }\mathit{t}=0.01\mathbf{s}$	MATLAB/Simulink, $\mathbf{\Delta }\mathit{t}\ne \mathit{const}$
  $\mathit{t}\left[\mathbf{s}\right]$	$\dot{\mathit{q}}$[A]	$\dot{\mathit{q}}$[A]	$\dot{\mathit{q}}$[A]
  0.000000	0.000000	0.000000	0.000000
  1.000000	0.502298	0.502287	0.497267
  2.000000	0.502046	0.502034	0.497013
  3.000000	9.494228	9.898286	10.295847
  4.000000	17.585402	18.310948	20.144906
  5.000000	0.491754	0.491532	0.485902
  6.000000	0.491507	0.491285	0.485620
  7.000000	9.294929	9.686345	10.099166
  8.000000	17.216256	17.918876	19.692078
  9.000000	0.481431	0.481007	0.474862
  10.000000	0.481189	0.480765	0.474468

  ).

When considering the comparison of simulations for a $0.1$ $\mathrm{s}$ step size, both qualitative and quantitative discrepancies are clearly visible (Figure 31). Above all, the curve generated in MATLAB/Simulink is smooth, while in the case of HAN, sudden changes in current value appear in the second, third, sixth, and seventh second. These changes have a periodic character and are directly related to the implementation of the method. The switch of datasets is realised every second, so to ensure the continuity of the solution, new initial values are introduced to its consecutive parts. As these initial conditions are obtained directly from the previous part of the solution for the same time stamp at which the next part of the solution starts, vertical shifts appear in the plot. A similar situation is expected to take place at every dataset change, yet might not be visible due to the different slopes of the plot at other time stamps. On the other hand, the maximum amperage in the fourth and eighth seconds is higher in MATLAB/Simulink. The fact that the HAN method predominantly utilises an analytical approach lets us conclude that in this aspect, DynPy is more reliable.

The plot in Figure 32 compares the results from both environments with the step size set to $0.01$ $\mathrm{s}$. In this scenario, only qualitative difference is exposed—both graphs are equally smooth and generally take the same shape. As previously, maximum current values at the fourth and eighth second are smaller in comparison with MATLAB/Simulink.

Finally, Figure 33 shows the comparison between both constant step size HAN solutions and variable step size output from MATLAB/Simulink. The graph obtained from MATLAB/Simulink performs marginally better qualitatively than the $0.1$ $\mathrm{s}$ HAN curve. In both cases, small, general waviness appears, as well as sudden rises in amperage at the third and seventh seconds (which also happen at the second and sixth seconds in the case of the $0.1$ $\mathrm{s}$ step size HAN solution). The $0.01$ $\mathrm{s}$ step size HAN method utilised in DynPy performs the best from the compared approaches, ensuring the highest smoothness and reliability of results.

Regardless of the selected step size, from a quantitative point of view, the hybrid analytical–numerical method utilised with the DynPy library showed better performance than the strictly numerical approach. When considering the step size, when selected properly, it also ensures high smoothness, resulting in an overall high quality of the obtained results.

5. Conclusions

This paper presents a general investigation of the DynPy library and its application for scientific and engineering activities. An overview of the developed tool was presented, with deeper insight into its modules and functionalities responsible for the provision of analytical and numerical solutions to systems described with linear and nonlinear ordinary differential equations. Nonlinear (more realistic representation of physical systems) and linear (idealised for theoretical or simplified analysis) models can be adopted and solved in an efficient way.

The solution for a model under consideration obtained with the DynPy library was presented in hybrid analytical–numerical (HAN) and fully analytical forms, where in the case of the latter, Heaviside forcing inputs were included. The HAN methodology, which enables the possibility of setting the integration step size, was chosen for comparison with MATLAB/Simulink results. The results obtained from the HAN method were compared with numerical simulations conducted in MATLAB/Simulink environment.

The analysis of results for the DC motor case leads to the conclusion that for cases with a variable integration step and a constant step of $0.01$ $\mathrm{s}$ in MATLAB/Simulink, the obtained results are consistent with the results obtained using the DynPy library with the specified constant step of $0.1$ $\mathrm{s}$. This confirms the correctness of the implementation in the developed library. However, in the case of a constant step of $0.1$ $\mathrm{s}$ in MATLAB/Simulink, significant quantitative and qualitative discrepancies were observed between the results—both in terms of values and waveform shapes. Deviating values indicate that under specific conditions, the DynPy library may respond better to changes in system dynamics than the MATLAB/Simulink tool. Specifically, the hybrid analytical–numerical method not only proved to return proper results in opposition to simulation with the same step size in MATLAB/Simulink but also presented high accuracy with $0.01$ $\mathrm{s}$ step simulations in the stated reference environment.

Summarizing the simulations of the RLC circuit, the main conclusion is that regardless of the selected step size, thanks to the solution being based on an analytical approach, the results are quantitatively better than the strictly numerical ones obtained from each analysed MATLAB/Simulink case. Qualitative discrepancies are different between approaches subjected to comparison. It is clearly visible that the DynPy HAN solution with a properly set step size provides smoother curves than MATLAB/Simulink with variable steps. Hence, the analytical solution updated numerically with appropriate frequency provides the most reliable results.

The comparisons of simulation results indicate that the numerical algorithms in the DynPy library have been implemented correctly, and the results obtained are accurate. The paper itself was generated strictly with the use of DynPy’s functionalities—all steps from building the model, through finding an analytical solution, conducting the simulations and their visualisation, to generating editable .tex files were carried out in the same environment. Thus, not only the accuracy and reliability of the tool were proven but its versatility and usefulness as well. Conclusions drawn from this study will enable further optimisation and development of the library for various engineering applications. Some of the future research directions include the development of a more user-friendly interface for hybrid analytical–numerical methods, the expansion of built-in analytical solvers with more types of step functions, and general refining and maintenance. It is worth mentioning that the library undergoes constant development, and numerous new features and fixes were added during and after the writing of this article.

By no means do the authors of this article imply the overall superiority of the DynPy library over well-known and field-tested software, but the main advantages have to be pointed out for comparison. Firstly, it is written in the Python programming language. Although it requires the user to obtain at least a basic knowledge of its structure and syntax, it also assures the fact that insight into developed algorithms and their working principles is possible, as well as local editing, in contrast to previously described software. Moreover, a large user community of the language, in general, ensures additional support on top of the one possible to obtain from DynPy’s built-in guides and functioning help methods. Thanks to the utilisation of the object-oriented paradigm, a unique way to build dynamic systems was implemented. Said paradigm also enables the user to create self-defined basic elements, without the need of downloading additional toolboxes or add-ins. Such an approach results in a low level of initial theoretical, mathematical, and modelling knowledge needed to use the tool. At no point when building the system, manual definition of equations of motion or governing equations is necessary (but still possible if needed), since they are derived directly from the user-provided composition of the model. Of course, verification of the obtained results is necessary, as with any other tool. Another advantage lies in the scope of solving possibilities. The same model can be solved both analytically (if possible from the mathematical point of view) and numerically in the same environment. Lastly, DynPy also works as a text editing tool. Thanks to the implementation of self-written and PyLaTeX inherited functionalities, fully formatted papers like this one are possible to generate without the need for workspace change, so the model preparation, solution finding, simulations, and reporting are all in one place.