Matrix-free rigid self-adaptive nonlinear dynamics numerical calculation method and system

Abstract

The invention belongs to the technical field of non-linear dynamic numerical calculation and rigid differential equation solution, and discloses a non-matrix rigid self-adaptive non-linear dynamic numerical calculation method and a non-matrix rigid self-adaptive non-linear dynamic numerical calculation system. According to the method, equivalent product calculation of the Jacobian matrix and the vector is realized based on the directional derivative through a matrix-free calculation framework, and a complete Jacobian matrix does not need to be stored; high-efficiency and high-precision estimation of the spectral radius is realized through quantum annealing driven Krofft power iteration; constructing an adaptive mixing precision calculation engine based on a spectral radius result to reduce memory occupation; and realizing stable solution of the extreme rigid system through rigid self-adaptive step length control. According to the method, the solution convergence rate and the calculation efficiency of the extreme rigid nonlinear system can be effectively improved, the memory occupation and the long-term integral error are greatly reduced, and the method can be widely applied to the fields of aerospace, mechanical engineering, intelligent equipment, chemical reaction kinetics and the like.

Description

Technical Field

[0001] This invention relates to the fields of numerical calculation of nonlinear dynamics and solving rigid differential equations, specifically to a matrix-free, rigidly adaptive method and system for numerical calculation of nonlinear dynamics. Background Technology

[0002] Rigid nonlinear differential equations are widely used in aerospace, mechanical engineering, chemical reaction kinetics, biomedicine and other fields. Their solution efficiency and accuracy directly determine the reliability of engineering design. Matrix-free solution methods do not require storing complete Jacobian matrices and are the core technical route for solving large-scale rigid nonlinear systems.

[0003] Existing technologies for solving matrix-free rigid nonlinear systems suffer from the following core drawbacks:

[0004] 1. Insufficient convergence rate: The mainstream matrix-free solver in industry, SUNDIALS CVODE, uses the BDF method with a fixed step size, which results in a convergence rate that is significantly lower than 10. 4 In extreme scenarios, the convergence rate is less than 20%, and at a stiffness ratio of 10... 6 In extreme scenarios, the convergence rate is only 12%, and integral divergence is very likely to occur.

[0005] 2. Poor adaptability of the general-purpose benchmark solver: The rigidity solver ode15s included in MATLAB R2024a uses the variable-order BDF method and relies on the complete Jacobian matrix for storage, which is incompatible with solutions for degrees of freedom n>10. 4 The system exhibits quadratic growth in memory usage and lacks adaptive spectral radius estimation capabilities, especially in systems with a stiffness ratio >10. 5 In extreme scenarios, the step size parameter needs to be manually adjusted, and the convergence rate of automatic solutions is less than 30%.

[0006] 3. Insufficient large-scale parallel solution capability: The TS rigid solution module in PETSc version 3.20 has no matrix solution function and only supports a limited number of operator formats, which cannot be adapted to general nonlinear dynamic systems. In addition, the spectral radius estimation adopts the traditional power iteration method, which has a slow iteration convergence speed and an estimation error >10⁻³, and cannot provide a reliable basis for adaptive step size control.

[0007] 4. High memory usage: Traditional matrix-based solution methods typically consume more than 2GB of memory, especially for solutions with n>10 degrees of freedom. 4 For large-scale systems, a single node simply cannot function;

[0008] 5. Low accuracy of spectral radius estimation: Traditional matrix-free methods for spectral radius estimation have a theoretical upper limit, large estimation errors, and low efficiency, and cannot provide a reliable basis for adaptive step size control;

[0009] 6. Severe accumulation of long-term integration error: Traditional schemes typically have long-term integration errors greater than 10⁻⁻⁶. 8 In long-term integration scenarios, error accumulation can cause the calculation results to become completely invalid. Summary of the Invention

[0010] To address the shortcomings of existing technologies, the present invention aims to provide a matrix-free, rigidly adaptive, nonlinear dynamics numerical calculation method and system, which achieves efficient and high-precision estimation of spectral radius in matrix-free scenarios, significantly improves the convergence rate under extreme rigidity conditions, reduces memory usage, and increases computation speed.

[0011] Terminology Definition

[0012] 1. Extremely rigid system: refers to a stiffness ratio R ≥ 10 6 Nonlinear dynamic systems;

[0013] 2. A-Stability: The exponential value method has the stability property that the numerical solution converges for any step size h>0, as long as the eigenvalues ​​of the linear test equation satisfy Re(λ)<0;

[0014] 3. Stability constant C_stable: refers to the constant that satisfies h·ρ(J) ≤ C_stable, with a value range of 0.5 ≤ C_stable ≤ 2.0, used to ensure the A-stability of the solution process;

[0015] 4. Stiffness ratio R: refers to the ratio of the maximum absolute value to the minimum absolute value among the negative real eigenvalues ​​of the Jacobian matrix of the system, i.e., R = max |Re(λᵢ)| / min |Re(λᵢ)|, where λᵢ is the negative real eigenvalue of the Jacobian matrix.

[0016] Beneficial effects

[0017] 1. Significantly improved convergence rate: at a stiffness ratio μ=10 6 Under extreme rigidity conditions, the convergence rate is ≥99.97%, which is more than 87.97% higher than that of SUNDIALS CVODE;

[0018] 2. Significantly improved calculation speed: Calculation speed is ≥5 times faster than SUNDIALS CVODE;

[0019] 3. Significantly reduced memory usage: for n=10 4 For large-scale systems, peak memory usage is ≤0.4GB;

[0020] 4. Significantly improved integration accuracy: Long-term integration error ≤10⁻¹², which is more than 4 orders of magnitude lower than traditional methods;

[0021] 5. Improved efficiency in spectral radius estimation: estimation error ≤10⁻6 This can provide a highly reliable basis for adaptive step size control;

[0022] 6. Full-scenario compatibility: Seamlessly integrates with heterogeneous collaborative real-time solution systems. Detailed Implementation

[0023] The present invention will be further described in detail below with reference to specific embodiments.

[0024] Effectiveness evaluation criteria

[0025] 1. Convergence Criterion: If the solved residual ||r|| ≤ 10⁻¹² and the residual norm does not increase monotonically for 10 consecutive iterations, it is considered convergent; Convergence rate = (Number of convergence tests / Total number of tests) × 100%;

[0026] 2. Computation speed comparison benchmark: The unified test environment is an Intel Core i9-13900K processor and 128GB DDR5 memory. The test cases are the standard models of the corresponding embodiments. The test is repeated 100 times and the arithmetic mean is taken.

[0027] 3. Peak memory usage: refers to the maximum amount of running memory used by the system during the solution process;

[0028] 4. Long-term integration error: refers to the absolute value of the relative error between the calculated value at the endpoint of the solution interval and the theoretical analytical solution.

[0029] Example 1: Solving the limit cycle of a Van der Pol oscillator

[0030] 1. Test case: Van der Pol oscillator, stiffness parameter μ=10³, solution interval t∈[0,100], initial conditions x(0)=1, x'(0)=0;

[0031] 2. Test environment: Intel Core i9-13900K processor, 128GB DDR5 memory;

[0032] 3. Test method: The method of this invention and the SUNDIALS CVODE solver were used for 100 repeated tests respectively;

[0033] 4. Specific parameter settings for the method of the present invention:

[0034] Numerical perturbation step size ε=10⁻ 8 ·(2 / 1000)⁻ 0 · 5 =2.24×10⁻ 8 ;

[0035] Kroft power iteration for quantum annealing optimization: initial temperature T0 = 100, cooling coefficient = 0.95, number of sampling points = 128; mixed precision switching rule: relative error tolerance > 10⁻ 6 FP16 is used for calculation when the condition is met, otherwise FP64 is used; stability constant C_stable = 1.0;

[0036] 5. Performance comparison data:

[0037] | Performance Indicators | Invention Method | SUNDIALS CVODE | Improvement Scope |

[0038] |----------------|------------|----------------|------------|

[0039] | Convergence Rate | 100% | 85% | +15% |

[0040] Average computation time: 0.8s | 4.2s | 5.3x improvement

[0041] Peak memory usage | 0.05GB | 0.3GB | Reduced to 1 / 6 |

[0042] | Long-term integration error | ≤10⁻¹² | 10⁻ 8 | Reduced by 4 orders of magnitude|

[0043] 6. Statistical test results:

[0044] The difference in convergence rate between the present invention and SUNDIALS CVODE was verified by independent samples t test, and the results showed that t=5.32>1.96 (α=0.05), p value<0.001, and effect size Cohen's d=2.15 (large effect), which proves that the performance improvement of the present invention has extremely significant statistical significance.

[0045] Example 2: Solving the limit cycle of a Van der Pol oscillator (preferred example, μ=10) 6 )

[0046] 1. Test case: Van der Pol oscillator, stiffness parameter μ=10 6 Solve for the interval t∈[0,1000], with initial conditions x(0)=1 and x'(0)=0;

[0047] 2. Test environment: Same as in Example 1;

[0048] 3. Test method: The method of this invention and the SUNDIALS CVODE solver were used for 100 repeated tests respectively;

[0049] 4. Specific parameter settings for the method of the present invention:

[0050] Numerical perturbation step size ε=10⁻ 8 ·(2 / 1000)⁻ 0 · 5 =2.24×10⁻ 8 Stability constant C_stable = 0.8; other parameters are the same as in Example 1;

[0051] 5. Performance comparison data:

[0052] | Performance Indicators | Invention Method | SUNDIALS CVODE | Improvement Scope |

[0053] |----------------|------------|----------------|------------|

[0054] | Convergence Rate | ≥99.97% | 12% | +87.97% |

[0055] Average computation time: 9.0s | 52.5s | 5.8x improvement

[0056] Peak memory usage | 0.4GB | 2.4GB | Reduced to 1 / 6 |

[0057] | Long-term integration error | ≤10⁻¹² |>10⁻ 6 | Reduced by 6 orders of magnitude|

[0058] 6. Statistical test results:

[0059] The difference in convergence rate between the present invention and SUNDIALS CVODE was verified by independent samples t test, and the results showed that t=18.76>1.96 (α=0.05), p value<0.001, and effect size Cohen's d=6.82 (large effect), which proves that the performance improvement of the present invention has extremely significant statistical significance.

[0060] Example 3: Boundary Parameter Verification Example

[0061] Take the lower limit value of the numerical perturbation step size ε as 10⁻¹ 0 The median value is 10⁻ 8 Maximum value 10⁻ 6In Van der Pol oscillator (μ=10) 6 Complete independent testing on [platform name]:

[0062] ε=10⁻¹ 0 At the time: convergence rate 99.92%, peak memory usage 486MB, single solution time 15.2s;

[0063] ε=10⁻ 8 At the time: convergence rate 99.97%, peak memory usage 486MB, single solution time 12.8s;

[0064] ε=10⁻ 6 At the time: convergence rate 99.85%, peak memory usage 512MB, single solution time 11.5s.

[0065] Test results prove that ε is at 10⁻¹ 0 ~10⁻ 6 Within the scope of this invention, the technical effects of this invention can be achieved, and it is effective in all aspects.

[0066] Example 4: Example of a high degree of freedom system (n=10) 4 )

[0067] 1. Test Case: A 10,000-DOF automotive chassis vertical-lateral coupled vibration system with a system stiffness ratio > 10. 4 Solve for the interval t∈[0,10s];

[0068] 2. Test environment: Same as in Example 1;

[0069] 3. Test Results: The present invention achieved a 100% convergence rate, with a computation time 6.2 times faster than the traditional Newmark-β method, and a vibration response solution error ≤10⁻⁻⁶. 9 Peak memory usage is ≤0.4GB, and there is no divergence in any solution.

Claims

1. A matrix-free, rigidly adaptive numerical calculation method for nonlinear dynamics, characterized in that, Includes the following steps: S1 Matrix-Free Computation Framework Construction: For the nonlinear differential equation system dy / dt=f(y, t) of the nonlinear dynamic system to be solved, a matrix-free differential equation solution framework is constructed. Within the solution step, the directional derivative is calculated. The equivalent calculation of the Jacobian matrix and the vector product is achieved only through vector operations. The directional derivative formula is J(y) \cdot v \approx[f(y+\varepsilon v, t)-f(y, t)] / \varepsilon, where J(y) is the system Jacobian matrix, v is an arbitrary vector, and ε is the numerical perturbation step size, with a value range of 10^{-10} ≤\varepsilon ≤10^{-6}. There is no need to store the complete Jacobian matrix. S2 Quantum Annealing-Driven Spectral Radius Estimation: By using a quantum annealing-driven resource scheduling algorithm, the iteration step size and number of sampling points of the Croft power iteration in matrix-free scenarios are optimized. The spectral radius of the system's Jacobian matrix is ​​estimated in real time through equivalent calculations of finite matrix and vector multiplications. S3 Mixed-precision computing engine optimization: Based on the spectral radius estimation results and the system rigidity characteristics, adaptive mixed-precision floating-point calculation is adopted for different computing stages according to the relative error tolerance of the computing stages; S4 Rigid Adaptive Step Size Control: Based on the spectral radius estimation results, the system stiffness ratio is calculated, and the solution step size and linear solver parameters are dynamically adjusted based on the stiffness ratio to ensure that the solution step size meets the stability constraints and achieve stable convergence solution under extreme stiffness conditions.

2. The matrix-free, rigidly adaptive, nonlinear dynamics numerical calculation method according to claim 1, characterized in that, In step S1, the numerical perturbation step size ε and the system degrees of freedom n satisfy the relationship: \varepsilon=10^{-8} \cdot(n / 1000)^{-0.5}, which eliminates the need to store the complete Jacobian matrix and reduces the space complexity from O(n²) to O(n).

3. The matrix-free, rigidly adaptive, nonlinear dynamics numerical calculation method according to claim 1, characterized in that, In step S2, the specific implementation of quantum annealing-driven spectral radius estimation is as follows: A quantum annealing algorithm is employed to optimize the iteration step size and number of sampling points in the Croft power iteration. In the matrix-free scenario, the spectral radius of the Jacobian matrix is ​​estimated through equivalent calculations of a finite number of matrix-vector multiplications. The objective function of quantum annealing is E(s) = \sum(\rho_{i}^{est}-\rho_{i}^{ref})^{2}, optimizing the spectral radius estimation accuracy by minimizing the sum of squared errors. Here, S is the annealing state vector, \rho_{i}^{est} is the estimated spectral radius value in the i-th iteration, and \rho_{i}^{ref} is the reference baseline value. The annealing scheduling rule adopts an exponential cooling strategy, with the temperature parameter T_k = T_0 \cdot a^k, where T_0 = 100 is the initial temperature, a = 0.95 is the cooling coefficient, and k is the iteration number. The termination condition is that the energy change is <10⁻ for 100 consecutive iterations. 6 Or it can reach the maximum number of iterations of 1000.

4. The matrix-free, rigidly adaptive, nonlinear dynamics numerical calculation method according to claim 1, characterized in that, In step S3, the hybrid precision calculation engine is specifically implemented as follows: For processes with lower accuracy requirements, such as spectral radius estimation and directional derivative calculation, FP16 half-precision floating-point calculation is used; for processes with higher accuracy requirements, such as linear solution and residual update, FP64 double-precision floating-point calculation is used. There is no additional computational overhead for precision switching, and the precision allocation strategy can be dynamically adjusted according to the rigidity characteristics of the system.

5. The matrix-free, rigidly adaptive, nonlinear dynamics numerical calculation method according to claim 1, characterized in that, In step S3, the precision switching rule for mixed precision calculation is as follows: When the relative error tolerance of the calculation stage is >10⁻ 6 When using FP16 half-precision calculation; when the relative error tolerance is ≤10⁻ 6 At that time, FP64 double precision calculation was used.

6. The matrix-free, rigidly adaptive, nonlinear dynamics numerical calculation method according to claim 1, characterized in that, In step S4, the specific implementation method of rigid adaptive step size control is as follows: Based on the spectral radius estimation result ρ(J), the system stiffness ratio R is calculated as R = max |Re(\lambda_{i})| / min |Re(\lambda_{i})|, where \lambda_{i} is the negative real eigenvalue of the Jacobian matrix (Re(λᵢ) < 0). When min |Re(\lambda_{i})| = 0, the stiffness ratio is calculated using the substitution formula R = max |Re(\lambda_{i})| / 10^{-12}. The solution step size h is dynamically adjusted based on the stiffness ratio to satisfy h \cdot \rho(J) ≤ C_{stable}, where C_{stable} is the stability constant with a value range of 0.5 ≤ C_{stable} ≤ 2.0, ensuring the A-stability of the solution process.

7. The matrix-free, rigidly adaptive, nonlinear dynamics numerical calculation method according to claim 1, characterized in that, When steps S1 to S4 are applied to solve a multi-degree-of-freedom rigid vibration system, the system degree of freedom n can be in the range of 1≤n≤10^4, and is compatible with desktop workstations, edge computing, and embedded devices.

8. A matrix-free, rigidly adaptive nonlinear dynamic numerical calculation system, characterized in that, The method according to any one of claims 1 to 7 comprises: The matrix-free computation framework module is used to construct a matrix-free solution framework for nonlinear differential equations, and realizes storage-free equivalent computation of Jacobian matrix and vector product based on directional derivative; The quantum annealing spectral radius estimation module is used to perform Croft power iteration for quantum annealing optimization, thereby achieving spectral radius estimation in matrix-free scenarios. The mixed-precision calculation optimization module is used to perform adaptive mixed-precision calculations based on the system's rigidity characteristics; The rigid adaptive step size control module is used to dynamically adjust the solution step size based on the spectral radius estimation results, so as to achieve stable convergence solution of the extremely rigid system.