CN114154110B - A numerical integration method for solving initial value problems of stiff ordinary differential equations in simulation - Google Patents
A numerical integration method for solving initial value problems of stiff ordinary differential equations in simulation Download PDFInfo
- Publication number
- CN114154110B CN114154110B CN202111479747.6A CN202111479747A CN114154110B CN 114154110 B CN114154110 B CN 114154110B CN 202111479747 A CN202111479747 A CN 202111479747A CN 114154110 B CN114154110 B CN 114154110B
- Authority
- CN
- China
- Prior art keywords
- algorithm
- solution
- solving
- integral
- iteration
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Classifications
-
- G—PHYSICS
- G06—COMPUTING OR CALCULATING; COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
- G06F17/13—Differential equations
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- Computational Mathematics (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Pure & Applied Mathematics (AREA)
- Data Mining & Analysis (AREA)
- Theoretical Computer Science (AREA)
- Operations Research (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- General Engineering & Computer Science (AREA)
- Algebra (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明公开了一种适用于刚性常微分方程初值问题的数值积分求解方法，在保证求解精度基础上提出了d‑LM方法，进一步提高了隐式运算稳定性。在填补国内空白的基础上，其运算效率接近国外同类产品；与国外商用系统仿真软件工具相比，本发明的数值方法更新颖，对1990年代之后的新的数值求解方法进行了实现与应用，数值求解过程更加稳定。The present invention discloses a numerical integration solution method suitable for initial value problems of rigid ordinary differential equations. On the basis of ensuring the solution accuracy, the d-LM method is proposed, which further improves the implicit operation stability. On the basis of filling the domestic gap, its operation efficiency is close to that of similar foreign products; compared with foreign commercial system simulation software tools, the numerical method of the present invention is more novel, and the new numerical solution method after the 1990s is realized and applied, and the numerical solution process is more stable.
Description
Technical Field
The invention belongs to the technical field of simulation, and relates to a numerical integration solving method for an initial value problem of a rigid ordinary differential equation in simulation.
Background
Many phenomena in engineering technology, such as operation of a control system, operation of an electric power system, motion of an aircraft, a process of chemical reaction and the like, are initial value problems (IVP, initialValue Problem) of ordinary differential equations (ODE, ordinary Differential Equation) in mathematical models, and many partial differential equations can be approximately simplified into initial value problems of ordinary differential equations to solve. Generally, the numerical solution of ordinary differential equations is relatively mature, and the theory is also perfect, with many mature solution algorithms being available.
However, there are a class of ordinary differential equations that suffer considerable difficulties in their solution, and the components of such ordinary differential equation solutions vary very rapidly and very slowly. From the point of view of numerical solution, a smaller integration step is used when the solution changes rapidly, and a larger integration step can be used when the component which changes rapidly tends to be stable. However, a great deal of theory and practice prove that no matter whether the solution changes steadily, the integral step length cannot be amplified, otherwise, the numerical solution is unstable and diverges. This phenomenon can be very difficult and annoying to solve numerically when studying system dynamics, especially when the system is large-scale or involves multiple disciplines.
This property of the ordinary differential equation is called stiffness (stiff) and this class of problem is called stiffness problem. The rigidity problem is also called a pathological equation or a bad condition equation in some documents, and the phenomenon is frequently encountered in practical engineering, for example, the inertia of a controlled object in a control system is large, a large time scale is provided, and a controller is sensitive in response and has a small time scale; for flight dynamics, its gesture motion is faster and the centroid motion is relatively slower; in a thermal fluid system, the pressure fluctuation is fast, the pressure is propagated in the fluid at sonic speed, the heat transfer process is slower, and the variable scale is obviously not on the same magnitude; for complex power systems, the time scale between different elements can also vary significantly, subject to parasitic capacitance.
Along with the third industrial revolution wave represented by high and new technologies such as atomic energy, aerospace, computers and the like at the beginning of fifties of the 20 th century, the numerical calculation of the system dynamics problem gradually goes from mathematical theory to engineering practice, and a plurality of open source programs and commercial simulation tools are further created.
The open source solver, which is relatively widely used, includes: the BDF algorithm is adopted and the version CVODE is still updated frequently so far; RADAU5 using the implicit Runge-Kutta method; dop853 and Dopri using the explicit finger-Kutta method; SDIRK adopting a diagonal implicit Runge-Kutta method; MEBDF based on a modified generalized BDF algorithm; DIFSUB and LSODE, which are designed and widely used based on the classical BDF algorithm; the open source computation package of Python, scipy.integration.ode (with calls LSODE and ODEPACK at the bottom); deSolve engineering packages based on an explicit range-Kutta method in the R language (according to the description of the official network, the bottom layer is still an open source ODEPACK computing package for calling the FORTRAN language, a package based on an explicit Euler/range-Kutta method and a package RADAU based on an implicit range-Kutta method); julia suites "differential equations. Jl" for differential equation high performance solver, and the like.
In addition to open source solvers, there are many very excellent commercial simulation solvers that have also gradually emerged in the general public after the seventies of the last century. Among all commercial system simulation tools, MATLAB/simulink is generally considered to have the best comprehensive performance, and an ODE solver of the MATLAB/simulink comprises almost all conventional algorithms which can be seen, such as BDF (binary coded field) method, modified BDF method, adams method, explicit or implicit Runge-Kutta method and the like; in addition, dymola, simulationX, motion, adams and AMESim are also very good dynamics solving tools, and because models processed by the tools are mainly multidisciplinary or multisome dynamics and other problems, the models are generally more rigid, so that solvers of the software basically adopt an implicit multi-step method, and BDF and a modified method thereof are the most common.
Mworks in China is also a simulation tool for Modelica model dynamics solving, but a solver is still realized by calling some foreign algorithm packages. Compared with the accumulation of the foreign prior art, china has obvious defects in the following two aspects: on the one hand, domestic ODE-IVP solvers which are autonomous and made in China are not available, and most of foreign program packages (including business and open sources) are called, such as CVODE, dassl, radau, MEBDF, LSODE and the like, and the domestic ODE is completely blank in the field; on the other hand, the foreign open source package formation years are early, and most of them do not include new methods in recent years. In addition, foreign business tools are widely used in the China market, the solver kernel is not greatly changed around 2000, and a new algorithm is not incorporated in recent years as in an open source package.
Disclosure of Invention
The invention solves the technical problem of providing a numerical integration solving method for the rigid ordinary differential equation initial value problem in simulation, and providing bottom support for solving system-level simulation (especially large-scale and multidisciplinary system simulation).
The invention is realized by the following technical scheme:
A numerical integral solving method of a rigid ordinary differential equation initial value problem in simulation comprises the following operations:
1) Aiming at the problem to be solved, selecting necessary parameters including an integral solving algorithm, an initial integral step length and a convergence residual error;
wherein the integral solution algorithm is selected from the BDF algorithm, IBDF algorithm, SDISRK algorithm, and SDISRK algorithm;
2) Judging according to the selected solving algorithm: if a multi-step method is adopted, judging that the current solver adopts a BDF or IBDF method, and continuing to execute the step 3 after judging;
if a single-step method is adopted, judging whether the current solver adopts SDISRK algorithm or SDISRK algorithm, and jumping to execute the step 5 after judging is completed;
3) Constructing a target nonlinear equation set according to the selected solving algorithm;
4) Respectively executing a plurality of prediction step algorithms to obtain prediction values of a plurality of real solutions, respectively substituting the prediction values into a target nonlinear equation set, and selecting a group of solution vectors with the smallest error from the prediction values as solution vectors of a prediction step;
5) Taking the solution vector obtained in the prediction step as an iteration initial value, and solving a least square solution of a target nonlinear equation set;
firstly constructing an iteration formula when solving, and introducing damping coefficients of iteration step sizes; updating the iteration step length and judging whether convergence requirements are met or not; if yes, outputting a result as a real solution vector, otherwise, continuously repeating the iterative process to obtain a new solution vector until the convergence error requirement is met;
6) Performing a variable step integration strategy: respectively executing a step length optimizing mechanism for providing an integral step length, a switching mechanism for judging and changing the integral step length and a variable step length numerical value back-inserting mechanism for ensuring the correct calculation result;
After the execution of the variable step integration strategy is finished, judging whether the error meets the requirement; if yes, executing 7);
if not, returning to the step 1), and resetting the integral step according to the step-changing strategy;
7) Advancing the time axis forward by one time step, and judging whether the simulation needs to be ended or not;
executing 8) if the result is satisfied; if not, returning to the step 1), and resetting the integral step according to the step-changing strategy;
8) And finishing the simulation calculation, and storing and displaying the calculation result.
Compared with the prior art, the invention has the following beneficial technical effects:
The numerical integration solving method of the rigid ordinary differential equation initial value problem in the simulation provides a solver suitable for solving the rigid and non-rigid ordinary differential equation initial value problem (ODE-IVP) for the standard foreign mainstream business and open source system simulation solver, and provides bottom support for solving the system-level simulation (especially large-scale and multidisciplinary system simulation).
The invention provides a damping Newton method based on classical Levenberg-Marquardt optimization algorithm, and provides a weighted-Levenberg-Marquardt optimization method (d-LM method for short); the algorithm can achieve both running efficiency and stability, and the method is equivalent to the foreign commercial tool in efficiency by comparison with similar calculation programs in the foreign commercial tool MATLAB.
In the prediction-correction method, the solution is predicted by a plurality of relatively simple explicit algorithms, so as to give a point closer to the true solution, ensure that the implicit solution in the next step is easier to converge, and consider the stability of the solution for protection, but slightly sacrifice the efficiency; although some calculation efficiency is lost in the prediction step link, the processing method can provide a better iteration initial value for the correction step, further reduce the iteration number m, and find that the processing is obviously improved on the global efficiency through comprehensive balance and a large number of tests.
The invention also adopts a new variable step length strategy capable of greatly improving the integral solving efficiency, and the operation efficiency is at the same level as that of the foreign similar products on the basis of filling the domestic blank; the d-LM method is further provided on the basis of guaranteeing the solving precision, and the operation stability is improved. Compared with foreign commercial system simulation software tools, the numerical method is more novel, a new numerical solution method after 1990 is realized and applied, and the numerical solution process is more stable.
Drawings
Figure 1 is a schematic overall flow diagram of the present invention,
FIGS. 2-5 are enlarged views of the first portion-fourth portion of FIG. 1, respectively;
FIG. 6 is a schematic diagram of a Butcher matrix of the 3-level 4-order single-diagonal implicit octark method;
FIG. 7 is a schematic diagram of a Butcher matrix of the 4-level 4-order single-diagonal implicit octyl RK method;
FIG. 8 is a schematic workflow diagram of variable step integration;
FIG. 9 is a schematic diagram of an organization architecture of the solver core code;
FIG. 10 is a schematic of the calculation result of the d-LM algorithm in the solver;
FIG. 11 is a graph showing the results of a calculation using fsolve commands in MATLAB;
FIG. 12 is a comparison of the calculated time consumption of the d-LM method and fsolve commands;
FIG. 13 is a result of solving a stiffness equation using simulink business software;
FIG. 14 is a comparison of the results of the present solver and simulink commercial software calculations;
FIG. 15 is a comparison of the computational efficiency of the present solver with a foreign business tool;
FIG. 16 is a numerical solution for y1 and y2 and the error of the solution for both;
FIG. 17 is a numerical solution for y3 and y4 and the error of the solution for both;
FIG. 18 is a schematic diagram of a Lorenz system three-dimensional evolution trajectory;
FIG. 19 is a schematic diagram of a pipeline model;
FIG. 20 is a schematic diagram of the pressures at the left and right ends of a pipeline of the pipeline model;
FIG. 21 is a graph showing the comparison of the present invention with a simulink calculation under pipeline model and pressure boundary conditions.
Detailed Description
The invention is described in further detail below in connection with examples, which are intended to be illustrative rather than limiting.
To implement an efficient, stable solver in simulation, three main concerns are:
Firstly, a proper calculation method is selected, so that the method has better stability, otherwise, the method is not easy to converge; meanwhile, the numerical solution precision is considered, and the selected numerical algorithm order cannot be too low;
Second, it is sufficient to consider whether the selected algorithm is suitable for typical stiffness problems, including but not limited to: rigid oscillation problems, nonlinear delay differentiation problems, rigid problems with discrete events, etc., which are typical of those that exist extensively in power system simulations;
Finally, in the process of realizing the algorithm, the running efficiency of the code is improved as much as possible, otherwise, the method has no engineering significance.
To this end, the present invention proposes the following methods and strategies to solve.
1. ODE solver developed based on linear multi-step method and Runge-Kutta method
Based on the theories of the linear multi-step method and the range-Kutta method, the solver selects four numerical algorithms, namely BDF, improved BDF and two 4-order single diagonal implicit octaves RK, which have good stability and can meet the requirements on accuracy, as a bottom core processing method. The numerical integral solver suitable for the rigid ordinary differential equation initial value problem can solve the ODE-IVP problem efficiently and stably, and particularly has a good effect on the rigid problem. The solver provides a variable step length mechanism for the four solving methods, and can fully exert the advantages of a large stable domain.
2. Novel method for solving nonlinear equation set least square problem
When solving the initial value problem of the ordinary differential equation by using the hidden format algorithm, the least square problem of the nonlinear equation set needs to be calculated. The solver fuses damping Newton method based on classical Levenberg-Marquardt optimization algorithm, and proposes a dammed-Levenberg-Marquardt optimization method, namely d-LM method. The algorithm can achieve both running efficiency and stability, and the method is equivalent to the foreign commercial tool in efficiency by comparison with similar calculation programs in the foreign commercial tool MATLAB.
3. Optimization strategy for predictive step algorithm
The linear multi-step method is a hidden-format solving method, and the most common method for solving an implicit equation is the "prediction-correction" method, which is also called the PC method class. The general step of the "predict-correct" method is P (EC) m, where P represents the prediction (Predictor), E represents the derivative value of the predicted step result, C represents the correction (Corrector), and m represents multiple iterations.
In order to ensure the convergence and the solving efficiency of m implicit iterations, a prediction result close to a true solution needs to be given as far as possible in a prediction step, and 5 calculation methods are provided by the solver, and 5 groups of prediction results are obtained by respectively applying the 5 calculation methods, wherein a group closest to the true solution is preferably selected as the input of the correction step. Although some calculation efficiency is lost in the prediction step link, the processing method can provide a better iteration initial value for the correction step, further reduce the iteration number m, and find that the processing is obviously improved on the global efficiency through comprehensive balance and a large number of tests.
4. Variable step length new strategy capable of greatly improving integral solving efficiency
Conventional variable step size solving strategies generally attempt to solve for the next value y (h) by taking h as an integral step size, and then reduce the step size (e.g., h/2 or h/3) according to a given ratio to solve for the next value y (h/2) or y (h/3) again. Then, judging whether a certain norm of the y (h)-y(h/2) or the y (h)-y(h/3) is smaller than a preset error epsilon, if so, considering that the current step size is proper, and advancing the time axis forward; if not, the current time step is not suitable, and a new integration step h new needs to be recalculated according to a certain step reduction strategy.
According to the invention, a judgment strategy is additionally added on the basis of the strategy, the concepts of step length adjustment and step length stabilization are introduced in the step length changing process, and when a plurality of integration step lengths h and errors epsilon are continuously kept in a range appointed in advance in the solving process, the current solved model is considered to enter a step length stable state, and the integration step length of the previous step is directly used without solving and adjusting the integration step length in the state; otherwise, the current solved model is considered to enter a step length adjustment state, and then the solution is carried out according to a traditional step length changing solving strategy. The processing method has the greatest advantages that on the premise of ensuring the accuracy of the calculation result, the calculation times are greatly reduced, and the solving efficiency is further improved.
The simulation solving method of the present invention is specifically described below.
The invention provides a numerical integral solving method of an initial value problem of a rigid ordinary differential equation in simulation, which aims to be equivalent to solving the initial value problem of an N-dimensional Chang Weifen equation of a real number domain, and has the expression as follows:
Where y is the derivative variable, y' represents the derivative of y, t is time, t 0 and t end represent the start and end integration times, N is the system of equations dimension (i.e., the dimension of y), and y 0 represents the value of y at time t 0.
Based on the linear multistep method and the theory of the Runge-Kutta method, four numerical algorithms with good stability and accuracy meeting the requirements are selected as the bottom core processing method:
1. classical BDF (3-5 th order) methods, also known as GEAR methods;
2. the BDF method is improved, namely IBDF (3-5 order) algorithm;
the 3.3-level 4-order single-diagonal implicit octaringe-Kutta method is SDISRK algorithm;
The 4.4-level 4-order single-diagonal implicit octaringe-Kutta method is SDISRK algorithm.
On the basis of the method, a set of numerical integral solver suitable for the initial value problem of the rigid ordinary differential equation is constructed, so that the ODE-IVP problem can be solved efficiently and stably, and the method has a good effect especially for the rigid problem; and a variable step length mechanism is provided for the four solving methods, so that the advantages of a large stable domain can be fully exerted.
Referring to fig. 1-5, a numerical integral solution method for an initial value problem of a rigid ordinary differential equation in simulation comprises the following operations:
Step 1: selecting necessary parameters such as an integral solving algorithm, an initial integral step length, a convergence residual error and the like, and continuously executing the step 2 after finishing;
wherein the integral solution algorithm is selected from the BDF algorithm, IBDF algorithm, SDISRK algorithm, and SDISRK algorithm;
step 2: judging according to the selected solving algorithm:
If a multi-step method is adopted, judging that the current solver adopts a BDF or IBDF method (an improved BDF method), and continuing to execute the step 3 after judging;
if a single-step method is adopted, the current solver is judged to adopt SDISRK or SDISRK methods, and the step 5 is carried out after the judgment is completed;
Step 3: constructing a target nonlinear equation set according to the selection of the solving algorithm in the previous step, and continuously executing the step 4 after finishing;
step 4: executing a preferred strategy of a prediction step result, respectively executing a plurality of prediction step algorithms to obtain a plurality of prediction values of real solutions, respectively substituting the prediction values into a target nonlinear equation set, and selecting a group of solution vectors with the smallest error as the solution vectors of the prediction step; after finishing, continuing to execute the step 5;
Step 5: taking the solution vector obtained in the prediction step as an iteration initial value, and solving a least square solution of a target nonlinear equation set by using a d-LM method;
firstly constructing an iteration formula when solving, and introducing damping coefficients of iteration step sizes; updating the iteration step length and judging whether convergence requirements are met or not; if yes, outputting a result as a real solution vector, otherwise, continuously repeating the iterative process to obtain a new solution vector until the convergence error requirement is met;
step 6: performing a variable step integration strategy: respectively executing a step length optimizing mechanism for providing an integral step length, a switching mechanism for judging and changing the integral step length and a variable step length numerical value back-inserting mechanism for ensuring the correct calculation result;
judging whether the error meets the requirement; if yes, executing step 7;
if not, returning to the step 1, and resetting the integral step according to the step-changing strategy.
Step 7: advancing the time axis forward by one time step, and judging whether the simulation needs to be ended or not;
if yes, executing step 8;
if not, returning to the step 1, and resetting the integral step according to the step-changing strategy.
Step 8: and finishing the simulation calculation, and storing and displaying the calculation result.
A more specific step description is given below.
A. In the step 2 of the solution, the solver provides four solving algorithms, namely two linear multi-step methods and two single-step high-precision methods, which are respectively named as a method A to a method D.
Method a (3-5 th order conventional BDF):
Method B (3-5 order modified BDF):
method C (3-level 4-order single diagonal implicit octange-Kutta):
The coefficients of the 3-level 4-order single-diagonal implicit oct method are generally given in the form of a Butcher matrix, as shown in fig. 6.
Method D (4-level 4-order single diagonal implicit octange-Kutta):
the 4-level 4-order single-diagonal implicit octreotide method expression differs from method C only in terms of the coefficients of the Butcher matrix. The Butcher matrix of the 4-level 4-order single-diagonal implicit octark method is shown in fig. 7.
The general equation is obtained by using a method A or a method B, the calculation amount of the solving process of the two methods is small, and the solving efficiency is high.
For some problems such as rigid oscillation, random, discrete-continuous mixing or complex nonlinearity, an attempt can be made to use a method C or a method D, the calculation amount of the solving process of the two methods is large, the solving efficiency is relatively low, but the advantage is that the numerical stability is good (the method C and the method D are A-stable, the method A or the method B is A (alpha) -stable), and meanwhile, the two methods belong to a pungent format (SYMPLECTIC SCHEMES), and have smaller numerical dissipation compared with the method A and the method B, so that the numerical fidelity is better.
B. construction of a set of target nonlinear equations
The four solving algorithms provided by the invention are all implicit solving, so that a target nonlinear equation set is required to be constructed.
The target nonlinear equation set expression of method a is:
the target nonlinear equation set expression of method B is:
In the nonlinear equation set, y n+k is an unknown number to be solved, the rest are known data, and F (y n+k) is a target nonlinear equation set.
The target nonlinear equation system of the method C and the method D solves the coefficients of the hidden RK, and the expression is as follows:
In the nonlinear equation set, K i is an unknown number to be solved, the rest are known data, and F (K i) is a target nonlinear equation set.
C. preferred strategy for predicting step outcome
Prediction step calculation is used in method a and method B. In the prediction step, the solver provides 5 algorithms for making predictions of the true solution:
1. The explicit Euler method is used as a prediction step to obtain a predicted value x 1 of a real solution;
2. displaying a 2-order Runge-Kutta method as a prediction step to obtain a predicted value x 2 of a real solution;
3. Displaying a 3-order Runge-Kutta method as a prediction step to obtain a predicted value x 3 of a real solution;
4. displaying a 4-order Runge-Kutta method as a prediction step to obtain a predicted value x 4 of a real solution;
5. And using the result of the previous step as a prediction step to obtain a predicted value x 5 of the true solution.
The above 5 solution vectors are respectively substituted into a 'target nonlinear equation set', and a set of solution vectors x pre with the smallest error is selected from the solution vectors as an input to a 'correction step'.
D. solving a target nonlinear equation set by using d-LM method
And taking the solution vector x pre obtained in the prediction step as an iteration initial value, and solving a least square solution of a 'target nonlinear equation set' by using a d-LM method, wherein a final result is a true solution vector.
The solving process of the d-LM method is as follows:
1) Introducing a target nonlinear equation set, and constructing an iterative formula based on a Levenberg-Marquardt method
(JTJ+ηI)hLM＝-JTF,g＝JTF,η≥0 (10)
Wherein F is a target nonlinear equation set to be solved, J is a Jacobian matrix of the target nonlinear equation set F, J T is a transpose matrix of J, eta is an adjustment coefficient between a gradient descent method (first-order optimization algorithm) and a Newton iteration method (second-order optimization algorithm), I is an identity matrix, rho is a trust domain radius, h LM is a convergence step length of each iteration step, and v is a process variable for calculating the trust domain radius.
2) Calculating damping coefficient
The most prominent advantage of newton's method is the fast convergence speed with local second order convergence, but the initial point of basic newton's method needs to be sufficiently "close" to the minimum point, otherwise it may cause the algorithm to not converge. In order to ensure the stability of the calculation process, an iteration step damping coefficient is introduced by using a damping Newton method.
Wherein lambda is the iteration step damping coefficient to be solved, h LM is the convergence step of each iteration step, and x is the solution vector at the current moment.
3) Updating the iteration step length, judging whether the convergence requirement is met, if so, outputting a result, otherwise, continuously repeating the iteration process to continuously obtain a new solution vector x new until the convergence error requirement is met.
hLM＝(JTJ+ηI)-1(-JTF) (15)
xnew＝x+λhLM (16)
E. Variable step integration strategy (see FIG. 8)
The four solution methods A-D provided by the invention allow the use of a variable step integration strategy, and the strategy can obviously improve the integration solution efficiency; the variable step integration strategy mainly has the following three mechanisms:
the first is an integral step length optimizing mechanism; in this mechanism, the step size is decremented to an exponentially decreasing relationship:
Wherein dt new is the updated integral step length, dt is the current integral step length, error is the current calculation error, epsilon is the preset allowable error, y p is the solution vector obtained by solving the p-step method, and z p+1 is the solution vector obtained by solving the p+1-step method;
Y p and z p+1 are calculated respectively, when error is less than or equal to 1.0, the calculation result meets the given allowable error epsilon requirement, the current calculation result can be accepted, and the current integration step is finished;
the second is a judging and switching mechanism between two states of step length adjustment and step length stabilization
Each integration step should theoretically be given an initial integration step and then iterated through equations (17) and (18) to get an appropriate integration step. Such a compensation adjustment mechanism is the "step adjustment" state. If each integral step is calculated according to the mechanism, the calculation result is necessarily correct, however, the solving efficiency is necessarily very slow, and a step stable state is also provided for further improving the solving efficiency.
When a plurality of continuous (5 continuous) integral step sizes are calculated according to the step size adjustment state, the next time step size can be calculated directly by adopting the average value of the plurality of previous integral step sizes without iteration, but the monitoring is needed by using a formula (18), so that the condition that 0.75 error is less than or equal to 1 is met, the measurement directly considers that the current integral step is calculated, and the time axis can be directly advanced.
The third is a variable step value back-insertion mechanism that is activated only in a multi-step process, and that is not required to be performed in a single step process. When the variable step length solution is adopted, as each integral step length is different, the two linear multistep methods (methods A and B) provided by the solver require that the integral step lengths of the first steps are the same, otherwise, numerical divergence is caused, so that a numerical value back-interpolation method is introduced, wherein Interp _1D is a one-dimensional interpolation function:
The three mechanisms have the following beneficial effects on the solver:
the step length optimizing mechanism provides proper integral step length for the solver, so that the result is not distorted due to overlarge calculation times and the whole solving efficiency is not dragged due to overlarge calculation times;
the judging and switching mechanism of the two states of step adjustment and step stabilization can reasonably judge whether the integral step is required to be changed or not, if not, the integral step is not adjusted as much as possible, the extra calculation amount brought by step adjustment is saved, and the solving efficiency is further improved;
the variable step size numerical value back-insertion mechanism ensures the accuracy of calculation results in a linear multi-step method and avoids result distortion caused by step size change.
The three mechanisms support the variable step integration strategy together, which is indispensable.
Specifically, the solver of the invention adopts C++ as a development language, relies on an open source linear algebraic computation package Eigen (V3.3.9), and adopts an Intel mathematical core function Library (MKL, math Kernel Library) to optimize and accelerate, and the basic code organization architecture is shown in FIG. 9.
Specific examples are given below.
Example 1: efficiency comparison of d-LM method and fsolve commands
BDF and IBDF in the solver are two very efficient solving methods, solving a nonlinear equation set is an indispensable important link, and the solver innovatively provides a d-LM method; in order to verify the solving efficiency, MATLAB is selected as a comparison verification reference, and the difference between the d-LM method based on C++ development and fsolve commands for solving the nonlinear equation set in MATLAB, which are given by the solver, is compared.
The present example selects the following nonlinear equation set as the test object.
The initial value of the solution vector of the equation is [0, 0], d-LM and fsolve of MATLAB of the solver are adopted to solve, and the calculation results are shown in FIG. 10 and FIG. 11 respectively.
From the above calculations, it can be seen that the same solution is obtained for the same mathematical model, i.e. x= [10.0449,4.9796, -0.0653], and the correctness is verified.
In addition, it can be seen that fsolve steps in MATLAB need to iterate 9 steps to enable the equation set to achieve 1e-12 solving accuracy, and the d-LM method provided by the solver can achieve 1e-12 solving accuracy only by 7 steps, so that fewer iteration times mean higher efficiency possibility.
FIG. 12 shows the solution time-consuming statistics of the d-LM method and fsolve of the present invention when repeated 20 solution calculations. As can be seen from the calculation results, the d-LM method and fsolve are at the same level, but the fsolve command is inefficient at the beginning of calculation, and the efficiency is improved after repeated for a plurality of times.
Example 2: example of typical stiffness problem
The following types of reaction velocity equations are often encountered in chemical reactions, where y 1,y2,y3 represents the concentration of the reactants, and this model is a typical stiffness problem, often used for simulation verification of various numerical solutions.
The results shown in FIG. 13 can be obtained by solving the ordinary differential equation above using MATLAB/simulink business tool. Taking the result as a comparison standard, solving the ordinary differential equation by using a 4-step improved BDF algorithm in the solver, and obtaining an error curve shown in figure 14 by making a difference between the obtained result and the simulink result, so that the consistency of the two results is very good.
The variable step convergence error in the current example is 1e-5, so the error of the resulting solution is also substantially around this order. If a more accurate solution is desired, the convergence error of the variable step is further tightened, but this certainly affects the calculation efficiency, which is a pair of contradictions in the solution process and needs to be comprehensively considered.
In addition to comparing accuracy, the present invention also compares the efficiency differences of different solutions. FIG. 15 shows a comparison of the solution efficiency of a dominant business tool abroad compared to the time consumed by the present solver.
As can be seen, the computation time of Simulink and SimulationX is distributed between 0.4 and 2.0s, and the BDF method of the solver is mainly distributed between 0.7 and 1.0s. Overall, the present solver is substantially comparable in computational efficiency to the mainstream business tool abroad. The SDISRK method is an implicit Runge-Kutta method, which is more time-consuming to calculate due to the larger calculation amount of the algorithm.
Example 3: typical Hamilton System oscillation problem
The Hamilton system is a typical kinetic system and is widely used in describing the system dynamics in the fields of celestial motion, mechanical motion, etc.
The Kepler equation is a typical Hamilton system, expressed as:
the Hamilton function is:
when the integral initial value is y 1(0)＝0,y2(0)＝1,y3(0)＝1,y4 (0) =0, the mathematical expression of the analytical solution is:
y1＝-sin(t),y2＝cos(t),y3＝cos(t),y4＝sin(t)
The oscillation equations described by formulas (22) and (23) are subjected to numerical integration solution by using the solver, and the obtained results are shown in fig. 16 and 17.
It can be seen that the numerical solution and the theoretical analytical solution have good consistency, and the errors of the two are of the order of 1e-5, which is related to the preset of the solver, but more accurate numerical solutions require more calculation time.
Example 4: lorenz chaotic system
The Lorenz system is a typical chaotic system, the dynamics process of the Lorenz system is complex, and the Lorenz system is a typical means for checking the correctness of a differential equation solving algorithm. The expression of the differential equation is:
The solver is utilized to carry out numerical solution on the Lorenz chaotic system, and the result is drawn into a three-dimensional evolution track which is compared with the MATLAB/simulink calculation result, see FIG. 18, so that the chaotic attractors reflected by the calculation results of the two are completely consistent, and the correctness of the solver can be proved.
Example 5: R-C-R one-dimensional fluid pipeline system
In the multidisciplinary system simulation, the flow resistance (R) and the flow capacity (C) are mainly considered in the pipeline model, the flow resistance mainly represents the resistance of the pipeline, and the flow capacity mainly represents the compressibility of the pipeline. When the pipeline is longer, the pipeline can be processed in a segmented way, and the corresponding model can be represented by R-C-R- & gtC-R.
For a simple circular pipeline with uniform cross section, the circular pipeline can be simply equivalent to an R-C-R model (shown in fig. 19), and the kinetic model is as follows:
wherein P in,out is the inlet and outlet pressure of the pipeline, The flow resistance is R, ρ, η, A, V, A and a speed of sound of fluid medium. In calculation, medium density ρ=1000 kg/m 3, pipeline length 12.73m, pipeline inner diameter 10mm and medium sound velocity a=1400 m/s are taken.
The calculated pressures at the left and right ends of the pipeline are shown in fig. 20, and under the same pipeline model and pressure boundary conditions, the calculation results of the solver and the simulink are calculated and compared, and are shown in fig. 21. It can be seen that the calculation result of the solver is completely consistent with the calculation result of simulink.
The embodiments given above are preferred examples for realizing the present invention, and the present invention is not limited to the above-described embodiments. Any immaterial additions and substitutions made by those skilled in the art according to the technical features of the technical scheme of the invention are all within the protection scope of the invention.
Claims (3)
1. The numerical integral solving method of the rigid ordinary differential equation initial value problem in the simulation is characterized by comprising the following operations:
1) Aiming at the simulation problem to be solved, including a chemical reaction speed equation and a fluid pipeline model, selecting necessary parameters including an integral solving algorithm, an initial integral step length and a convergence residual error;
wherein the integral solution algorithm is selected from the BDF algorithm, IBDF algorithm, SDISRK algorithm, and SDISRK algorithm;
2) Judging according to the selected solving algorithm: if a multi-step method is adopted, judging that the current solver adopts a BDF or IBDF method, and continuing to execute the step 3 after judging;
if a single-step method is adopted, judging whether the current solver adopts SDISRK algorithm or SDISRK algorithm, and jumping to execute the step 5 after judging is completed;
3) Constructing a target nonlinear equation set according to the selected solving algorithm;
4) Respectively executing a plurality of prediction step algorithms to obtain prediction values of a plurality of real solutions, respectively substituting the prediction values into a target nonlinear equation set, and selecting a group of solution vectors with the smallest error from the prediction values as solution vectors of a prediction step;
5) Taking the solution vector obtained in the prediction step as an iteration initial value, and solving a least square solution of a target nonlinear equation set;
firstly constructing an iteration formula when solving, and introducing damping coefficients of iteration step sizes; updating the iteration step length and judging whether convergence requirements are met or not; if yes, outputting a result as a real solution vector, otherwise, continuously repeating the iterative process to obtain a new solution vector until the convergence error requirement is met;
6) Performing a variable step integration strategy: respectively executing a step length optimizing mechanism for providing an integral step length, a switching mechanism for judging and changing the integral step length and a variable step length numerical value back-inserting mechanism for ensuring the correct calculation result;
After the execution of the variable step integration strategy is finished, judging whether the error meets the requirement; if yes, executing 7);
if not, returning to the step 1), and resetting the integral step according to the step-changing strategy;
7) Advancing the time axis forward by one time step, and judging whether the simulation needs to be ended or not;
executing 8) if the result is satisfied; if not, returning to the step 1), and resetting the integral step according to the step-changing strategy;
8) Ending the simulation calculation, and storing and displaying the calculation result;
the construction of the target nonlinear equation set is as follows:
the target nonlinear equation set expression of the BDF algorithm is:
The target nonlinear equation set expression of IBDF algorithm is:
In the nonlinear equation set, y n+k is an unknown number to be solved, the rest is known data, and F (y n+k) is a target nonlinear equation set;
The target nonlinear equation sets of SDISRK algorithm and SDISRK algorithm solve for coefficients in the implicit RK, whose expression is:
Wherein K i is an unknown number to be solved, the rest are known data, and F (K i) is a target nonlinear equation set;
the least squares solution of the set of target nonlinear equations is solved using the following operations:
501 Introducing a target nonlinear equation set, and constructing an iteration formula based on a Levenberg-Marquardt method:
(JTJ+ηI)hLM＝-JTF,g＝JTF,η≥0
if ρ > 0
Η=η×max {1/3,1- (2ρ -1) 3 }; v=2
Otherwise
η＝η*ν；ν＝2*ν
Wherein F is a target nonlinear equation set to be solved, J is a Jacobian matrix of the target nonlinear equation set F, J T is a transpose matrix of J, eta is an adjustment coefficient between a gradient descent method and a Newton iteration method, I is a unit matrix, rho is a trust domain radius, h LM is a convergence step length of each iteration step, and v is a process variable for calculating the trust domain radius;
502 Introducing an iteration step damping coefficient:
Wherein lambda is the iteration step damping coefficient to be solved, h LM is the convergence step length of each iteration step, and x is the solution vector at the current moment;
503 Updating the iteration step length, repeating the iteration process to continuously obtain a new solution vector x new until the convergence error requirement is met, and outputting a result:
hLM＝(JTJ+ηI)-1(-JTF)
xnew＝x+λhLM；
the execution of each mechanism in the variable step integration strategy is as follows:
601 A step size is reduced to an exponential decreasing relationship) optimizing mechanism:
Wherein dt new is the updated integral step length, dt is the current integral step length, error is the current calculation error, epsilon is the preset allowable error, y p is the solution vector obtained by solving the p-step method, and z p+1 is the solution vector obtained by solving the p+1-step method;
If error is less than or equal to 1.0, the calculation result meets the requirement of a given allowable error epsilon, the current calculation result can be accepted, and the current integration step is finished;
If error >1.0, then the current calculation result does not meet the given allowable error requirement, and a new dt new should be recalculated;
602 Judging the switching mechanism for changing the integral step length, which is a mechanism for judging and switching between two states of step length adjustment and step length stabilization:
The step length adjusting state is as follows: giving an initial integration step length on each integration step length, and then obtaining an integration step length through iteration of the formula (17) and the formula (18);
the step length steady state is as follows: when a plurality of continuous integral step sizes are all calculated according to the step size adjusting state, the next time step size is directly calculated by adopting the average value of the plurality of previous integral step sizes without iteration, and monitoring is carried out by using the step (18); if the condition that the error is less than or equal to 0.75 and less than or equal to 1 is met in the monitoring process, the current integration step is directly considered to finish calculation, and the time axis is directly pushed forward;
603 Variable step size numerical value back-insertion mechanism:
the mechanism is activated and validated only in a multi-step method, and a numerical value back-interpolation method is introduced, wherein Interp _1D is a one-dimensional interpolation function:
2. The method for solving the numerical integral of the rigid ordinary differential equation initial value problem in the simulation according to claim 1, wherein the rigid oscillation, random, discrete-continuous mixture or complex nonlinear problem is adopted by SDISRK algorithm or SDISRK algorithm when the integral solving algorithm is selected; the BDF algorithm or IBDF algorithm is used for the remaining equations.
3. The method for solving the numerical integral of the rigid ordinary differential equation initial value problem in simulation according to claim 1, wherein the prediction step algorithm is used for performing the prediction of the true solution by:
401 An explicit Euler method is used as a prediction step to obtain a predicted value x 1 of a real solution;
402 Displaying a 2-order Runge-Kutta method as a prediction step to obtain a predicted value x 2 of a real solution;
403 Displaying 3-order Runge-Kutta method as a prediction step to obtain a predicted value x 3 of a true solution;
404 Displaying 4-order Runge-Kutta method as a prediction step to obtain a predicted value x 4 of a true solution;
405 Using the result of the previous step as a prediction step to obtain a predicted value x 5 of a real solution;
Substituting the above 5 solution vectors into a target nonlinear equation set respectively, and selecting a group of solution vectors x pre with the smallest error from the solution vectors as input for correction steps:
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111479747.6A CN114154110B (en) | 2021-12-06 | 2021-12-06 | A numerical integration method for solving initial value problems of stiff ordinary differential equations in simulation |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111479747.6A CN114154110B (en) | 2021-12-06 | 2021-12-06 | A numerical integration method for solving initial value problems of stiff ordinary differential equations in simulation |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN114154110A CN114154110A (en) | 2022-03-08 |
| CN114154110B true CN114154110B (en) | 2024-11-26 |
Family
ID=80453136
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202111479747.6A Active CN114154110B (en) | 2021-12-06 | 2021-12-06 | A numerical integration method for solving initial value problems of stiff ordinary differential equations in simulation |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN114154110B (en) |
Families Citing this family (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2025001243A1 (en) * | 2023-06-25 | 2025-01-02 | 华为云计算技术有限公司 | Method for solving linear equations and related device |
| CN119090900A (en) * | 2024-08-19 | 2024-12-06 | 四川大学 | Skin cancer lesion segmentation method based on discrete neural memory ordinary differential equation |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103646152A (en) * | 2013-12-23 | 2014-03-19 | 南方电网科学研究院有限责任公司 | An Electromagnetic Transient Simulation Method of Power System Based on Matrix Exponent |
| CN105404610A (en) * | 2015-10-23 | 2016-03-16 | 三峡大学 | Electromagnetic transient calculation method based on two-stage three-order single diagonally implicit Runge-Kutta method |
Family Cites Families (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CA2807300C (en) * | 2010-09-20 | 2017-01-03 | Exxonmobil Upstream Research Company | Flexible and adaptive formulations for complex reservoir simulations |
| CN107290977B (en) * | 2017-06-07 | 2019-09-27 | 清华大学 | Backward discrete state event-driven simulation of power electronic method, equipment and medium |
| CN108416132B (en) * | 2018-02-28 | 2021-11-12 | 东南大学 | Automatic variable-step-size simulation acceleration method for distributed photovoltaic cluster |
| CN109211276B (en) * | 2018-10-30 | 2022-06-03 | 东南大学 | SINS initial alignment method based on GPR and improved SRCKF |
-
2021
- 2021-12-06 CN CN202111479747.6A patent/CN114154110B/en active Active
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103646152A (en) * | 2013-12-23 | 2014-03-19 | 南方电网科学研究院有限责任公司 | An Electromagnetic Transient Simulation Method of Power System Based on Matrix Exponent |
| CN105404610A (en) * | 2015-10-23 | 2016-03-16 | 三峡大学 | Electromagnetic transient calculation method based on two-stage three-order single diagonally implicit Runge-Kutta method |
Also Published As
| Publication number | Publication date |
|---|---|
| CN114154110A (en) | 2022-03-08 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Park et al. | Unstructured grid adaptation: status, potential impacts, and recommended investments towards CFD 2030 | |
| CN109190233B (en) | A structure topology optimization method | |
| CN107016155B (en) | Convergence estimation for nonlinear PDE and linear solver | |
| CN114154110B (en) | A numerical integration method for solving initial value problems of stiff ordinary differential equations in simulation | |
| Gu et al. | Co-simulation of algebraically coupled dynamic subsystems without disclosure of proprietary subsystem models | |
| CN117436370A (en) | Overdetermined matrix equation parallel method and system for fluid mechanics grid generation | |
| CN118504151B (en) | Structural mode calculation method based on embedded physical information graph neural network | |
| US20050256683A1 (en) | Method and arrangement for designing a technical system | |
| Yang et al. | A new active learning method for reliability analysis based on local optimization and adaptive parallelization strategy | |
| Yu et al. | On the efficiency of the advancing-front surface mesh generation algorithm | |
| JP2011154439A (en) | Optimization processing program, method, and apparatus | |
| Ferreau et al. | A parallel active-set strategy to solve sparse parametric quadratic programs arising in MPC | |
| Morales et al. | A finite element method for active vibration control of uncertain structures | |
| CN120600175A (en) | A topology optimization method based on physical information neural network (PINN) | |
| Luna-Ortiz et al. | An input/output model reduction-based optimization scheme for large-scale systems | |
| CN118378572B (en) | Intelligent optimization design method and device for liquid rocket engine flow regulator | |
| CN117634298A (en) | Neural network model generalization enhancement method for predicting compressible two-gas Riemannian solutions | |
| Senecal | Efficient coupling algorithms and reduced-order methods for high-fidelity multiphysics simulations of nuclear reactors | |
| Ugolotti et al. | Differentiated ML-based modeling of structured grids for gradient-based optimization | |
| Alvarruiz et al. | A toolkit for water distribution systems’ simulation using the loop method and high performance computing | |
| Persson et al. | How to compare performance of robust design optimization algorithms, including a novel method | |
| Clarich et al. | Formulations for Robust Design and Inverse Robust Design | |
| CN119601145B (en) | A gas-solid coupling vibration simulation method for ceramic matrix composites | |
| Wang et al. | Towards the efficient calculation of quantity of interest from steady Euler equations I: a dual-consistent DWR-based h-adaptive Newton-GMG solver | |
| CN117435308B (en) | A Modelica model simulation method and system based on parallel computing algorithm |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |