This paper presents an advanced six-step linear multistep method of order five, based on second derivative non-hybrid block backward differentiation formula (BDF), developed for numerical solution of stiff systems of ordinary differential equations. The development of the multistep collocation approach is carried out using the matrix inversion techniques. The block structure facilitates simultaneous computation of multiple solution points, improving computational efficiency. Moreover, the power series is adopted as basis function for driving the discrete and continuous formulations. The analysis of the method such as consistency, zero-stability, order and error constants is presented, confirming its suitability for stiff systems. Numerical experiments on standard are tested on stiff and non-stiff ordinary differential equations showed the new method outperforms existing method in terms of accuracy.

Modern energy and electronic systems comprise a large number of nonlinear elements operating in dynamic modes. Accurate modeling of such systems requires proper consideration of nonlinear current-voltage characteristics (CVC) and transient processes, which is especially important when designing rectifiers, inverters, control systems and other power electronics devices. Existing electrical circuit simulators do not always provide users with the necessary flexibility, scalability, or compatibility with enterprise safety standards, and may have legal restrictions. Custom effective modeling methods for such systems allow for creating specialized software solutions to analyze dynamic modes of electrical circuits, free from the limitations of commercial simulators and tailored to specific engineering and scientific tasks. The purpose of the work is to develop a method for numerical modeling of dynamic modes of electrical circuits comprising semiconductor diodes or other elements with nonlinear CVC, described by nodal equations. The scientific novelty lies in the development of a homotopy approach to overcoming high condition number of the Jacobian matrix when solving nonlinear differential-algebraic equations (DAEs) of electrical circuits, based on deformation of CVC with adaptation of the deformation parameter depending on the residual norm and condition number of the matrix; in the development of a method for extracting linearly independent differential equations from DAEs of electrical circuits without preliminary analysis of their topology; in the development of a universal stamp of a nonlinear element, based on linearization of the functional in the vicinity of the current approximation, allowing for integration of diode models with various CVCs into nodal equations. Materials and methods. Theoretical electrical engineering methods were used in the work, including the modified nodal potential method. The proposed numerical modeling method involves integration of elements with nonlinear CVC into nodal equations; extraction of differential and algebraic parts from DAEs, based on singular value decomposition of the matrix standing before the derivative vector; transformation of DAEs into a nonlinear system of algebraic equations using backward differentiation formulas (BDF) with variable time step. Initial points for BDF were determined by the diagonally implicit Runge–Kutta method of second order accuracy. Numerical solution of the obtained nonlinear equations was performed by the damped Newton–Raphson method. To reduce the condition number of the Jacobian matrix in transient modes, when the spread of differential conductivities reaches 12 orders and higher, a homotopy approach was proposed, consisting of gradual deformation of the diode CVC from a smoothed to the original curve during convergence, while maintaining a given value of the condition number. Results. To demonstrate the proposed solutions, computer simulation of a bridge rectifier operating on an active-inductive load with two types of diode CVC was performed: piecewise-linear and smooth, corresponding to the Shockley equation with series resistance. The deformation parameter and damping coefficient were adaptively changed depending on the residual norm of the functional and the condition number of the Jacobian matrix. Comparison of simulation results with different methods of specifying diode CVC showed that differences appear predominantly in transient processes of switching diode operation modes. It has been found that to ensure convergence of numerical solution in diode switching modes, characterized by high condition number of the Jacobian matrix, the homotopy approach is more effective than diagonal regularization. The proposed method for numerical modeling of dynamic modes of electrical circuits with nonlinear elements has a natural algorithmic structure, allowing for simple software implementation. Conclusions. 1. The most universal diode stamp, obtained on the basis of linearization of the functional derived from the CVC equation in the vicinity of the current approximation, has been identified. 2. A method for extracting linearly independent differential equations from DAEs of electrical circuits without preliminary analysis of circuit topology has been proposed. 3. A method for calculating the Jacobian matrix for solving nonlinear DAE has been proposed. 4. To ensure convergence of numerical solution with high condition number of the Jacobian matrix, it is preferable to apply the homotopy approach.

Superconvergence analysis for nonlinear Ginzburg-Landau equation with backward differential formula finite element method

Numerical method for Fokker–Planck equations based on backward differentiation formulas

ABSTRACT This paper introduces a robust and efficient numerical framework for solving nonlinear time‐fractional partial integro‐differential equations (NLTFPIDEs). Temporal discretization is performed using the weighted and shifted Grünwald–Letnikov formula, incorporating the fractional trapezoidal rule and the second‐order backward differentiation formula (BDF2). Spatial discretization leverages Chebyshev nodes as discretization points, with the Lagrange‐collocation method applied to approximate partial derivatives. For irregular computational domains, the framework utilizes the finite block method (FBM) in two dimensions. Nonlinearities in the equations are handled through the quasilinearization technique. A comprehensive stability and convergence analysis using the energy method confirms the reliability of the proposed schemes. Numerical experiments validate the theoretical findings, highlighting the accuracy and efficiency of the approach.

Abstract In scientific machine learning (SciML), a key challenge is learning unknown, evolving physical processes and making predictions across spatio-temporal scales. For example, in real-world manufacturing problems like additive manufacturing, users adjust known machine settings while unknown environmental parameters simultaneously fluctuate. To make reliable predictions, it is desired for a model to not only capture long-range spatio-temporal interactions from data but also adapt to new and unknown environments; traditional machine learning models excel at the first task but often lack physical interpretability and struggle to generalize under varying environmental conditions. To tackle these challenges, we propose the attention-based spatio-temporal neural operator (ASNO), a novel architecture that combines separable attention mechanisms for spatial and temporal interactions and adapts to unseen physical parameters. Inspired by the backward differentiation formula, ASNO learns a transformer for temporal prediction and extrapolation and an attention-based neural operator for handling varying external loads, enhancing interpretability by isolating historical state contributions and external forces, enabling the discovery of underlying physical laws and generalizability to unseen physical environments. Empirical results on SciML benchmarks demonstrate that ASNO outperforms existing models, establishing its potential for engineering applications, physics discovery, and interpretable machine learning.

The Bratu's equation is a strongly nonlinear second-order ordinary differential equation that arises in electrospinning process and models temperature distribution within a flame in combustion theory. Bratu-type equations are used to simulate the ignition of flammable gases and flame propagation. In this paper, two methods are proposed to obtain highly accurate and reliable approximate solutions of Bratu-type boundary value problems. The first technique is a power series method which is based on the generalised Cauchy product that simplifies the difficulty associated with the nonlinear terms. Subsequently, explicit recurrence relations for the expansion coefficients of the series solutions are obtained. The second approach uses a twelfth-order second derivative backward differentiation formula that is implemented as a boundary value method. This numerical method is referred to as second derivative backward differentiation boundary value method. Three examples are given to illustrate the effectiveness, reliability, and accuracy of the proposed methods. The results obtained from both methods are in excellent agreement with the known exact solution. Comparison of the approximate and exact solutions shows that the proposed methods are reliable and accurate in solving a class of strongly nonlinear boundary value problems of Bratu-type.

ABSTRACT The time fractional cable equation (TFCE) extends the classical cable model to describe the anomalous diffusion of ionic motion in the nervous system. This paper presents an efficient numerical method for solving the TFCE, combining the orthogonal spline collocation (OSC) method with the alternating direction implicit (ADI) method for spatial discretization. For time discretization, it employs a second‐order backward differentiation formula (BDF2) along with the classical L1 approximation for the Caputo fractional derivative. The stability and convergence of the proposed method are rigorously analyzed. Finally, numerical experiments validate the theoretical results and demonstrate its accuracy and efficiency.

This study presents a systematic benchmarking of numerical methods for solving ordinary differential equations (ODEs) applied to damped single-degree-of-freedom (SDOF) vibration systems. Ten solvers—including Runge–Kutta variants, Adams–Bashforth–Moulton, Rosenbrock, and Backward Differentiation Formula (BDF)—were evaluated under both non-stiff and stiff conditions by varying mass, damping, and stiffness parameters. Analytical solutions were used as references to quantify global error, convergence behavior, and computational efficiency. High-order adaptive solvers, such as Verner’s and Runge–Kutta 7/8, consistently achieved the highest accuracy, reducing global error by up to 15% compared with the classical Runge–Kutta (ODE45) method. Implicit methods, including Rosenbrock and BDF, demonstrated superior stability in stiff and highly damped cases. In contrast, low-order approaches, particularly the Trapezoidal rule, exhibited the largest errors, exceeding 30% in oscillatory regimes. The results confirm that solver performance is problem-dependent, emphasizing that no single algorithm is universally optimal. Beyond the technical contributions, this study introduces a pedagogical framework that allows engineering students to visualize solver trade-offs, quantify numerical accuracy, and interpret computational efficiency. The educational integration strengthens conceptual understanding of numerical methods and supports data-driven solver selection in vibration analysis and related engineering applications.

Abstract In this paper, we employ an adaptive time-stepping strategy to conduct a comprehensive theoretical analysis and numerical simulations of the magnetohydrodynamic (MHD) system. Our objective is to formulate a numerical scheme based on the variable time-steps second-order backward difference formula (VBDF2), which exhibits linear, decoupled, and semi-discrete properties. Firstly, considering the complexity of handling the nonlinear terms in the MHD equations, we will introduce a method that combines the exponential scalar auxiliary variable with the zero-energy contribution. This allows us to construct a new auxiliary variable, which is of crucial importance in dealing with the nonlinear terms. Subsequently, by reformulating the MHD system using this new auxiliary variable, we have successfully derived an equivalent system that is unconditionally stable in terms of energy. Secondly, based on the adaptive time-stepping strategy, we establish a semi-discrete MHD system to guarantee the stability of the system. It is worth noting that, under the favorable condition of being unrestricted by the time-step size, we present a rigorous error analysis. Finally, we carry out several numerical simulations to verify the theoretical results by using an effective adaptive time-step strategy.

ABSTRACT In this paper, we study the numerical solution of semilinear parabolic equations with nonlinear vanishing delay, employing a linearized variable‐time‐step three‐step backward differentiation formula (vBDF3) for time discretization and the virtual element method (VEM) for spatial discretization. The proposed scheme is flexible for variable time steps and spatial meshes. Based on the technique of the discrete orthogonal convolution kernels, we establish the unconditional optimal ‐norm error estimate for the fully discrete solution through the error splitting technique. Finally, numerical experiments are given to demonstrate the high‐order accuracy and unconditional convergence of the linearized scheme.

First-order Ordinary Differential Equations (ODEs) are often characterized by stiffness, especially in models that describe complex real-world processes. This work presents a four-point, fixed-coefficient, diagonally implicit block backward differentiation formula (4BBDF) of second order, developed to address the numerical challenges associated with stiffness. The formulation incorporates a diagonal matrix into the Lagrange interpolation polynomial and is constructed using Maple to ensure accuracy and stability. Newton’s method is used to handle the nonlinear systems that arise. The proposed 4BBDF method is mathematically verified to be consistent, zero-stable, A-stable, and of second-order accuracy. Its implementation in C++ shows improved computational efficiency, reducing the number of steps required by approximately 50% when compared to existing methods. These results indicate that the proposed scheme is a reliable and effective tool for solving stiffness in ODE.

In this paper, we propose a linear stabilized first-order backward differentiation formula scheme (LsBDF1) and a linear stabilized Crank–Nicolson scheme (LsCN) with nonuniform time steps, which together with two additional stabilization terms, preserve energy stability and the maximum bound principle (MBP) for ternary Allen–Cahn equations. Under mild constraints on the time steps and two stabilizing parameters, we prove that the LsBDF1 scheme and the LsCN scheme preserve the discrete MBP and unconditional energy stability. Furthermore, the optimal L ∞ error estimates are rigorously analysed. Various numerical experiments are conducted to verify the theoretical results.

Most physical systems exhibit nonlinear behavior while in motion, making their resolution challenging due to nonlinearity, dynamic effects, and sensitivity to parameters such as frequency and amplitude. Traditional analytical and numerical approaches can address these challenges but offer high computational costs, particularly in solving the system of free vibrations produced by the tapered beam. Predicting the behavior of this model is complicated, due to its high sensitivity and nonlinearity. Previously, standard neural network models have been used to solve dynamical systems, but they lack efficiency in handling nonlinearity. In this paper, we propose a novel deep learning model that predicts the amplitude of vibrations of a tapered beam. The primary focus of this study is to address the nonlinearity of the model and accurately predict the amplitude of vibrations. To solve this issue, we introduce a deep neural network designed to manage both nonlinearity and dynamical effects, including amplitude. The approach is significant in terms of computational and time efficiency compared to traditional numerical methods. The proposed work provides comparative results generated by the deep neural network, the backward difference formula as an analytical technique, and the Adams–Bashforth–Moulton predictor–corrector method as a numerical approach. The results demonstrate that our model outperforms existing numerical and analytical techniques. With the help of mean square error, Thiel’s inequality coefficient, and mean absolute error, the accuracy of our model can be verified; the lower these values, the more accurate our model will be. In our proposed model, the values are 8.389× 10−9 for mean square error, 5.563×10−4 for Thiel’s inequality coefficient, and 0.347 for mean absolute error; all these values are close to zero, signifying the accuracy of our model. The conclusion confirms that our proposed approach, even with changeable hyperparameters, is more suitable and accurate than numerical and analytical methods.

ABSTRACT We investigate the stability of two families of three‐level two‐step schemes that extend the classical second‐order backward differentiation formula (BDF2) and second‐order Adams–Moulton (AM2) schemes. For a free parameter restricted to an appropriate range that covers the classical case, we show that both the generalized BDF2 and the generalized AM2 schemes are A‐stable. We also introduce the concept of uniform‐in‐time stability which characterizes a scheme's ability to inherit the uniform boundedness over all time of the solution of damped and forced equation with the force uniformly bounded in time. We then demonstrate that A‐stability and uniform‐in‐time stability are equivalent for an arbitrary three‐level two‐step schemes. Next, these two families of schemes are utilized to construct efficient and unconditionally stable IMEX schemes for systems that involve a damping term, a skew symmetric term, and a forcing term. These novel IMEX schemes are shown to be uniform‐in‐time energy stable in the sense that the norm of any numerical solution is bounded uniformly over all time, provided that the forcing term is uniformly bounded time, the skew symmetric term is dominated by the dissipative term, together with a mild time‐step restriction. Numerical experiments verify our theoretical results. They also indicate that the generalized schemes could be more accurate and/or more stable than the classical ones for suitable choice of the parameter.

A new 2-point diagonally implicit variable step size super class of block backward differentiation formula (2DVSSBBDF) for solving first order stiff initial value problems (IVPs) is developed. The method is derived by introducing a lower triangular matrix in the coefficient matrix of existing 2-point variable step size superclass of block backward differentiation formula for the integration of stiff IVPs. The order of the method is 4. The stability analysis indicates that the method is both zero and A-stable. The Numerical results obtained are compared with some existing built in Matlab ODEs solvers in particular ODE15s and ODE23s and the performance of the new scheme showed an advantage in accuracy and computation time over some existing algorithms. The new method can serve as an alternative and efficient method for solving stiff IVPs.

ABSTRACTThis paper is concerned with efficient and accurate numerical schemes for the Cahn‐Hilliard‐Navier‐Stokes phase field model of binary immiscible fluids. By introducing two Lagrange multipliers for each of the Cahn‐Hilliard and Navier‐Stokes parts, we reformulate the original model problem into an equivalent system that incorporates the energy evolution process. Such a nonlinear, coupled system is then discretized in time using first‐ and second‐order backward differentiation formulas, in which all nonlinear terms are treated explicitly and no extra stabilization term is imposed. The proposed dynamically regularized Lagrange multiplier (DRLM) schemes are mass‐conserving and unconditionally energy‐stable with respect to the original variables. In addition, the schemes are fully decoupled: Each time step involves solving two biharmonic‐type equations and two generalized linear Stokes systems, together with two nonlinear algebraic equations for the Lagrange multipliers. A key feature of the DRLM schemes is the introduction of the regularization parameters which ensure the unique determination of the Lagrange multipliers and mitigate the time step size constraint without affecting the accuracy of the numerical solution, especially when the interfacial width is small. Various numerical experiments are presented to illustrate the accuracy and robustness of the proposed DRLM schemes in terms of convergence, mass conservation, and energy stability.

Numerical approximations of shape gradients and their applications have recently caused much interest in the shape optimization community. In this paper, existing research results on shape optimization governed by steady problems (e.g., elliptic problems in Hiptmair et al. [BIT Numer. Math. 55 (2015) 459–485]) are extended to those governed by parabolic problems. Convergence analysis is presented for numerical approximations of shape gradients associated with a parabolic problem. Both the backward Euler scheme and the backward differentiation formula are employed for time discretization, and the Galerkin finite element method is used for spatial discretization of the parabolic state and adjoint problems. The error of the distributed shape gradient is shown to have a higher convergence order in the mesh-size than that of the boundary type. A priori error estimates with respect to the time step-size are also presented. Numerical examples are provided to illustrate the theoretical results.

ABSTRACTA fully discrete implicit scheme is proposed for the Swift‐Hohenberg model, combining the third‐order backward differentiation formula (BDF3) for the time discretization and the second‐order finite difference scheme for the space discretization. Applying the Brouwer fixed‐point theorem and the positive definiteness of the convolution coefficients of BDF3, the presented numerical algorithm is proved to be uniquely solvable and unconditionally energy stable. Further, the numerical solution is shown to be bounded in the maximum norm. The proposed scheme is rigorously proven to be convergent in the norm by the discrete orthogonal convolution kernel, which transforms the four‐level form into the three‐level gradient form for the approximation of the temporal derivative. Consequently, the error estimate for the numerical solution is established by utilization of the discrete Grönwall inequality. Numerical examples in 2D and 3D cases are provided to support the theoretical results.

Abstract In this paper, we successfully propose a highly efficient linear, decoupled, fully discrete, variable time step second‐order backward difference formula (VBDF2) finite element numerical scheme for the Cahn–Hilliard–Navier–Stokes (CHNS) model. First, by utilizing the invariant energy quadratization (IEQ) method, we derive the equivalent CHNS model, which is crucial for handling the nonlinear terms. Second, by employing the finite element method for spatial discretization and the VBDF2 approach for time discretization, we obtain the fully discrete numerical scheme. Then, we also demonstrate the unique solvability, mass conservation property, and stability of the IEQ‐VBDF2 scheme. Subsequently, we present the error estimates accompanied by rigorous theoretical derivations. Finally, we creatively design an adaptive variable‐time‐step strategy and conduct several numerical simulations to verify the theoretical results. This strategy boasts excellent performance. It can maintain a high level of computational accuracy while significantly enhancing computational efficiency, offering a more efficient and reliable solution for the relevant numerical computing field.