Solving method of ordinary differential equation solving circuit based on RK6-Adams6 hybrid algorithm

Abstract

The invention discloses a solving method of an ordinary differential equation solving circuit based on an RK6-Adams6 hybrid algorithm, and the method comprises the steps: taking an initial input parameter as the input of a function solving module, executing an ordinary differential equation function value solving task stage by stage through the function solving module, and transmitting a solved function value to a high-precision algorithm module, an RK6-Adams6 hybrid algorithm is iteratively executed by using the high-precision algorithm module, and then an iteration result is fed back to the function solving module, so that cyclic interactive iterative operation is performed stage by stage, and a final numerical solution of an ordinary differential equation is obtained; according to the technical scheme, the requirement for high-precision solution of a complex ODE system can be met, the global error fluctuation range is reduced in a self-adaptive dynamic step length adjustment mode, and the stability of the process of iteratively solving an ordinary differential equation is improved; in the iteration process, a function solving module is used for optimizing the operation delay of a key operation link according to the dependency relationship among a plurality of operation chains of the equation, and the solving speed of the equation is increased.

Description

Technical Field

[0001] The present application relates to the technical field of high-performance computing hardware design, and in particular to a method for solving ordinary differential equations based on an RK6-Adams6 hybrid algorithm. Background Art

[0002] Ordinary differential equations (ODEs) have been widely used in fields such as aerospace, neuroscience, and industrial simulation. They describe the dynamic patterns of dynamic systems and complex processes, providing a theoretical foundation for scientific exploration and technological development in multiple fields. However, since most ODEs cannot be directly and accurately solved analytically, numerical methods have become the core means of solving ODEs. With the continuous deepening of scientific and technological research, the demand for accuracy, flexibility, and speed in solving ODEs continues to grow. At the same time, the mathematical calculation problems faced in multiple fields have become increasingly complex.

[0003] Traditional numerical solvers are software-based, but their computational power and speed are insufficient for complex and computationally intensive ODE systems. Circuit-based numerical solvers are an effective means of addressing these issues. In such engineering applications, FPGAs are often used as the circuit implementation, significantly improving computational speed when handling large numbers of complex ODE systems. However, while hardware circuits must achieve rapid solution speeds, resource utilization still needs to be improved.

[0004] Numerical methods for solving ordinary differential equations are generally divided into single-step and multi-step methods. The former include the Euler method and the Runge-Kutta method, while the latter include the Hamming method and the Adams method. Existing numerical solutions for ODEs often rely on low-order or single numerical algorithms, which are unable to meet the accuracy requirements for high-order, high-performance ODE solutions.

[0005] Therefore, a circuit design method that balances fast solution efficiency and high solution accuracy is urgently needed. This method can not only meet the requirements of high-precision and fast solution of complex ODE systems, but also improve hardware resource utilization. This method can break through the accuracy and speed limitations at the circuit level and provide key technical support for high-real-time simulation. Summary of the Invention

[0006] The purpose of this application is to realize the circuit design of numerical solutions to high-order ordinary differential equations, while taking into account the requirements of fast solution efficiency and high solution accuracy, improving the utilization rate of hardware resources, so that the circuit design of numerical solutions can be applied to the solution of high-order ordinary differential equations in scenarios such as aerospace, neuroscience, and industrial simulation.

[0007] The technical solution of the present application is to provide a solution method for an ordinary differential equation solution circuit based on the RK6-Adams6 hybrid algorithm, wherein the solution circuit includes a data transmission control module, a data cache module, an operator state management module, an adaptive step size control module, a high-precision algorithm module, a function solution module, and a data path control module. The method includes:

[0008] S1: The host computer transmits the initial input parameters of the ordinary differential equation to the data cache module through the data transmission control module, and the data cache module starts to cache the initial input parameters;

[0009] S2: After the cache is completed, the main state machine of the operator state management module jumps to the RK6 main state, and the data cache module transmits the initial input parameters to the function solver module;

[0010] S3: The data path control module schedules the iteration process of the RK6 algorithm. According to the jump state of the RK6 slave state machine, the function solving module and the high-precision algorithm module perform cyclic interactive iterative operations in each stage to obtain the initial solution of the Adams6 algorithm. Specifically:

[0011] In a single iteration, the RK6 slave state machine jumps to the slave state in the order of RK1 to RK7. Each jump is an algorithm stage. In a single algorithm stage, the data path control module transmits the function value output by the function solving module to the high-precision algorithm module. The adaptive step size control module configures the corresponding step size according to the algorithm stage. The high-precision algorithm module calculates the slope corresponding to the function value based on the corresponding step size and the function value to obtain the initial solution, and feeds the calculated result parameters back to the function solving module. The function solving module updates the function value based on the current result parameters. This cycle iterates until a predetermined number of initial solutions are obtained.

[0012] S4: After the iteration is completed, the main state machine jumps to the Adams6 main state. The data path control module schedules the iteration process of the Adams6 algorithm. According to the jump state of the Adams6 slave state machine, the function solving module and the high-precision algorithm module are controlled to perform cyclic interactive iterative operations stage by stage to obtain the final numerical solution of the ordinary differential equation.

[0013] S5: After the iteration is completed, the operator state management module jumps to the idle state, the data path control module terminates the interactive task, and the data cache module uploads the final numerical solution of the ordinary differential equation to the host computer through the data transmission control module.

[0014] Furthermore, the operator state management module also includes an RK6 slave state machine, and the S3 specifically includes:

[0015] S31: The function solving module brings the initial input parameters into the ordinary differential equation and calculates the function value y of the ordinary differential equation n, RK6 slave state machine jumps from idle state to RK1 slave state;

[0016] S32: The data path control module transmits the current function value to the high-precision algorithm module. The adaptive step size control module configures the step size for the current algorithm stage. The high-precision algorithm module calculates the slope K2 based on the current step size and the function value. The data path control module transmits the result parameter corresponding to the slope K2 to the function solving module. The function solving module uses the current result parameter as a new initial value to update the function value. The RK6 slave state machine jumps from the RK1 slave state to the RK2 slave state.

[0017] S33: Perform interactive calculations stage by stage in the manner of S32 until a total of 7 slope values ​​from K1 to K7 are obtained, and RK6 jumps from the state machine to the RK7 slave state;

[0018] S34: The data path control module transmits the current function value to the high-precision algorithm module, the adaptive step size control module configures the step size, and the high-precision algorithm module calculates the initial solution y of the Adams6 algorithm based on the current step size, function value, and each slope value. n+1 , the data path control module transmits the initial solution to the function solving module, and the function solving module uses the initial solution as a new initial value to update the function value;

[0019] S35: The RK6 slave state machine jumps from the RK7 slave state to the RK1 slave state, and repeats the process from S32 to S34 until the RK6 algorithm iteration is completed and a predetermined number of initial solutions are obtained.

[0020] Furthermore, the calculation formula of the RK6 algorithm in S3 is expressed as:

[0021] K1=f(t n ,y n )

[0022]

[0023] Where K1 to K7 are the slope estimates of different algorithm stages within the time step of the RK6 algorithm, t n is the current time point, y n At time t n The function value at y n+1 is time t n+1 The function value at , h is the time step, and f(*) is the slope equation of time and function value.

[0024] Furthermore, in S32, the adaptive step size control module configures the step size for the current algorithm stage, specifically including:

[0025] Calculate the current configuration phase time t n The absolute value of the difference from the dynamic real-time time t|tn -t|, compare|t n -t| and configure the minimum stage step size h n-min The size of , if |t n -t|≥h n-min , then the configuration step size is the preset step size h corresponding to the current algorithm stage n , if |t n -t| <h n-min , then the configuration step size is adjusted to |t n -t|.

[0026] Furthermore, the operator state management module also includes an Adams6 slave state machine, and the S4 specifically includes:

[0027] S41: The data path control module transmits the initial solution cached in the data cache module to the high-precision algorithm module, and the Adams6 slave state machine jumps from the idle state to the Ad1 slave state;

[0028] S42: The data path control module transmits the current function value to the high-precision algorithm module. The adaptive step size control module configures the step size. The high-precision algorithm module calculates the intermediate prediction value of the ordinary differential equation based on the obtained initial solution, the current step size, and the function value. and will Feedback to the function solving module, the function solving module will Update the function value y as the new initial value n , Adams6 slave state machine jumps from state Ad1 to state Ad2;

[0029] S43: Interactive calculation is performed again in the manner of S42, based on the obtained initial solution, current step size, function value and Compute numerical solutions to ordinary differential equations;

[0030] S44: The Adams6 slave state machine jumps from state Ad2 to state Ad1, and repeats the process from S42 to S43 until the Adams6 algorithm iteration is completed and the final numerical solution of the ordinary differential equation is obtained. The Adams6 slave state machine jumps from the Ad2 slave state to the idle state.

[0031] Furthermore, the calculation formula of the Adams6 algorithm in S4 is expressed as:

[0032]

[0033] Where, is the intermediate prediction value of the Adams prediction formula, which is used in the subsequent correction steps. is the final numerical solution after correction.

[0034] Furthermore, the initial input parameters include the expected flight time of the trajectory, the three components of the initial trajectory position, the three components of the initial trajectory velocity, the pressure-related parameters, and the self-inertia force-related parameters. In S32, the result parameters are the same as the parameter types included in the initial input parameters. The function solving module uses the result parameters fed back by the high-precision algorithm module as new initial values, and brings them into the ordinary differential equation to calculate the ordinary differential equation function value used in the next algorithm stage.

[0035] Furthermore, the function solving module includes multiple operation chains, each operation chain is used to calculate different parameters of ordinary differential equations. The function solving module takes the operation chain with the longest time as the critical path, merges the equation blocks that meet the merging conditions in the critical path, and calculates other non-critical paths in parallel with the critical path. The intermediate results that meet the sharing conditions in the operation chain are shared by multiple paths to reduce the operation cycle and scheduling overhead. The merging conditions are that the input values ​​are the same, there is no dependency between each other, and the subsequent paths are crossed. The sharing condition is that the intermediate result is used by the equation blocks in the subsequent link for more than or equal to 2 times.

[0036] During the research and development process, the inventors of the present application discovered that the numerical solution would be unstable if the RK6-Adams6 hybrid algorithm was executed in a stage-by-stage interactive iterative operation manner. The reason is the problem of step size control during the solution process. Different numerical algorithms have different requirements for solution accuracy and stability. The RK6 algorithm tends to adopt a larger step size to improve computational efficiency, while the Adams6 algorithm needs to adopt a smaller step size to maintain numerical stability. Therefore, in the scenario where the two algorithms collaborate in computing, there is discontinuity in the step size switching between different stages, especially near the algorithm transition node, where the step size suddenly switches from a larger value to a smaller value or vice versa, which can easily lead to data instability and rapid propagation of errors during the iteration process, seriously affecting the overall stability and accuracy of the numerical solution of ordinary differential equations. Currently, there is still a lack of a mechanism for adaptive step-size adjustment for algorithm transition nodes, resulting in the inability of step-size configuration at key nodes to effectively adapt to dynamic changes. The present invention proposes an adaptive step-size control mechanism that combines the characteristics of the algorithm stage, which can dynamically adjust the step-size according to the distance between the current stage and the transition node, avoiding data instability and error propagation problems caused by sudden step-size changes when switching between different algorithms. The mechanism can efficiently complete the step-size configuration through only three basic operators: subtraction, absolute value, and comparison. It has the advantages of simple structure, fast response, and strong adaptability, and significantly improves the stability and accuracy of the numerical solutions of ordinary differential equations.

[0037] The beneficial effects of this application are:

[0038] First, the technical solution in this application is based on FPGA to set up an adaptive step control module, a high-precision algorithm module, and a function solving module. The high-precision algorithm module and the function solving module are used to execute the RK6-Adams6 hybrid algorithm through stage-by-stage interactive iterative operations to solve ordinary differential equations, which can meet the needs of high-precision solution of complex ODE systems. During the execution of the algorithm, the adaptive step control module is used to adaptively and dynamically adjust the step size according to the divided algorithm stages. Since the RK6 algorithm requires a larger step size to be configured, while the Adams6 algorithm requires a smaller step size to be configured, data instability and error propagation problems will occur at the time nodes of the transition between the two algorithms. To solve this problem, the technical solution of this application uses the method of comparing the minimum stage step size to configure the step size, calculates the absolute difference between the current configuration stage time and the dynamic real-time time, compares the absolute difference with the configured minimum stage step size, and then dynamically adjusts the step size according to the comparison result to avoid data instability and error propagation during node transition, narrow the global error fluctuation range, and improve the stability of the numerical solution of the iteratively solved ordinary differential equation.

[0039] Second, the technical solution in this application also sets up a data cache module based on FPGA. The data cache module adopts a storage reuse mechanism to cache all intermediate calculation results of the algorithm operation. At the same time, the high-precision algorithm module adopts a reconfigurable method to realize the configuration of efficient algorithm logic paths, reducing the number of operators and the logical complexity of the circuit, which can meet the needs of high-efficiency solutions.

[0040] Third, the computing speed is significantly improved: the function solving module optimizes the operation delay of key operation links based on the dependencies between multiple operation chains of the equation, merges the equation blocks that meet the merging conditions to reduce the operation cycle and scheduling overhead, and improves the equation solving speed to meet the needs of high-efficiency solution. BRIEF DESCRIPTION OF THE DRAWINGS

[0041] The advantages of the above and / or additional aspects of the present application will become apparent and readily understood from the description of the embodiments in conjunction with the following drawings, in which:

[0042] Figure 1 1 is a schematic diagram of a workflow of an ordinary differential equation solver based on the RK6-Adams6 hybrid algorithm according to an embodiment of the present application;

[0043] Figure 2 1 is a schematic diagram of a circuit structure for solving ordinary differential equations based on the RK6-Adams6 hybrid algorithm according to an embodiment of the present application;

[0044] Figure 3 This is a schematic diagram of the state jump of the operator management module according to an embodiment of the present application;

[0045] Figure 4 1. It is a schematic diagram of the task allocation and workflow of the RK6-Adams6 hybrid algorithm according to an embodiment of the present application;

[0046] Figure 5 This is a schematic diagram of the task allocation and computational work structure of a six-degree-of-freedom spatial motion equation function solving module according to an embodiment of the present application;

[0047] Figure 6 This is a schematic diagram of the logic path structure of the RK6-Adams6 hybrid algorithm according to an embodiment of the present application;

[0048] Figure 7 This is a schematic diagram of a circuit structure of an adaptive step size control module for solving the problem of adjacent phase nodes according to an embodiment of the present application;

[0049] Figure 8 It is a schematic diagram of the storage reuse mechanism structure of the data cache module according to an embodiment of the present application. DETAILED DESCRIPTION

[0050] In order to more clearly understand the above-mentioned objectives, features and advantages of the present application, the present application is further described in detail below in conjunction with the accompanying drawings and specific embodiments. It should be noted that the embodiments of the present application and the features therein can be combined with each other in the absence of conflict.

[0051] In the following description, many specific details are set forth to facilitate a full understanding of the present application. However, the present application may also be implemented in other ways different from those described herein. Therefore, the scope of protection of the present application is not limited to the specific embodiments disclosed below.

[0052] In this embodiment, taking the six-degree-of-freedom spatial motion equation as an example, the solution method of the ordinary differential equation solving circuit based on the RK6-Adams6 hybrid algorithm is explained. The six-degree-of-freedom spatial motion equation can describe the entire motion state of the carrier, which is usually expressed in the form of a group of ordinary differential equations. The FPGA-based ordinary differential equation solver of this application can adopt the RK6-Adams6 hybrid algorithm to realize the numerical solution operation of the six-degree-of-freedom spatial motion equation to meet its high-precision and fast solution requirements.

[0053] like Figure 1 、 Figure 2 and Figure 4As shown, this embodiment provides a solution method for an ordinary differential equation solving circuit based on the RK6-Adams6 hybrid algorithm. The ordinary differential equation solving circuit is designed based on an FPGA and includes a data transmission control module, a data cache module, an operator state management module, an adaptive step size control module, a high-precision algorithm module, a function solving module, and a data path control module. The operation flow of the solution method is as follows:

[0054] S1: The host computer transmits the initial input parameters of the ordinary differential equation to the data transmission control module through the high-speed interface communication protocol (PCIe 3.0×8 interface protocol). The main state machine of the operator state management module jumps from the idle state to the data uplink and downlink state, and the data cache module begins to cache the initial input parameters. The initial input parameters include nine parameters: the expected flight time of the trajectory, the three components of the initial trajectory position, the three components of the initial trajectory velocity, pressure-related parameters, and inertial force-related parameters.

[0055] like Figure 3 As shown, the operator state management module includes a master state machine and two slave state machines. The master state machine is used to control the overall iteration process of the RK6-Adams6 hybrid algorithm, and the two slave state machines are used to control the iteration process of the RK6 algorithm and the Adams6 algorithm respectively.

[0056] The main state machine includes idle state, data uplink and downlink state, RK6 main state (i.e., RK6 algorithm iteration state), and Adams 6 main state (i.e., Adams 6 algorithm iteration state). Among them, the idle state is used to control circuit reset and initialize data; the data uplink and downlink state is used to control data interaction between the solution circuit and the host computer; the RK6 main state is used to control the iterative loop operation of the RK6 algorithm. The RK6 algorithm includes 5 iterations, each iteration has 7 algorithm stages, and a total of 35 algorithm stages, which can provide 5 initial solutions for subsequent Adams6 algorithm iterations; the Adams6 main state is used to control the iterative loop operation of the Adams6 algorithm. The iteration of the Adams6 algorithm includes controlling explicit prediction and implicit correction iterative calculations.

[0057] The RK6 slave state machine includes an idle state, an RK1 slave state, an RK2 slave state, an RK3 slave state, an RK4 slave state, an RK5 slave state, an RK6 slave state, and an RK7 slave state. The idle state is used to control the circuit reset and start the operation. The RK1 slave state is used to control the operation of K2 in the RK6 operation formula. The RK2 slave state is used to control the operation of K3 in the RK6 operation formula. The RK3 slave state is used to control the operation of K4 in the RK6 operation formula. The RK4 slave state is used to control the operation of K5 in the RK6 operation formula. The RK5 slave state is used to control the operation of K6 in the RK6 operation formula. The RK6 slave state is used to control the operation of K7 in the RK6 operation formula. The RK7 slave state is used to control the y in the RK6 operation formula. n+1 Operation.

[0058] The Adams6 slave state machine includes the idle state, Ad1 slave state, and Ad2 slave state. The idle state is used to control the circuit reset and start the operation, and the Ad1 slave state is used to control the Adams6 operation formula. Predictive operation, Ad2 is used to control the y in the Adams6 operation formula from the state n+1 Correction operation.

[0059] In this embodiment, the master state machine and the two slave state machines are state control logic for the numerical calculation circuit (i.e., ordinary differential equation solving circuit), which is essentially a timing control module in the hardware circuit. Specifically, they can be implemented in FPGA using a finite state machine (FSM, Finite State Machine, which divides the system operation process into several states, performs state jumps according to input conditions under clock drive, can accurately control multi-stage algorithm iterations, and is suitable for data flow and control flow scheduling, such as the staged iterative control of RK6 and Adams6). They can be modeled and implemented using a hardware description language (such as Verilog or VHDL) to drive the orderly execution of the entire numerical solution process.

[0060] S2: After the data cache module completes the caching of all initial input parameters, the main state machine jumps from the data uplink and downlink state to the RK6 main state, and the data cache module transmits the initial input parameters to the function solving module.

[0061] Specifically, after the main state machine enters the RK6 main state, the data path control module obtains initial input parameters from the data cache module and transmits the initial input parameters to the function solving module.

[0062] S3: The data path control module starts scheduling the iteration process of the RK6 algorithm. According to the jump state of the RK6 slave state machine, the function solving module is controlled to execute the ordinary differential equation function value solving task stage by stage, and the solved function value is transmitted to the high-precision algorithm module. The high-precision algorithm module is controlled to iteratively calculate the initial solution of the Adams6 algorithm. The data cache module caches all intermediate calculation results. The calculation formula of the RK6 algorithm includes:

[0063] K1=f(t n ,y n )

[0064]

[0065] Where K1 to K7 are the slope estimates of different algorithm stages within the time step of the RK6 algorithm, t n is the current time point, y n At time t n The function value at (the function value is the state variable that changes with time in the ordinary differential equation), y n+1 is time t n+1 where h is the time step (i.e., the interval between adjacent time points), and f(*) is the slope equation of the function value and time, which is used to describe the change of the state variable in the selected ordinary differential equation over time.

[0066] In this embodiment, the RK6 algorithm and the Adams algorithm are combined to form a hybrid structure in which the RK6 algorithm provides the initial solution and the Adams algorithm completes the subsequent calculations, so as to meet the calculation requirements of high precision and high stability of the final numerical solution. Since Adams6 is a multi-step method, a linear combination of 6 slope values ​​is required during the operation, 5 of which are reused results of the first five steps, and the remaining value is calculated by the function value and time point output by the current function solving module; therefore, the RK6 single-step method is required to provide it with an initial starting value (i.e., an initial solution), wherein RK6 has 7 algorithm stages for each iteration, and 5 iterations total 35 stages, which can provide 5 starting values ​​for Adams6.

[0067] The data path control module schedules the iteration process of the RK6 algorithm to obtain the initial solution of the Adams6 algorithm. Specifically, the following steps are included:

[0068] S31: The function solving module includes multiple operation chains, each of which is used to calculate different parameters of the ordinary differential equation. The function solving module brings the initial input parameters into the operation chain for operation, and combines the final operation results and outputs them to obtain the function value y of the ordinary differential equation. n , cache the function value to the data cache module, and the RK6 slave state machine jumps from the idle state to the RK1 slave state;

[0069] The ordinary differential equations include the center-of-mass dynamics equation of the launching system and the center-of-mass dynamics equation of the projectile system. The center-of-mass dynamics equation is expressed as:

[0070]

[0071] Where x, y, and z are the three components of the initial trajectory position, m is the mass of the projectile, and v is x 、v y 、v z is the three components of relative velocity, G B is the direction cosine array from the missile system to the ground launch system, G V is the direction cosine array from the velocity system to the ground launch system, g ωe 、ω e 、g' r , r are gravitational parameters, ω ex 、ω ey 、ω ez ω e The three components of a ij The matrix is ​​the centrifugal inertia force related matrix, b ij The matrix is ​​the Coriolis inertial force correlation matrix, R 0x 、R 0y 、R 0z is the three components of position offset, P e is a scalar term related to gravitational potential energy, Y 1c , Z 1c is the coordinate component of the center of mass of the projectile in the y and z directions, x 1e is the position component of the projectile in the x direction in the ground launch system, ω Tz1 、ω Ty1 is the angular velocity of the projectile around the z and y axes, C x is the axial force coefficient, is the derivative of the axial force coefficient with respect to the angle of attack α, q is the dynamic pressure, S M is the reference area, β is the sideslip angle, and subscript 1 is the reference coordinate system number.

[0072] The dynamic equation around the center of mass is expressed as:

[0073]

[0074] Where, I x , I y , I z To simplify the three components of the inertia tensor, ω Tx1 、ω Ty1 、ω Tz1 are the three components of acceleration rotating around the center of mass, is the derivative of the lateral moment coefficient with respect to the sideslip angle β and the angle of attack α, l kis the projectile length, is the dimensionless angular velocity, is the damping moment coefficient, R' is the control arm coefficient, r c 、x c 、x g is the geometric size parameter, δ r , δ ψ 、 is the rudder deflection angle.

[0075] In this embodiment, the operation chain refers to a link composed of a series of equation blocks with close dependencies that need to be executed sequentially (each equation block can be regarded as an operation step, and the functions of these equation blocks are implemented based on the circuits in the hardware FPGA). For the ordinary differential equation of the six-degree-of-freedom spatial motion equation, there are many parameters that need to be calculated, and the operation must be performed together by multiple operation chains to obtain the final result.

[0076] In this embodiment, the function solving module can sort out multiple operation chains of ordinary differential equations based on expressions, parallelize independent operation chains, and share the intermediate operation results of each equation block across multiple paths, avoiding repeated calculations, reducing operation cycles and scheduling overhead, optimizing the delay of key operation chains, and compressing the total execution time while ensuring dependencies, thereby improving solution efficiency. Specifically, this includes:

[0077] The function solver module includes multiple operation chains, each used to calculate different parameters of ordinary differential equations. The function solver uses the longest operation chain as the critical path (the critical path determines the shortest solution time; other non-critical paths can be run in parallel, but only the critical path cannot be skipped). Equation blocks in this critical path that meet the merging conditions are merged to form compound operation equation blocks (i.e., equation blocks that can output multiple parameter results). This reduces operation cycles and scheduling overhead. The merging conditions refer to the same input values, no dependencies between them, and cross-links (cross-links refer to the fact that the output values ​​of the operation chains containing multiple equation blocks are merged into the same equation block). Other non-critical paths are run in parallel with the critical path (i.e., the links between the input and output nodes of the non-critical path are aligned in time with the critical path and run in parallel). Intermediate results in the operation chain that meet the sharing conditions are shared across multiple paths to significantly reduce repeated calculations and optimize scheduling delays. The sharing condition is that the intermediate result is used by equation blocks in subsequent chains at least twice.

[0078] like Figure 5As shown, for the above-mentioned six-degree-of-freedom spatial motion equation, the function solving module can perform operations based on the optimized operation chain, split and sort out the equation operation tasks, and use adders / subtractors, dividers, multipliers, fixed-float / floating-point converters, open operators, and trigonometric function solvers to complete the function solving tasks. According to the dependencies between multiple operation chains, multiple equation blocks are calculated in parallel to finally obtain the function value. This process optimizes the operation delay of the key operation link and improves the solution speed. The dependency between each equation block and the number of cycles of pipeline operation are marked in the figure. It can be seen from the figure that in the pipeline operation process of the function solving module, multiple parameters can be calculated in parallel, and each parameter operation process can reuse multiple operators, reducing the pipeline calculation length and the waste of operator resources. Taking the center of mass dynamics equation as an example for parameter operation, the specific process is as follows:

[0079] Task 1: Calculate the atmospheric velocity v and convert the three components of the relative velocity v x 、v y 、v z As the input value for calculation, the calculation equation block is: The calculation of this equation block calls 3 multipliers, 2 adders and 1 square root operator, which requires 46 cycles of operation. Figure 5 This corresponds to the output result of the 53rd cycle of the pipeline.

[0080] Task 2: Calculate the velocity inclination angle θ and the track yaw angle σ. The velocity inclination angle θ and the atmospheric velocity v are calculated in parallel. Both use the three components of the relative velocity as input values. The velocity inclination angle θ is calculated first, and the track yaw angle σ needs to be calculated after obtaining the atmospheric velocity v. The calculations of the velocity inclination angle θ and the track yaw angle σ meet the merging conditions and are merged. The calculation equations for the two are as follows:

[0081]

[0082] In this equation block, the velocity inclination angle θ is calculated first, calling a divider and a CORDIC trigonometric function solver. Figure 5 The corresponding pipeline output is the 74th cycle. After the atmospheric velocity v is solved, the equation block is input to calculate the track yaw angle σ, calling a divider and a CORDIC trigonometric function solver. Figure 5 The corresponding pipeline is output in the 118th cycle; at the same time, in order to facilitate subsequent calculations, after outputting the velocity inclination θ and the track yaw angle σ, the equation block calls the CORDIC trigonometric function solver again to convert the two parameters into trigonometric functions in advance, corresponding to Figure 5 The middle pipeline calculates cycles 106 and 150.

[0083] Task 3, calculate the mass m, according to the set fixed coefficient mass second consumption As well as the input initial projectile mass m0 and trajectory expected flight time t, the calculation equation block is: Due to the existence of a fixed coefficient, compared to the start time of Task 1 and Task 2, Task 3 starts calculation in the third cycle of pipeline calculation, calling a multiplier and an adder (subtraction) in the calculation. Figure 5 Pipeline output at cycle 27.

[0084] Task 4, centrifugal inertia coefficient term (a ij matrix) and the centrifugal inertia coefficient term (b ij The matrix is ​​calculated in parallel. Due to the long calculation cycle of mass m, the three components of displacement z, y, z and the three components of position offset R are simultaneously calculated in the calculation stage. 0x 、R 0y 、R 0z Parallel calculation of x+R 0x 、y+R 0y 、z+R 0z Three parameters, after Task 3 completes the calculation of mass m, start the parallel calculation of centrifugal inertia force coefficient term and centrifugal inertia force coefficient term;

[0085] Among them, a ij and b i The matrix coefficients of x have been mapped in the equation block circuit to facilitate the calculation of the spatial motion equations for various initial parameter conditions. The equation block for the centrifugal inertia coefficient term requires the calculation of 19 multipliers, 9 adders, and a ij Matrix and x+R 0x 、y+R 0y 、z+R 0z In the multiplication calculation, parallel operation is used, and the 9 coefficients of the matrix are multiplied with the three parameters at the same time to reduce the delay cycle. Figure 5 The centrifugal inertia coefficient is output in the 67th cycle of pipeline calculation; It is the differential operation value of the velocity component of the center of mass dynamics equation. The operation of the equation block of the centrifugal inertia coefficient term is the same as that of the centrifugal inertia coefficient term. It requires 13 multipliers and 6 adders to complete the operation of the equation block. Figure 5 The pipeline calculates the output results in the 35th cycle.

[0086] Task 5, gravitational parameter g ωe 、ω e 、g' r , r calculation needs to be completed in stages:

[0087] The first step is to calculate the geocentric radius coefficient r, using x+R 0x 、y+R 0y 、z+R 0zThe three parameters are input values, and the calculation equation block is: It is necessary to call 2 adders, 3 multipliers and 1 square root operator. Since r is used by multiple subsequent equation blocks, it is shared by multiple paths. After the calculation is completed, in order to facilitate the calculation of subsequent steps, the coefficient is reciprocally calculated, corresponding to Figure 5 Output results in the 78th cycle.

[0088] The second step is to calculate the gravitational term related parameter block. The calculation equation block is:

[0089]

[0090] This step will split the calculation into front and back, because the calculation result of the first step of the center of gravity coefficient r only affects Therefore, after the first step is started, another one is calculated in parallel, and the calculation is completed in the 65th cycle of the corresponding pipeline, while the previous part is completed in the 86th cycle. After the current part is completed, the gravity term calculation only needs to call one adder, and the calculation is completed in 4 cycles, corresponding to Figure 5 Output the results in the 90th cycle.

[0091] from Figure 5 It can be seen that Task 4 and Task 5 are in different operation chains. When they are executed, they are calculated step by step according to their respective dependent equation blocks. Since there is no dependency between Task 4 and Task 5, they can be executed in parallel, reducing the pipeline calculation length.

[0092] Task 6, the calculation of the slip angle β and the power angle α, requires the calculation results of the velocity inclination angle θ and the track yaw angle σ in Task 2. The calculation equation block is as follows Figure 5 As shown, the slip angle β needs to be solved by calling 4 adders, 7 multipliers and 6 CORDIC trigonometric functions.

[0093] Task 7, C matrix (G B Matrix) calculation, complete the calculation of the 9 coefficient values ​​in the matrix, and use the calculation results of the velocity inclination and track yaw angle in Task 2 as input values. The calculation equation block is as follows Figure 5 As shown, the pitch angle Yaw angle ψ, roll angle γ. In this part of the calculation, the roll angle γ is mapped to 0, that is, sinγ=0, cosγ=1, which is not shown in the figure. The 9 coefficient values ​​are calculated in parallel. The same items are called for calculation to reduce the operation cycle and resource consumption. For example, the sinψsinγ item in the matrix G B12 and G B22 The coefficients are all calculated. This part requires 12 multipliers and 4 adders (subtraction). Figure 5 The pipeline outputs the result in cycle 222.

[0094] Similarly, the B matrix (G V The calculation of matrix) is consistent with that of C matrix, and the input parameters are the same, where the roll angle γ is mapped to 0, and the B matrix is ​​mapped to Figure 5 The pipeline outputs the result in the 156th cycle.

[0095] Task 8, Aerodynamic R and Thrust F c Calculation is performed based on the B and C matrix calculation results of Task 7. The calculation equation block is as follows Figure 5 As shown, where the aerodynamic R term and thrust F c The calculation of the term is essentially the multiplication of matrix and vector, that is, the multiplication of 3*3 matrix and 3*1 vector. The same circuit structure can be designed to complete the corresponding calculation by inputting different variables. In the calculation of 3*3 matrix and 3*1 vector, 12 variables are designed as configuration items to meet the calculation of different coefficient variables of aerodynamic R term and thrust term. In this structure, 9 multipliers and 6 adders are required. c Item in Figure 5 The output result of the corresponding pipeline 231 cycle, R item is Figure 5 The corresponding pipeline output results in the 166th cycle.

[0096] Task 9, the final calculation of the center of mass dynamic equation involves adding the centrifugal inertia force term and the Coriolis inertia force term from Task 4, the gravity term from Task 5, and the aerodynamic R term and thrust term from Task 8, and dividing by the mass from Task 3 to obtain v x 、v y 、v z The differential of the differential.

[0097] In this embodiment, Figure 5 The calculation of the ballistic height term is also involved, which does not affect the center of mass dynamics equation and the numerical solution process of ordinary differential equations.

[0098] S32: In the RK1 slave state, the data path control module transmits the current function value to the high-precision algorithm module. The adaptive step size control module configures the step size for the current algorithm stage. The high-precision algorithm module calculates the slope K2 based on the current step size and function value. The data path control module transmits the result parameter corresponding to the slope K2 to the function solution module. The function solution module uses the current result parameter as the new initial value to update the function value y n , RK6 slave state machine jumps from RK1 slave state to RK2 slave state;

[0099] In this embodiment, during the interactive operation between the function solving module and the high-precision algorithm module, the function solving module calculates the function value y using the initial input parameters or the result parameters fed back by the subsequent high-precision algorithm module. n , then the function value y nProvided to the high-precision algorithm module, the high-precision algorithm module converts the current time and function value y n Substitute the slope equation f(*) to calculate the slope, and then feed the function value corresponding to the current slope back to the function solving module. For example, in the RK1 slave state, the function solving module uses the initial input parameters to calculate the function value y n , and then transmit it to the high-precision algorithm module, which converts the function value y n Bring in , calculate K2, and then K2 corresponding Feedback to the function solving module, the function solving module will As the new initial value, calculate the new y n .

[0100] It should be noted that the adaptive step size control module needs to dynamically adjust the step size according to the divided algorithm stages to achieve the task of adjusting the step size at different stages. Because the RK6 algorithm uses a larger step size and the Adams6 algorithm uses a smaller step size, step size adjustment is prone to data instability and error propagation near the transition time node between the two algorithms, resulting in poor numerical solution stability of the ordinary differential equation. To address this problem near the transition node, the module design adopts a method of comparing the minimum stage step size to configure the step size.

[0101] Specifically, the adaptive step size control module configures the step size for the current algorithm stage, which includes the following process:

[0102] Calculate the current configuration phase time t n The absolute value of the difference from the dynamic real-time time t|t n -t|, compare the absolute value|t n -t| and configure the minimum stage step size h n-min The size of , if |t n -t|≥h n-min , indicating that the current configuration stage is not close to the transition node, then the configuration step size is the preset step size h corresponding to the current algorithm stage n , if |t n -t| <h n-min , indicating that the current configuration stage is approaching the transition node, the configuration step is adjusted to |t n -t|.

[0103] like Figure 7 As shown in the figure, the adaptive step size control module combines the reconfigurable method in hardware implementation, splitting the process of configuring the step size into one subtraction stage, one absolute value operation stage, and one comparison stage. Only one subtractor, one absolute value operator, and one comparator need to be called to complete the operation, so as to achieve the task of adjusting the step size applicable to different algorithm stages.

[0104] S33: In the RK2 slave state, the data path control module transmits the updated function value to the high-precision algorithm module, the adaptive step size control module configures the step size for the current new algorithm stage, and the high-precision algorithm module calculates the slope K3 based on the current step size and the updated function value. The data path control module transmits the result parameter corresponding to the slope K3 to the function solver module, and the function solver module updates the function value again. In this way, interactive operations are performed stage by stage until a total of 7 slope values ​​from K1 to K7 are obtained. The RK6 slave state machine jumps from the RK6 slave state to the RK7 slave state.

[0105] S34: In the RK7 slave state, the data path control module transmits the current function value to the high-precision algorithm module. The adaptive step size control module configures the step size for the current algorithm stage. The high-precision algorithm module calculates the initial solution y of the Adams6 algorithm based on the current step size, function value and various slope values. n+1 , the data path control module will initially decode n+1 Transmitted to the function solving module, the function solving module will be the initial solution y n+1 Update the function value as the new initial value;

[0106] S35: RK6 slave state machine jumps from RK7 slave state to RK1 slave state, repeats the process from S32 to S34 until the RK6 algorithm iteration is completed and a predetermined number of initial solutions are obtained. RK6 slave state machine jumps from RK7 slave state to idle state, wherein the initial solution includes y n-1 、y n-2 、y n-3 、y n-4 、y n-5 .

[0107] like Figure 3 As shown, during the iterative execution of the RK6 algorithm, when the 5 starting values ​​y are not met n-1 、y n-2 、y n-3 、y n-4 、y n-5 When the condition is met, RK6 starts from the state machine Figure 3 The "①" method iterates between the RK1-RK7 states. Conversely, when the conditions for providing 5 starting values ​​are met, Figure 3 The “②” method jumps to the idle state, ends the RK6 slave state machine jump task, and the main state machine transfers from RK6 to Adams6 state.

[0108] It should be noted that all intermediate calculation results generated in the above iterative process are stored in real time in the cache RAM inside the data cache module. The cached data includes K1, K2, K3, K4, K5, K6, K7, y n 、y n-1、y n-2 、y n-3 、y n-4 、y n-5 Etc., to provide continuous and accurate input data for the algorithm iteration module in the next stage.

[0109] like Figure 8 As shown in the figure, the data cache module includes 6 groups of RAMs for caching data information. Its storage reuse mechanism is based on the polling arbiter to arbitrate data storage and access. In the order of RAM1 to RAM6, the priority is refreshed once each iteration. During the iteration of the RK6 algorithm (i.e., the first 5 rounds of iteration), the write operation is allocated to the RAM with the highest priority according to the change of priority. The y obtained by the 5 iterations of the RK6 algorithm are respectively n-5 、y n-4 、y n-3 、y n-2 、y n-1 Stored in RAM1, RAM2, RAM3, RAM4, RAM5, after the Adams6 algorithm iteration starts (that is, after the sixth round of iteration starts), while assigning write operations to the highest priority RAM according to the priority change, the remaining groups of RAM are assigned read operations, so that the startup values ​​in RAM1 to RAM5 are read by the data path control module, and the calculated intermediate prediction values Stored in RAM6, in the next round of iteration, the write operation is assigned to RAM1 with the highest priority according to the change in priority, and the read operation is assigned to the remaining groups of RAM, so that the values ​​in RAM2 to RAM6 are read by the data path control module, and the calculated corrected numerical solution is stored in RAM1, and so on, to achieve the purpose of data refresh.

[0110] The storage reuse mechanism ensures that data remains valid after being written to the storage block by providing a life cycle for the data in the storage block until the data is invalidated and overwritten by new data, avoiding frequent transfers, reducing the logic complexity of the circuit, and improving resource utilization.

[0111] like Figure 8 As shown, since 5 groups of y are completed in sequence during the iteration of the RK6 algorithm n+1 The data of the value is stored and advanced according to the iteration stage when the Adams6 algorithm is iteratively run, so the data cache module requires 6 RAMs to complete the cache data information.

[0112] S4: After the iteration is completed, the main state machine jumps to the Adams6 main state. The data path control module starts to schedule the iteration process of the Adams6 algorithm. According to the jump state of the Adams6 slave state machine, the function solving module is controlled to execute the ordinary differential equation function value solving task stage by stage, and the solved function value is transmitted to the high-precision algorithm module. The high-precision algorithm module is controlled to iteratively calculate the final numerical solution of the ordinary differential equation. The data cache module caches all intermediate calculation results. The calculation formula of the Adams6 algorithm includes:

[0113]

[0114] Where, is the intermediate prediction value of the Adams prediction formula, which is used in the subsequent correction steps. is the final numerical solution after correction.

[0115] The datapath control module schedules the iteration process of the Adams6 algorithm to obtain the final numerical solution of the ordinary differential equation. Specifically, it includes the following steps:

[0116] S41: The data path control module transmits the initial solution cached in the data cache module to the high-precision algorithm module, and the Adams6 slave state machine jumps from the idle state to the Ad1 slave state;

[0117] S42: In the Ad1 slave state, the data path control module transmits the current function value to the high-precision algorithm module, the adaptive step size control module configures the step size for the current algorithm stage, and the high-precision algorithm module calculates the intermediate prediction value of the ordinary differential equation based on the obtained initial solution, current step size and function value And the intermediate prediction value Feedback to the function solving module, the function solving module will be the intermediate prediction value Update the function value y as the new initial value n , Adams6 slave state machine jumps from state Ad1 to state Ad2;

[0118] S43: In the Ad2 slave state, the data path control module transmits the current function value to the high-precision algorithm module, the adaptive step size control module configures the step size for the current algorithm stage, and the high-precision algorithm module calculates the value based on the obtained initial solution, current step size, function value and intermediate prediction value. Calculate the numerical solution of the ordinary differential equation and feed the current numerical solution back to the function solving module. The function solving module uses the current numerical solution as the new initial value to update the function value y n ;

[0119] S44: The Adams6 slave state machine jumps from state Ad2 to state Ad1, and repeats the process from S42 to S43 until the Adams6 algorithm iteration is completed and the final numerical solution of the ordinary differential equation is obtained. The Adams6 slave state machine jumps from the Ad2 slave state to the idle state.

[0120] like Figure 3 As shown in the figure, during the iteration of Adams6 algorithm, when the iteration condition is not met, Adams6 starts from the state machine and stops at Figure 3 "③" mode loops and iterates between the Ad1-Ad2 states. Conversely, when the iteration condition is met, Figure 3 The "④" method jumps to the idle state, ends the slave state machine task, and the master state machine is transferred from Adams6 to the idle state; among them, satisfying the iteration condition means reaching the preset number of iterations, and the number of iterations is set according to the trajectory expected flight time parameter in the initial input parameters. For example, when the trajectory expected flight time t = 1s and the step size h is 0.01, the number of iterations is set to 100 times.

[0121] like Figure 6 As shown in the figure, the high-precision algorithm module splits the operation process of the RK6-Adams6 hybrid algorithm into 7 multiplication stages and 5 addition stages, and realizes the parallel execution of tasks corresponding to different algorithm stages by reusing operators.

[0122] S5: After the iteration is completed, the operator state management module jumps to the idle state, the data path control module terminates the interactive task, and the data cache module uploads the final numerical solution of the ordinary differential equation to the host computer through the data transmission control module.

[0123] In this embodiment, when the ordinary differential equation solving task reaches a predetermined number of iterations, the iteration ends, the data cache module outputs the final numerical solution of the ordinary differential equation to the data transmission control module, and the data transmission control module uploads the final numerical solution of the ordinary differential equation to the host computer through the PCIe protocol.

[0124] In this embodiment, the solution circuit adopts a data-driven pipeline logic structure, which can split the calculation process into several consecutive small stages. Each small stage is executed using a different circuit module. For example, a single algorithm stage in RK6 is divided into small stages such as reading data, calculating slope, multiplying coefficients, weighted summation, and caching results. This allows multiple data to participate in the calculation at different stages at the same time, thereby improving calculation efficiency and reducing delays. The encoding format of the numerical data processed by the solution circuit is the IEEE754 standard 64-bit double-precision floating point number.

[0125] It should be noted that before executing this solution method, it is necessary to configure a data transmission control module, a data cache module, an operator state management module, an adaptive step size control module, a high-precision algorithm module, a function solving module and a data path control module on the FPGA-based lower computer hardware circuit. Among them, the function solving module pre-constructs the corresponding hardware operation path according to the form of the selected ordinary differential equation (such as the six-degree-of-freedom equation) so that it can smoothly execute the numerical solution process after receiving the input. In the high-precision algorithm module, a reconfigurable method is used based on the RK6-Adams6 hybrid algorithm to realize the design of efficient algorithm logic paths (that is, before running the high-precision algorithm module, a set of switchable reconfigurable logic path structures will be pre-designed, and then the operator state management module will dynamically schedule and configure path resources for each algorithm stage during runtime). The operation processes of different algorithms are uniformly planned, and efficient splitting and combination methods are used to dynamically reconstruct the logic paths in different algorithm stages (the logic paths include the seven-stage operation paths of the RK6 algorithm and the two-stage operation paths of prediction and correction of the Adams algorithm. The operation formulas of each stage are different, and the required paths are different). This meets the requirements of high-precision numerical solutions, improves resource utilization, and reduces the number of operators. The reconstructed logic path only requires 7 multipliers and 5 adders to complete the core computing tasks of all algorithms, significantly reducing circuit resource consumption.

[0126] The reconstructed logical path structure is as follows Figure 6 As shown, the logic path of the RK6 algorithm can cover the Adams6 algorithm. The algorithm logic path corresponds to the RK6 or Adams sub-state operation from the state machine. Figure 6 The variable values ​​of the input data ports of the RK6-Adams6 hybrid algorithm at different algorithm stages are shown in Table 1. Figure 6 The ports in the .

[0127] Table 1 Numerical values ​​of input variables for reconstructing the structure under different states

[0128]

[0129] The steps in this application can be adjusted in order, combined, and deleted according to actual needs.

[0130] The units in the device of the present application can be combined, divided and deleted according to actual needs.

[0131] Although the present application is disclosed in detail with reference to the accompanying drawings, it should be understood that these descriptions are merely exemplary and are not intended to limit the application of the present application. The scope of protection of the present application is defined by the appended claims and may include various modifications, alterations and equivalents made to the invention without departing from the scope and spirit of the present application.

Claims

1. A method for solving ordinary differential equations based on the RK6-Adams6 hybrid algorithm, characterized in that: The solution circuit includes a data transmission control module, a data cache module, an operator state management module, an adaptive step size control module, a high-precision algorithm module, a function solution module and a data path control module. The method includes: S1: The host computer transmits the initial input parameters of the ordinary differential equation to the data cache module through the data transmission control module, and the data cache module starts to cache the initial input parameters; S2: After the cache is completed, the main state machine of the operator state management module jumps to the RK6 main state, and the data cache module transmits the initial input parameters to the function solver module; S3: The data path control module schedules the iteration process of the RK6 algorithm. According to the jump state of the RK6 slave state machine, the function solving module and the high-precision algorithm module perform cyclic interactive iterative operations in each stage to obtain the initial solution of the Adams6 algorithm. Specifically: In a single iteration, the RK6 slave state machine jumps to the slave state in the order of RK1 to RK7. Each jump is an algorithm stage. In a single algorithm stage, the data path control module transmits the function value output by the function solving module to the high-precision algorithm module. The adaptive step size control module configures the corresponding step size according to the algorithm stage. The high-precision algorithm module calculates the slope corresponding to the function value based on the corresponding step size and the function value to obtain the initial solution, and feeds the calculated result parameters back to the function solving module. The function solving module updates the function value based on the current result parameters. This cycle iterates until a predetermined number of initial solutions are obtained. S4: After the iteration is completed, the main state machine jumps to the Adams6 main state. The data path control module schedules the iteration process of the Adams6 algorithm. According to the jump state of the Adams6 slave state machine, the function solving module and the high-precision algorithm module are controlled to perform cyclic interactive iterative operations stage by stage to obtain the final numerical solution of the ordinary differential equation. S5: After the iteration is completed, the operator state management module jumps to the idle state, the data path control module terminates the interactive task, and the data cache module uploads the final numerical solution of the ordinary differential equation to the host computer through the data transmission control module.

2. The method for solving an ordinary differential equation circuit based on the RK6-Adams6 hybrid algorithm according to claim 1, characterized in that: The arithmetic unit state management module also includes an RK6 slave state machine, and the S3 specifically includes: S31: The function solving module brings the initial input parameters into the ordinary differential equation and calculates the function value y of the ordinary differential equation n , RK6 slave state machine jumps from idle state to RK1 slave state; S32: The data path control module transmits the current function value to the high-precision algorithm module. The adaptive step size control module configures the step size for the current algorithm stage. The high-precision algorithm module calculates the slope K2 based on the current step size and the function value. The data path control module transmits the result parameter corresponding to the slope K2 to the function solving module. The function solving module uses the current result parameter as a new initial value to update the function value. The RK6 slave state machine jumps from the RK1 slave state to the RK2 slave state. S33: Perform interactive calculations stage by stage in the manner of S32 until a total of 7 slope values ​​from K1 to K7 are obtained, and RK6 jumps from the state machine to the RK7 slave state; S34: The data path control module transmits the current function value to the high-precision algorithm module, the adaptive step size control module configures the step size, and the high-precision algorithm module calculates the initial solution y of the Adams6 algorithm based on the current step size, function value, and each slope value. n+1 , the data path control module transmits the initial solution to the function solving module, and the function solving module uses the initial solution as a new initial value to update the function value; S35: The RK6 slave state machine jumps from the RK7 slave state to the RK1 slave state, and repeats the process from S32 to S34 until the RK6 algorithm iteration is completed and a predetermined number of initial solutions are obtained.

3. The method for solving an ordinary differential equation circuit based on the RK6-Adams6 hybrid algorithm according to claim 2, characterized in that: The calculation formula of the RK6 algorithm in S3 is expressed as: K1=f(t n ,y n ) Where K1 to K7 are the slope estimates of different algorithm stages within the time step of the RK6 algorithm, t n is the current time point, y n At time t n The function value at y n+1 is time t n+1 The function value at , h is the time step, and f(*) is the slope equation of time and function value.

4. The method for solving an ordinary differential equation circuit based on the RK6-Adams6 hybrid algorithm according to claim 2, wherein: In S32, the adaptive step size control module configures the step size for the current algorithm stage, specifically including: Calculate the current configuration phase time t n The absolute value of the difference from the dynamic real-time time t|t n -t|, compare|t n -t| and configure the minimum stage step size h n-min The size of , if |t n -t|≥h n-min , then the configuration step size is the preset step size h corresponding to the current algorithm stage n , if |t n -t| <h n-min , then the configuration step size is adjusted to |t n -t|.

5. The method for solving an ordinary differential equation circuit based on the RK6-Adams6 hybrid algorithm according to claim 1, wherein: The operator state management module also includes an Adams6 slave state machine, and the S4 specifically includes: S41: The data path control module transmits the initial solution cached in the data cache module to the high-precision algorithm module, and the Adams6 slave state machine jumps from the idle state to the Ad1 slave state; S42: The data path control module transmits the current function value to the high-precision algorithm module. The adaptive step size control module configures the step size. The high-precision algorithm module calculates the intermediate prediction value of the ordinary differential equation based on the obtained initial solution, the current step size, and the function value. and will Feedback to the function solving module, the function solving module will Update the function value y as the new initial value n , Adams6 slave state machine jumps from state Ad1 to state Ad2; S43: Interactive calculation is performed again in the manner of S42, based on the obtained initial solution, current step size, function value and Compute numerical solutions to ordinary differential equations; S44: The Adams6 slave state machine jumps from state Ad2 to state Ad1, and repeats the process from S42 to S43 until the Adams6 algorithm iteration is completed and the final numerical solution of the ordinary differential equation is obtained. The Adams6 slave state machine jumps from the Ad2 slave state to the idle state.

6. The method for solving an ordinary differential equation circuit based on the RK6-Adams6 hybrid algorithm according to claim 5, characterized in that: The calculation formula of the Adams6 algorithm in S4 is expressed as: Where, is the intermediate prediction value of the Adams prediction formula, which is used in the subsequent correction steps. is the final numerical solution after correction.

7. The method for solving an ordinary differential equation circuit based on the RK6-Adams6 hybrid algorithm according to claim 5, characterized in that: The initial input parameters include the expected flight time of the trajectory, the three components of the initial trajectory position, the three components of the initial trajectory velocity, the pressure-related parameters and the self-inertia force-related parameters. In S32, the result parameters are of the same type as the parameters included in the initial input parameters. The function solving module uses the result parameters fed back by the high-precision algorithm module as new initial values, and brings them into the ordinary differential equation to calculate the ordinary differential equation function value used in the next algorithm stage.

8. The method for solving ordinary differential equations based on the RK6-Adams6 hybrid algorithm according to claim 2, wherein: The function solving module includes multiple operation chains, each of which is used to calculate different parameters of ordinary differential equations. The function solving module takes the operation chain with the longest time as the critical path, merges the equation blocks that meet the merging conditions in the critical path, calculates other non-critical paths in parallel with the critical path, and performs multi-path sharing on the intermediate results that meet the sharing conditions in the operation chain to reduce the operation cycle and scheduling overhead. The merging conditions are that the input values ​​are the same, there is no dependency between each other, and the subsequent paths cross, and the sharing condition is that the intermediate result is used by the equation blocks in the subsequent link for more than or equal to 2 times.

Images