High-precision adaptive numerical simulation method and system for solving nonlinear integral-differential equation
Through a unified intelligent control framework, high-precision adaptive numerical simulation of nonlinear integral-differential equations is achieved, which solves the problem of low efficiency of grid adaptation, singular integration and time stepping strategies in existing technologies, improves simulation accuracy and efficiency, and is applicable to fields such as high-energy physics, biology and financial mathematics.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- Filing Date
- 2025-12-16
- Publication Date
- 2026-03-24
AI Technical Summary
Existing simulation software suffers from problems such as insufficient mesh adaptation capability, insufficient singular integral processing capability, low efficiency of time stepping strategy, and numerical instability induced by interpolation method when solving nonlinear integral-differential equations, making it difficult to achieve a balance between simulation accuracy, efficiency, and stability.
By adopting a multi-level adaptive strategy, grid management, integral calculation, time evolution, and interpolation reconstruction are placed under a unified intelligent control framework. Through the coordinated control of time step adaptation, integral accuracy adaptation, and grid space adaptation, the dynamic optimal allocation of computing resources in terms of time scale, integral accuracy, and spatial dimension is achieved.
It significantly improves computational accuracy, efficiency, and stability, solves the "three difficulties" in simulation, and has the ability to solve nonlinear integral-differential equations that can be widely applied in different fields.
Smart Images
Figure CN121723679A_ABST
Abstract
Description
Technical Field
[0001] This invention relates to the field of computer simulation and numerical computation technology, and in particular to a high-precision adaptive numerical simulation method and system for solving nonlinear integral-differential equations. Background Technology
[0002] In scientific computing and engineering simulations in fields such as high-energy physics, fluid mechanics, biology, financial mathematics, and image processing, nonlinear integral-differential equations are key mathematical tools for describing complex systems with memory effects, nonlocal interactions, or long-range correlations. They are particularly suitable for modeling complex physical systems with singular integral kernels, long-range correlations, and memory effects.
[0003] The general form of a nonlinear integral-differential equation is: ; in, The function to be solved is K, which is the integration kernel and often exhibits singularities at certain points (such as when |xy|→0). Analytical solutions to these equations are usually difficult to obtain, and numerical solutions are the only feasible technical approach for system performance prediction, reliability assessment, and optimal design.
[0004] In existing simulation software, solving such equations generally adopts a "discretization + numerical evolution" framework, but its implementation faces the following technical bottlenecks: (1) Lack of mesh adaptive capability: Using a pre-set fixed mesh (such as a uniform mesh or a logarithmic mesh) cannot adaptively adjust according to the dynamic changes of the physical field (such as the temperature field or the stress field) during the simulation. In regions with maximum physical gradients (such as heat sources or crack tips), insufficient mesh resolution will lead to the loss of key information and become the main source of global error.
[0005] (2) Insufficient singular integral handling capability: The integral term is calculated using low-order numerical integration algorithms (such as the trapezoidal rule and Simpson's rule). When the kernel function is singular in the near field, the calculation accuracy drops sharply. In practice, empirical truncation parameters are often forced to be introduced, which destroys the self-consistency of the model and causes the simulation results to deviate from the actual physical behavior.
[0006] (3) The time stepping strategy is inefficient: using a fixed time step cannot adapt to the dynamic changes in the physical process. Redundant calculations are generated in the phase of gradual change, which wastes computing power; in the phase of drastic change, the stability constraints require extremely small steps, which leads to a surge in computing costs and low overall simulation efficiency.
[0007] (4) Numerical instability induced by interpolation methods: Low-order interpolation (such as linear interpolation) is used to reconstruct physical field data at non-grid points. Non-physical numerical oscillations are easily induced in regions with abrupt gradient changes or discontinuities. Errors are amplified globally through integration, which may lead to simulation instability or distortion of results.
[0008] In summary, the existing technologies suffer from defects in their respective algorithms, which are also independent of each other. This leads to a core contradiction in simulation: the difficulty in achieving a balance between accuracy, efficiency, and stability. Therefore, there is an urgent need to provide a high-precision adaptive numerical simulation method and system for solving nonlinear integral-differential equations to address the aforementioned technical problems. Summary of the Invention
[0009] To address the problems existing in the prior art, this invention provides a high-precision adaptive numerical simulation method and system for solving nonlinear integral-differential equations. It aims to achieve a synergistic improvement in computational accuracy, efficiency, and stability by placing grid management, integral calculation, time evolution, and interpolation reconstruction under a unified intelligent control framework through a multi-level adaptive strategy.
[0010] The technical solution of this invention is implemented as follows: A high-precision adaptive numerical simulation method for solving nonlinear integral-differential equations includes the following steps: S1. Set simulation conditions, including the parameters of the objective equation, initial conditions, computational domain, initial trial step size, and ODE tolerance; construct a grid space based on the computational domain; the grid space includes a number of discrete grid points; any point between grid points is denoted as a non-grid point. The non-grid points and the previously defined discrete grid points together constitute a continuous computational domain. ODE stands for Ordinary Differential Equation.
[0011] S2, Based on the solution at the grid points Construct a global cubic spline interpolation function to estimate the position of any point in the computational domain. function value Where n represents the number of iterations and T represents time; This represents the initial conditions; S3. If there are updated grid points, the adaptive numerical simulation integration algorithm is used to calculate the integral term on the right-hand side of the objective equation. Here, "updated" specifically refers to the solution function value of the new time layer obtained by numerical methods at a fixed spatial location (i.e., the grid point coordinates remain unchanged).
[0012] S4. Using the embedded Runge-Kutta method, calculate from arrive The trial solution is obtained, and the estimate of the local truncation error is also obtained. This indicates the trial step size; S5. If the estimated value is greater than the ODE tolerance, then reduce the trial step size. If necessary, repeat steps S3 and S4; otherwise, accept the trial solution as the solution for the next iteration and update the time. At the same time, the trial step size is increased based on the estimated value. By intelligently and dynamically controlling the step size, computational efficiency is effectively improved. While maintaining the same or higher accuracy, the total computation time can be reduced by 30% to 70% compared to the fixed step size used in traditional solutions.
[0013] S6. Dynamic mesh adaptive processing, including checking the gradient and curvature distribution of the solution and adjusting the mesh points; wherein, for regions where the solution changes drastically, the mesh points are densified; for regions where the solution changes smoothly, the mesh points are sparsed; here, "densification" and "sparseness" are also understood as increasing and decreasing.
[0014] The solution of the old grid space is mapped to the adjusted new grid space using the global cubic spline interpolation function, so as to achieve adaptive optimal allocation of computing resources in the spatial direction. S7. Repeat steps S2-S6 until the iteration termination condition is met; the iteration termination condition is either reaching the preset number of iterations or reaching the preset time node, either one is sufficient.
[0015] S8. Output the solution within the computational domain. .
[0016] This invention is the first to place time step adaptation, integral accuracy adaptation, and grid space adaptation within a unified framework. Through a collaborative control mechanism, it achieves dynamic optimal allocation of computing resources in terms of time scale, integral accuracy, and spatial dimension, effectively solving the "triple challenge" of accuracy, efficiency, and stability.
[0017] Before constructing the global cubic spline interpolation function, at the current time step There are N+1 grid points in the grid space. The numerical solution for the grid points has been obtained. The goal is to construct a definition that is present throughout the entire computational domain. Global cubic spline interpolation function on Satisfying: (1) in each subinterval superior, It is a cubic polynomial; (2) exist The point passes through a given data point, i.e. (3) (4) Has continuous first and second derivatives over the entire interval; and In this case, natural boundary conditions are usually adopted, that is... and The construction process includes the following steps: S2-1. Calculate the interval length, i.e., calculate the distance between adjacent grid points. ; S2-2, Construct a system of equations to solve for the second derivative. Let... Let be the second derivative of the spline function at the nodes (an unknown quantity to be determined). For internal nodes... Based on the condition that the first derivative of a spline function is continuous, the following equation can be derived: ; Combined with natural boundary conditions and This constitutes a question about This is a system of N+1 linear equations. Since the coefficient matrix is tridiagonal and strictly diagonally dominant, the system of equations has a unique solution and can be solved directly using an efficient tridiagonal matrix algorithm (such as the Thomas algorithm).
[0018] S2-3, Calculate the spline coefficients for each subinterval. For each subinterval The corresponding cubic polynomial The format is: ;in, .
[0019] After solving for all coefficients, the global cubic spline interpolation function... Defined as a combination of piecewise functions: In this invention, the global cubic spline interpolation function... It can be used for integral calculations in S3 and for obtaining function values for new grid points in S6, i.e., through... To obtain the function value at any point x in space.
[0020] Furthermore, the grid space employs a logarithmic-linear hybrid grid; wherein, for the computational threshold Set a scale transition point ; In small-scale intervals Within, a logarithmic coordinate grid is used, wherein the grid points are positioned according to... Uniform distribution ensures high resolution even when x is extremely small, accurately capturing the singularity of the integral kernel and rapid changes in the solution; in large-scale regions... Within, a uniform linear grid is used, wherein the grid points are positioned according to... It is uniformly distributed to efficiently handle regions where the solution changes relatively gently. Furthermore, the calculation process of the adaptive numerical simulation integration algorithm is as follows: S3-1. Input integration parameters, including the integrand. Integration interval and global tolerance Maximum recursion depth Minimum interval length ; S3-2, Order Let this be the current integration interval. As the initial local tolerance ; S3-3, Perform the following integral recursion process: S3-31. Calculate the integrand using the rules of higher-order integration. Integral value within the current integration interval and error estimation Specifically, higher-order integration rules include the Gauss-Kronrod 7-15 point method or the Simpson integration algorithm, etc.
[0021] S3-32. Recursion Termination Determination: Determine whether any of the following termination conditions are met: The length of the current integration interval is less than or equal to The current recursion depth is greater than or equal to ; If satisfied, then retain the current integration interval. and Otherwise, continue executing S3-33; S3-33. Divide the current integration interval into two sub-intervals, namely the left interval and the right interval, and assign a local tolerance value to each sub-interval. ;like Divide into sub-intervals left= and right= , It should be noted that if c is a singular point of the integrand, then the current integration interval is preferentially divided, that is, new subintervals are created so that the singular point lies within one of these subintervals. Specifically, if the integrand... If point c is undefined, tends towards infinity, or changes drastically, direct integration will lead to decreased accuracy or even calculation failure, then point c is a singular point.
[0022] S3-34. During recursion, the sub-interval is used as the new current integration interval, and the integral value and error estimate of the sub-interval are calculated; wherein, the integral value and error estimate of the left interval are respectively and The integral value and error estimate of the right interval are respectively and The calculation results of the two sub-intervals are combined to obtain the integral value and error estimate of the parent interval, respectively. and ; S3-35, Output Integration Interval The integral value is the right-hand integral term.
[0023] The above methods can be used to calculate integral terms efficiently and with high accuracy. Robust handling of singular kernels, where Ω represents the computational threshold y, i.e., the non-original mesh point in the computational domain. For kernel function, when and When very close, It may tend towards infinity, resulting in the singularity of the integrand. The merging process of S3-35 is a layer-by-layer callback, and each layer... All must be with Compare once to ensure that the error estimate of [a,b] is less than or equal to .
[0024] Since the endpoints in the subinterval partitioning process during integration are not necessarily the discrete points corresponding to the grid points, a global cubic spline interpolation function is used to calculate the integral value.
[0025] An adaptive integral algorithm is used to robustly handle singular kernels, and cubic spline interpolation is combined to provide high-order accurate function estimation, which significantly improves the solution accuracy, especially near singular points and in regions with large gradients.
[0026] Furthermore, the estimated value is calculated as follows: Evaluate values using a shared set of intermediate functions. Construct two approximate solutions of different orders; where, ; ,s represents the series of the Runge-Kutta method; This represents the time node offset coefficient for the i-th level; Indicates intermediate step coefficients, used for calculation and Weighted combination; This indicates the calculation process of the right-hand integral term; The approximate solution includes higher-order trial solutions. and low-order trial solutions An error estimation vector is constructed using two approximate solutions. ,Right now: ; in, and All are intermediate coefficients, including those mentioned above. and All of these can be obtained by querying the Butcher table; The estimated value is ;in, for The j-th component; ATOL and RTOL represent the absolute error tolerance and relative error tolerance, respectively. In S5, the ODE tolerance is expressed as: ;like ,but As the proposed solution.
[0027] Furthermore, in step S6, the mapping step is as follows: S6-1, for each of the grid points Calculate the grid demand density value ,Right now: ; The second derivative of the global cubic spline interpolation function represents the curvature of the grid points. Sensitivity coefficient; S6-2, For any two adjacent grid points forming an interval Calculate the average demand density ,Right now: ; S6-3, For any one of the grid points The average demand density and preset baseline density Comparison: If Then in the interval Insert at least one new grid point to achieve the encryption; if Then the interval Merge the grid points with adjacent grid points that also meet the sparsity condition, and delete all grid points within the merged new grid to achieve the sparsity. This represents the preset threshold, and ;or ; S6-4, The final grid point sequence is represented as follows: , The location of the new grid point is represented by the global cubic spline interpolation function. The solution, i.e. This means completing the mapping from the old grid space to the new grid space.
[0028] Furthermore, in the iterations, the dynamic mesh adaptive processing of S6 is only executed periodically; wherein, when At that time, perform grid adjustment. This indicates the threshold adjustment, and n represents the nth iteration.
[0029] By performing mesh adaptation periodically rather than every step, the computational overhead caused by frequent mesh reconstruction and data structure updates is significantly reduced while ensuring computational accuracy, thus achieving a balance between accuracy and efficiency.
[0030] Furthermore, in S5, if ,but ;like ,but ; in, This represents the estimated value. This indicates the ODE tolerance; and These represent the trial step sizes before and after shrinking / increasing, respectively; and These are used to limit the maximum and minimum values of the trial step size, respectively. This is the error estimate corresponding to the estimated value; and These are the step size growth factor and the step size reduction factor, respectively.
[0031] This invention also provides a high-precision adaptive numerical simulation system for solving nonlinear integral-differential equations, which applies the high-precision adaptive numerical simulation method for solving nonlinear integral-differential equations as described above; it includes the following modules: The initialization module is used to execute S1; The spline interpolation module is used to perform S2; The adaptive integral module is used to execute S3; An adaptive step size control module is used to execute S3 to S5; The dynamic mesh adaptive module is used to execute S6.
[0032] Compared with the prior art, the present invention achieves the following beneficial effects: (1) This invention provides a high-precision adaptive numerical simulation method and system for solving nonlinear integral-differential equations, realizing the collaborative optimal allocation of multi-dimensional computing resources. Specifically, for the first time, time step adaptation, integral accuracy adaptation, and grid space adaptation are placed in a unified framework. Through a collaborative control mechanism, the dynamic optimal allocation of computing resources in terms of time scale, integral accuracy, and spatial dimension is realized, effectively solving the "three difficulties" of accuracy, efficiency, and stability in high-precision numerical solutions.
[0033] (2) An adaptive integral algorithm is adopted, which robustly handles the singularity of the kernel function through singularity identification and interval partitioning strategies. At the same time, combined with the cubic spline interpolation function, it provides continuous and high-order precision function value estimation for integral calculation.
[0034] (3) By comparing error estimation and ODE tolerance, intelligent dynamic control of step size is achieved. Under the premise of meeting global accuracy requirements, the total number of time steps is reduced to the maximum extent, which significantly improves computational efficiency.
[0035] (4) By dynamically adjusting the density of grid points in the grid space through the degree of change of the grid points, the computing resources are adaptively and optimally configured in the spatial direction, thus avoiding the waste of resources caused by uniform grids.
[0036] (5) The present invention provides a framework solution with high versatility. By replacing the integrator kernel K and the right-hand side F, it can be widely applied to solving nonlinear integral-differential equations in different fields, such as the BK equation in high-energy physics, the FKPP equation in biology, and the jump-diffusion model in financial mathematics, etc., and has extremely high application value and promotion prospects.
[0037] (6) The modular design makes the system easy to implement, debug and expand. Unified error control and a wealth of optional configuration parameters (such as tolerance and mesh type) enable it to flexibly adapt to different accuracy and efficiency requirements and have good engineering feasibility. Attached Figure Description
[0038] Figure 1 This is a flowchart illustrating a high-precision adaptive numerical simulation method for solving nonlinear integral-differential equations provided in an embodiment of the present invention. Figure 2 This is a flowchart illustrating the adaptive numerical simulation integration algorithm provided in this embodiment of the invention. Detailed Implementation
[0039] To make the objectives, technical solutions, and advantages of this invention clearer, the technical solutions in the embodiments of this invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only a part of the embodiments of this invention, and not all of them. All other embodiments obtained by those skilled in the art based on the embodiments of this invention without creative effort are within the scope of protection of this invention.
[0040] like Figures 1 to 2 As shown, this embodiment provides a high-precision adaptive numerical simulation method for solving nonlinear integral-differential equations, including the following steps: S1. Set simulation conditions, including the parameters of the objective equation, initial conditions, computational domain, initial trial step size, and ODE tolerance; construct a mesh space based on the computational domain; in this embodiment, a logarithmic-linear hybrid mesh is used; wherein, for the computational threshold... Set a scale transition point In small-scale intervals Within, a logarithmic coordinate grid is used, where the grid points are arranged according to... Uniform distribution ensures high resolution even when x is extremely small, accurately capturing the singularity of the integral kernel and rapid changes in the solution; in large-scale regions... Within, a uniform linear grid is used, and the grid points are arranged according to... It is uniformly distributed to efficiently handle regions where the solution changes relatively gently.
[0041] The grid space comprises a discrete set of grid points; any point between the grid points is designated as a non-grid point. The non-grid points and the existing discrete grid points together constitute a continuous computational domain.
[0042] S2, Solution based on grid points Construct a global cubic spline interpolation function to estimate the function value at any point x in the computational domain. Where n represents the number of iterations and T represents time; This indicates the initial conditions.
[0043] Before constructing the global cubic spline interpolation function, at the current time step There are N+1 grid points in the grid space. The numerical solution for the grid points has been obtained. The goal is to construct a definition that is present throughout the entire computational domain. Global cubic spline interpolation function on Satisfying: (1) in each subinterval superior, It is a cubic polynomial; (2) exist The point passes through a given data point, i.e. (3) (4) Has continuous first and second derivatives over the entire interval; and In this case, natural boundary conditions are usually adopted, that is... and The construction process includes the following steps: S2-1. Calculate the interval length, i.e., calculate the distance between adjacent grid points. ; S2-2, Construct a system of equations to solve for the second derivative. Let... Let be the second derivative of the spline function at the nodes (an unknown quantity to be determined). For internal nodes... Based on the condition that the first derivative of a spline function is continuous, the following equation can be derived: ; Combined with natural boundary conditions and This constitutes a question about This is a system of N+1 linear equations. Since the coefficient matrix is tridiagonal and strictly diagonally dominant, the system of equations has a unique solution and can be solved directly using an efficient tridiagonal matrix algorithm (such as the Thomas algorithm).
[0044] S2-3, Calculate the spline coefficients for each subinterval. For each subinterval The corresponding cubic polynomial The format is: ;in, .
[0045] After solving for all coefficients, the global cubic spline interpolation function... Defined as a combination of piecewise functions: In this invention, the global cubic spline interpolation function... It can be used for integral calculations in S3 and for obtaining function values for new grid points in S6, i.e., through... To obtain the function value at any point x in space.
[0046] S3. If there are updated grid points, use the adaptive numerical simulation integration algorithm to calculate the integral term on the right-hand side of the objective equation. Here, "updated" specifically refers to the solution function value of the new time layer obtained by numerical methods at a fixed spatial location (i.e., the grid point coordinates remain unchanged).
[0047] Specifically, the calculation process of the adaptive numerical simulation integral algorithm is as follows: S3-1. Input integration parameters, including the integrand. Integration interval and global tolerance Maximum recursion depth Minimum interval length These fixed (static) parameters, as well as the fixed parameters discussed later, are also part of the simulation conditions.
[0048] S3-2, Order Let this be the current integration interval. As the initial local tolerance ; S3-3, Perform the following integral recursion process: S3-31. Calculate the integrand using the rules of higher-order integrals. The integral value within the current integration interval and error estimation In this embodiment, the higher-order integration rule is the Gauss-Kronrod 7-15 point method.
[0049] S3-32. Recursion Termination Determination: Determine whether any of the following termination conditions are met: The length of the current integration interval is less than or equal to The current recursion depth is greater than or equal to ; If satisfied, then retain the current integration interval. and Otherwise, continue executing S3-33; S3-33. Divide the current integration interval into two sub-intervals, namely the left interval and the right interval, and assign a local tolerance value to each sub-interval. ;like Divide into sub-intervals left= and right= , It should be noted that if c is a singular point of the integrand, then the current integration interval is preferentially divided, that is, new subintervals are created so that the singular point lies within one of these subintervals. Specifically, if the integrand... If point c is undefined, tends towards infinity, or changes drastically, direct integration will lead to decreased accuracy or even calculation failure, then point c is a singular point.
[0050] S3-34. During recursion, the sub-interval is used as the new current integration interval, and the integral value and error estimate of the sub-interval are calculated; where the integral value and error estimate of the left interval are respectively... and The integral value and error estimate of the right interval are respectively and The calculation results of the two sub-intervals are combined to obtain the integral value and error estimate of the parent interval, respectively. and ; S3-35, Output Integration Interval The integral value is the integral term on the right side.
[0051] The above methods can be used to calculate integral terms efficiently and with high accuracy. Robust handling of singular kernels, where Ω represents the computational threshold y, i.e., the non-original mesh point in the computational domain. For kernel functions, when x and y are very close, It may tend towards infinity, resulting in the singularity of the integrand. The merging process of S3-35 is a layer-by-layer callback, and each layer... All must be with Compare once to ensure that the error estimate of [a,b] is less than or equal to .
[0052] Since the endpoints in the subinterval partitioning process during integration are not necessarily the discrete points corresponding to the grid points, a global cubic spline interpolation function is used to calculate the integral value.
[0053] An adaptive integral algorithm is used to robustly handle singular kernels, and cubic spline interpolation is combined to provide high-order accurate function estimation, which significantly improves the solution accuracy, especially near singular points and in regions with large gradients.
[0054] S4. Using the embedded Runge-Kutta method, calculate from arrive The trial solution is obtained, and the estimate of the local truncation error is also obtained. This represents the trial step size; in this embodiment, the estimated value is calculated as follows: Evaluate values using a shared set of intermediate functions. Construct two approximate solutions of different orders; where, ; ,s represents the series of the Runge-Kutta method; This represents the time node offset coefficient for the i-th level; Indicates intermediate step coefficients, used for calculation and Weighted combination; This indicates the calculation process of the integral term on the right-hand side; Approximate solutions include higher-order trial solutions. and low-order trial solutions An error estimation vector is constructed using two approximate solutions. ,Right now: ; in, and All are intermediate coefficients, including those mentioned above. and All of these can be obtained by querying the Butcher table; The estimated value is ;in, for The j-th component; ATOL and RTOL represent the absolute error tolerance and relative error tolerance, respectively.
[0055] S5. If the estimated value is greater than the ODE tolerance, then reduce the trial step size. If necessary, repeat steps S3 and S4; otherwise, accept the trial solution as the solution for the next iteration and update the time. At the same time, increase the trial step size based on the estimated value. By intelligently and dynamically controlling the step size, computational efficiency is effectively improved. While maintaining the same or higher accuracy, the total computation time can be reduced by 30% to 70% compared to the fixed step size used in traditional solutions.
[0056] Specifically, if ,but ;like ,but ,and As a trial solution, that is, as the solution for the next iteration.
[0057] in, This represents an estimated value. Indicates ODE tolerance; and These represent the trial step size before and after shrinking / increasing, respectively; and These are used to limit the maximum and minimum values of the trial step size, respectively; E is the error estimate corresponding to the estimated value; and These are the step size growth factor and the step size reduction factor, respectively.
[0058] S6. Perform dynamic mesh adaptive processing periodically, including checking the gradient and curvature distribution of the solution and adjusting the mesh points; for regions where the solution changes drastically, densify the mesh points; for regions where the solution changes smoothly, sparse the mesh points; here, "densification" and "sparseness" can also be understood as increasing and decreasing.
[0059] By using a global cubic spline interpolation function, the solution of the old grid space is mapped to the adjusted new grid space, so as to achieve adaptive optimal allocation of computing resources in spatial direction.
[0060] Specifically, the mapping process is as follows: S6-1, for each grid point Calculate the grid demand density value ,Right now: ; The second derivative of the global cubic spline interpolation function represents the curvature of the grid points; Sensitivity coefficient; S6-2, For any interval formed by two adjacent grid points Calculate the average demand density ,Right now: ; S6-3, For any grid point Average demand density and preset baseline density Comparison: If Then in the interval Insert at least one new grid point to achieve encryption; if Then the interval Merge with adjacent intervals that meet the same sparsity condition, and delete all grid points within the merged new interval to achieve sparsity; This represents the preset threshold, and By adjusting these two thresholds, the trigger "thresholds" for mesh refinement (increasing resolution) and sparsity (reducing computational load) can be easily controlled, thereby enabling precise trade-offs and strategy control between computational accuracy and computational cost.
[0061] S6-4, The final grid point sequence is represented as follows: , The location of the new grid point is represented by the global cubic spline interpolation function. The solution, i.e. This means completing the mapping from the old grid space to the new grid space.
[0062] when At that time, perform grid adjustment. The threshold value is denoted by , and n represents the nth iteration. By performing mesh adaptation periodically rather than every step, the computational overhead caused by frequent mesh reconstruction and data structure updates is significantly reduced while maintaining computational accuracy, achieving a balance between accuracy and efficiency.
[0063] S7. Repeat steps S2-S6 until the iteration termination condition is met; the iteration termination condition is either reaching the preset number of iterations or reaching the preset time node, either one is sufficient.
[0064] S8. Output the solution within the computational domain. .
[0065] This invention is the first to place time step adaptation, integral accuracy adaptation, and grid space adaptation within a unified framework. Through a collaborative control mechanism, it achieves dynamic optimal allocation of computing resources in terms of time scale, integral accuracy, and spatial dimension, effectively solving the "triple challenge" of accuracy, efficiency, and stability.
[0066] This embodiment also provides a high-precision adaptive numerical simulation system for solving nonlinear integral-differential equations, including the following modules: The initialization module is used to execute S1; The spline interpolation module is used to perform S2; The adaptive integration module is used to perform the integration calculation of S3; An adaptive step size control module is used to execute S3 to S5; The dynamic mesh adaptive module is used to execute S6.
[0067] Taking the solution of the Balitsky-Kovchegov (BK) equation as an example, the BK equation is a core nonlinear integral-differential equation in high-energy physics describing gluon saturation. It has a complex form and the integral kernel has singularity, making it a typical example to verify the effectiveness of this invention.
[0068] During initialization, the physical parameters set include: initial velocity. =0, terminate speed =10, initial scattering amplitude Other numerical parameters include: The spatial range is That is, corresponding to the computational domain, the initial number of grid points is 1000; the initial number of grid points is 1000; the initial trial step size is... =0.1 (corresponding to) Integral tolerance (i.e. ODE tolerance An initial grid space with respect to r is generated using a log-linear hybrid grid, with transition points set as... (correspond The small-scale interval is... Large-scale intervals are .
[0069] After initialization, the iterative process is executed. During this process, when the adaptive integration module is called to calculate the right-hand integral term of the BK equation, the adaptive integration module automatically handles the... and The core singularity at that location. Within the adaptive step size control module, The value can be dynamically adjusted between 0.01 and 0.5 based on the evolution of the solution. In the dynamic mesh adaptive module, every 5 iterations (corresponding to...) Each iteration checks the curvature of the solution and adjusts the r-grid locally accordingly.
[0070] Compared with traditional methods that employ fixed step size RK4, trapezoidal integration, and linear interpolation, the embodiments of the present invention, in their evolution to... At that time, it has the following significant advantages: in terms of accuracy, the scattering amplitude The computational error was reduced by approximately 60%; in terms of efficiency, the total computation time was reduced by approximately 45%. In terms of stability, no numerical oscillations or divergences were observed throughout the entire computational domain.
[0071] Based on the disclosure and teachings of the foregoing specification, those skilled in the art can make changes and modifications to the above embodiments. Therefore, the present invention is not limited to the specific embodiments disclosed and described above, and some modifications and changes to the present invention should also fall within the protection scope of the claims of the present invention. Furthermore, although some specific terms are used in this specification, these terms are only for convenience of explanation and do not constitute any limitation on the present invention.
Claims
1. A high-precision adaptive numerical simulation method for solving nonlinear integral-differential equations, characterized in that, Includes the following steps: S1. Set simulation conditions, including parameters of the objective equation, initial conditions, computational domain, initial trial step size, and ODE tolerance; construct a grid space based on the computational domain; the grid space includes a discrete number of grid points; S2, Based on the solution at the grid points Construct a global cubic spline interpolation function to estimate the position of any point in the computational domain. function value Where n represents the number of iterations and T represents time; This represents the initial conditions; S3. If there are updated grid points, then use the adaptive numerical simulation integration algorithm to calculate the integral term on the right-hand side of the objective equation. S4. Using the embedded Runge-Kutta method, calculate from arrive The trial solution is obtained, and the estimate of the local truncation error is also obtained. This indicates the trial step size; S5. If the estimated value is greater than the ODE tolerance, then reduce the trial step size. If necessary, repeat steps S3 and S4; otherwise, accept the trial solution as the solution for the next iteration and update the time. At the same time, the trial step size is increased based on the estimated value. ; S6. Dynamic mesh adaptive processing, including checking the gradient and curvature distribution of the solution and adjusting the mesh points; wherein, for regions where the solution changes drastically, the mesh points are densified; for regions where the solution changes smoothly, the mesh points are sparsed. Using the global cubic spline interpolation function, the solution of the old grid space is mapped to the adjusted new grid space; S7. Repeat steps S2-S6 until the iteration termination condition is met; S8. Output the solution within the computational domain. .
2. The high-precision adaptive numerical simulation method for solving nonlinear integral-differential equations according to claim 1, characterized in that, The grid space employs a logarithmic-linear hybrid grid; wherein, for the computational threshold Set a scale transition point ; In small-scale intervals Within, a logarithmic coordinate grid is used, wherein the grid points are positioned according to... Uniform distribution; large-scale intervals Within, a uniform linear grid is used, wherein the grid points are positioned according to... It is uniformly distributed.
3. The high-precision adaptive numerical simulation method for solving nonlinear integral-differential equations according to claim 1, characterized in that, The calculation process of the adaptive numerical simulation integral algorithm is as follows: S3-1. Input integration parameters, including the integrand. Integration interval and global tolerance Maximum recursion depth Minimum interval length ; S3-2, Order Let this be the current integration interval. As the initial local tolerance ; S3-3, Perform the following integral recursion process: S3-31. Calculate the integrand using the rules of higher-order integration. Integral value within the current integration interval and error estimation ; S3-32. Recursion Termination Determination: Determine whether any of the following termination conditions are met: The length of the current integration interval is less than or equal to The current recursion depth is greater than or equal to ; If satisfied, then retain the current integration interval. and Otherwise, continue executing S3-33; S3-33. Divide the current integration interval into two sub-intervals, namely the left interval and the right interval, and assign a local tolerance value to each sub-interval. ; S3-34. During recursion, the sub-interval is used as the new current integration interval, and the integral value and error estimate of the sub-interval are calculated; wherein, the integral value and error estimate of the left interval are respectively and The integral value and error estimate of the right interval are respectively and The calculation results of the two sub-intervals are combined to obtain the integral value and error estimate of the parent interval, respectively. and ; S3-35, Output Integration Interval The integral value is the right-hand integral term.
4. A high-precision adaptive numerical simulation method for solving nonlinear integral-differential equations according to claim 1, characterized in that, The estimated value is calculated as follows: Evaluate values using a shared set of intermediate functions. Construct two approximate solutions of different orders; where, ; ,s represents the series of the Runge-Kutta method; This represents the time node offset coefficient for the i-th level; Indicates intermediate step coefficients, used for calculation and Weighted combination; This indicates the calculation process of the right-hand integral term; The approximate solution includes higher-order trial solutions. and low-order trial solutions An error estimation vector is constructed using two approximate solutions. ,Right now: ； in, and All are intermediate coefficients; The estimated value is ;in, for The j-th component; ATOL and RTOL represent the absolute error tolerance and relative error tolerance, respectively. In S5, the ODE tolerance is expressed as: ;like ,but As the proposed solution.
5. A high-precision adaptive numerical simulation method for solving nonlinear integral-differential equations according to claim 1, characterized in that, In step S6, the mapping step is as follows: S6-1, for each of the grid points Calculate the grid demand density value ,Right now: ; The second derivative of the global cubic spline interpolation function represents the curvature of the grid points. Sensitivity coefficient; S6-2, For any two adjacent grid points forming an interval Calculate the average demand density ,Right now: ; S6-3, For any one of the grid points The average demand density and preset baseline density Comparison: If Then in the interval Insert at least one new grid point to achieve the encryption; if Then the interval Merge the grid points with adjacent grid points that also meet the sparsity condition, and delete all grid points within the merged new grid to achieve the sparsity. This represents the preset threshold, and ; S6-4, The final grid point sequence is represented as follows: , The location of the new grid point is represented by the global cubic spline interpolation function. The solution, i.e. This means completing the mapping from the old grid space to the new grid space.
6. A high-precision adaptive numerical simulation method for solving nonlinear integral-differential equations according to claim 1, characterized in that, In the iteration, the dynamic mesh adaptive processing of S6 is only executed periodically; wherein, when At that time, perform grid adjustment. This indicates the threshold adjustment, and n represents the nth iteration.
7. A high-precision adaptive numerical simulation method for solving nonlinear integral-differential equations according to claim 1, characterized in that, In S5, if ,but ;like ,but ; in, This represents the estimated value. This indicates the ODE tolerance; and These represent the trial step sizes before and after shrinking / increasing, respectively; and These are used to limit the maximum and minimum values of the trial step size, respectively; E is the error estimate corresponding to the estimated value; and These are the step size growth factor and the step size reduction factor, respectively.
8. A high-precision adaptive numerical simulation system for solving nonlinear integral-differential equations, employing a high-precision adaptive numerical simulation method for solving nonlinear integral-differential equations as described in any one of claims 1 to 7; characterized in that, Includes the following modules: The initialization module is used to execute S1; The spline interpolation module is used to perform S2; The adaptive integral module is used to execute S3; An adaptive step size control module is used to execute S3 to S5; The dynamic mesh adaptive module is used to execute S6.
Citation Information
Cited By
A PINN intelligent simulation method for non-continuous ODE systems
CN122365969A