Dynamic system simulation and numerical solution method based on impulse response

Abstract

The invention relates to the field of computer simulation and scientific calculation, and provides a method for carrying out numerical solution on a dynamic system described by a generalized differential equation. According to the method, a system is divided into two subsystems which interact through a port function, a zero-input independent solution and an impulse response of the first subsystem are obtained in a computer, time domain convolution is carried out on the zero-input independent solution and the impulse response of the zero-input independent solution and the impulse response of the zero-input independent solution and the impulse response of the first subsystem and port input of a second subsystem, and a cross-subsystem effect is decomposed into an independent solution item, a historical effect item and an instant effect item; and forming an equivalent equation of the second subsystem and executing numerical integration to obtain a system time domain solution. The method is suitable for electrical, mechanical, hydraulic, control, structural, thermal, electromagnetic, fluid and other systems represented by generalized differential equations, and can be popularized to other time evolution models. According to the method, repeated solution of a linear part is replaced by impulse response, and the matrix operand is reduced while the spectral characteristics and precision are reserved, so that the simulation efficiency is remarkably improved, and the stability is improved.

Description

Technical Field

[0001] This invention relates to the fields of computer simulation and scientific computing technology, specifically to a method for dynamic system simulation and numerical solution based on impulse response. Background Technology

[0002] In modern engineering and scientific computing, the dynamic behavior of many physical systems can be described by time-domain equations. These equations originate from electrical systems, mechanical structures, hydraulic and pneumatic devices, thermal processes, fluid motion, chemical processes, and multi-physics coupled systems. After modeling, the time-domain equations of such systems can typically be expressed as ordinary differential equations or differential-algebraic equations. To perform dynamic simulations of these systems in a computer, numerical integration of these time-domain equations is required to obtain the system's response over time.

[0003] Based on the differences in system dimension and modeling granularity, simulation methods can be broadly classified into two categories: system-level one-dimensional simulation and field-level three-dimensional simulation. The former is characterized by lumped parameter modeling, focusing on describing the system's energy transfer, signal interaction, and control logic, and is commonly used in fields such as circuit and electronic systems, power and control systems, hydraulic servo systems, and chemical process systems. The latter typically discretizes the spatial continuum using methods such as finite element method, finite volume method, finite difference method, boundary element method, or spectral element method, and is used in applications such as structural dynamics, thermal processes, electromagnetic fields, and fluid simulation. Its discretization results also form a system of equations that evolve over time, which are solved in a computer using a numerical integrator.

[0004] In existing technologies, time-domain solutions for dynamic systems typically rely on direct numerical integration methods, such as explicit methods (e.g., Euler's method, central difference method) and implicit methods (e.g., trapezoidal rule, Newmark-β method, and Runge-Kutta family). Explicit methods have lower computational costs, but strict stability conditions, requiring extremely small step sizes for systems with high-frequency modes or rigid characteristics. Implicit methods, while superior in numerical stability, require solving a system of linear equations at each time step, causing computational costs to increase rapidly with system size. Traditional frameworks based on direct integration often struggle to simultaneously meet the combined requirements of efficiency, stability, and accuracy.

[0005] Besides direct integration methods, model order reduction methods are also common, including modal truncation, static condensation, and subspace projection. These methods accelerate simulation by reducing the system's degrees of freedom, but inevitably sacrifice high-frequency dynamic characteristics, resulting in incomplete spectrum, error accumulation over time, and potential introduction of non-physical energy exchange into coupled systems. For engineering applications that require maintaining complete dynamic characteristics, the applicability of order reduction methods is significantly limited.

[0006] In many engineering models, systems often contain both linear and nonlinear components. Existing numerical integration methods require repeated solutions to the linear component at each time step, even if its characteristics remain unchanged during simulation. This repetitive computation not only leads to efficiency losses but can also introduce accumulated errors due to multiple approximations of high-frequency components, reducing the stability of long-term integration. The computational burden increases further as system size, coupling degree, and time-varying parameter characteristics increase. Typical scenarios include electrical and electronic systems, structural dynamics, electromagnetic transients, thermo-fluid-structure interaction, fluid motion, and multibody systems. These systems commonly exhibit cross-port, cross-domain, or cross-module interactions, typically requiring synchronous solutions within a unified framework. This makes it difficult to fully utilize the pre-calculated linear structure and achieve efficient decomposition or reuse between modules. Existing time-step-based integration architectures repeatedly process the spectral characteristics of the linear component at each step, resulting in significant repetitive computation and limiting overall solution efficiency.

[0007] In summary, existing technologies generally suffer from problems such as repeatedly solving linear parts, inability to reuse structural information, difficulty in maintaining spectral integrity, and difficulty in achieving efficient large-scale simulation. These problems become increasingly apparent as the model size and coupling complexity increase. Summary of the Invention

[0008] This invention aims to address the common technical problems in existing dynamic system simulation methods, such as low computational efficiency, high overhead due to repeated solutions of linear parts, and difficulty in achieving a balance between stability and accuracy. It proposes a numerical solution method for dynamic systems based on impulse response, executed in a computer. The method divides the system to be solved into two subsystems coupled through port functions. The independent solution of the first subsystem under zero-input conditions and its impulse response to port inputs are obtained in the computer. The impulse response is then convolved with the port inputs of the second subsystem in the time domain to analytically decompose the cross-subsystem action into independent solution terms, historical action terms, and immediate action terms. Based on this, the computer constructs an equivalent equation containing these three terms and performs numerical integration on the equivalent equation, thereby avoiding repeated solutions of the linear subsystem during the time-domain solution process.

[0009] The method of this invention is applicable to various dynamic systems described by ordinary differential equations, differential algebraic equations, or time-continuous equations formed by spatial discretization. These systems can originate from circuit networks, hydraulic systems, mechanical or multibody systems, structural or thermal processes, electromagnetic fields, acoustics, fluid mechanics, control systems, or comprehensive models formed by the coupling of multiple physical domains.

[0010] In this invention, the impulse response of the first subsystem can be obtained in a computer through various means.

[0011] Firstly, an approximate impulse-type port input can be applied to the first subsystem in the time domain, and the corresponding port output can be recorded during the numerical integration process, thereby directly constructing the time-domain impulse response.

[0012] Secondly, based on the frequency response characterizing the coupling port relationship of the first subsystem, the time-domain impulse response can be obtained in the computer through inverse Fourier transform, and if necessary, causality correction can be performed on the frequency response or the time-domain impulse response generated therefrom to make the calculation results satisfy physical causality.

[0013] Third, modal decomposition can be performed on the first subsystem to obtain modal parameters in the computer to characterize its linear dynamic characteristics, and the impulse response can be reconstructed in the time domain based on the selected dominant mode. In this process, one or more residual vectors or equivalent correction vectors can be combined to compensate for the dynamic characteristics that did not participate in the modal reconstruction and improve the reconstruction accuracy.

[0014] The above-mentioned approaches can be applied to different types of system architectures and data sources, and provide high modeling flexibility and numerical stability while maintaining physical consistency.

[0015] The convolution decomposition mechanism of this invention transforms cross-subsystem coupled computation from solving a global matrix of size N into convolution operations at the port dimension, with a single-step computational complexity of approximately O(P²·L), where P is the number of ports and L is the effective length of the impulse response. In engineering systems, P is typically much smaller than N, thus this invention significantly reduces the computational load of simulations while preserving the system's spectral characteristics. Furthermore, to accommodate different time integration strategies, the computer can perform truncation, windowing, resampling, or interpolation processing on the impulse response sequence before convolution to suppress high-frequency ringing, reduce numerical noise, and ensure energy consistency.

[0016] Compared with traditional direct numerical integration methods, this invention preserves the inherent spectral characteristics of the system by obtaining the impulse response in advance, so that the time-domain solution maintains high consistency in the low, medium and high frequency ranges. Unlike existing network equivalence or model order reduction methods, this invention can significantly reduce the amount of computation without omitting dynamic characteristics or introducing approximate assumptions. Therefore, it still has good stability and scalability in rigid systems, large-scale models and multi-physics domain coupled scenarios.

[0017] In summary, this invention achieves an effective balance between efficiency, stability, and accuracy, providing a unified time-domain solution mechanism for various dynamic systems, and has significant engineering practical value and a wide range of applications. Attached Figure Description

[0018] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the contents shown in the accompanying drawings are briefly described below. It should be understood that the following drawings are only used to illustrate some embodiments of the present invention. For those skilled in the art, other equivalent structures or processes can be obtained based on these drawings without creative effort: Figure 1 This is a flowchart of the numerical solution method for dynamic systems according to an embodiment of the present invention; Figure 2 This is a structural diagram of the software functional modules used to implement the method according to an embodiment of the present invention; Figure 3 This is a schematic diagram of a computer device structure for implementing the method according to an embodiment of the present invention. Explanation of reference numerals in the attached figures

[0019] In the accompanying drawings of this invention, the same or similar reference numerals denote the same or similar components or steps. The following lists the main reference numerals and their corresponding component names to illustrate the structure and execution flow of this invention: 100 — Overall Flowchart of Numerical Solution Methods for Dynamic Systems 101 — System Partitioning Steps 102 — Steps for calculating the characteristics of a linear subsystem 103 — Steps for constructing equivalent equations 104 — Steps for Time-Domain Numerical Solution and Result Output 200 — Overall Structure of Functional Modules in Numerical Solving Software for Dynamic Systems 201 — System Module Division 202 — Linear Subsystem Characteristic Calculation Module 203 — Equivalent Equation Construction Module 204 — Time-Domain Numerical Solution and Result Output Module 300 — Overall Structure of Computer Device 301 — Processor 302 — System Bus 303 — Memory 304 — Storage Media 305 — Program Module Detailed Implementation

[0020] The embodiments of this invention are merely illustrative of the principles and implementation methods of the invention, and are not intended to limit the scope of protection. For those skilled in the art, any equivalent substitutions or modifications made to the following embodiments without departing from the spirit and essence of this invention should be considered to fall within the scope of protection of this invention.

[0021] The method of this invention is executed in a computer for numerical integration and simulation calculation of equations in a dynamic system, addressing the problems of insufficient simulation efficiency, limited numerical stability, and difficulty in achieving both high and low accuracy commonly found in existing technologies. This method belongs to the field of computer simulation and scientific computing technology and can be implemented as a software module, function library, or interface component in general simulation software, circuit and control simulation platforms, finite element analysis programs, and computational fluid dynamics (CFD) solvers, thereby embedding it into engineering-level simulation systems to perform time-domain numerical solutions.

[0022] This invention is applicable to various dynamic systems described by time-domain equations. Such systems, after modeling, can typically be represented as ordinary differential equations (ODEs) or differential-algebraic equations (DAEs). For models obtained by discretizing a spatial continuum using methods such as finite element method, finite volume method, finite difference method, boundary element method, or spectral element method, their semi-discrete equations can also be equivalently represented as ODEs or DAEs in time. Regardless of whether the system originates from electrical and electronic, control, mechanical, multibody, hydraulic and pneumatic, thermal, chemical processes, structural vibration, electromagnetic fields, acoustics, fluid mechanics, or a comprehensive model formed by the coupling of multiple physical domains, after parametric modeling or spatial discretization, its time-domain behavior can ultimately be uniformly represented as a set of equations that evolve over time.

[0023] In dynamic system simulation, modeling methods can be categorized into one-dimensional simulation and three-dimensional simulation based on whether the state variables depend on spatial coordinates. One-dimensional simulation uses lumped parameter modeling, where each component constructs the overall system through port quantities and algebraic constraints. Three-dimensional simulation uses spatial discretization to form semi-discrete equations, used for continuous medium analysis such as structural, thermal, electromagnetic, and fluid dynamics. Despite the different modeling approaches, the time-domain expressions of these models can ultimately be reduced to ODEs or DAEs. In this invention, this unified set of time-continuous equations is called generalized differential equations, and this method is based on this unified mathematical representation to achieve numerical solutions for dynamic systems in a computer.

[0024] In one implementation, to facilitate the description of the internal structure of the overall dynamic system, the system can be divided into several relatively independent parts, each called a subsystem. Subsystems are calculated using their own models, while different subsystems still interact with each other. These interactions are typically transmitted through port quantities. Port quantities can be pairs of related quantities representing power exchange, such as force and displacement, voltage and current, or pressure and flow rate, used to reflect the interaction between subsystems.

[0025] To express such cross-subsystem interactions within a unified structure, this invention abstracts the influence relationships corresponding to port quantities into port functions, enabling the interactions between different subsystems to be described within the same mathematical framework. Based on this representation, the generalized differential equation of the dynamic system can be written as: f(x, dx / dt, d²x / dt², …, t) = p(y, t) g(y, dy / dt, d²y / dt², …, t) = q(x, t) Where x and y represent the state variables of the two subsystems respectively, f(·) and g(·) are the generalized operators of each subsystem; p(·) and q(·) are port functions used to characterize the effect of the other subsystem on this subsystem, and can be linear, nonlinear or time-varying.

[0026] In one implementation, the generalized differential equations denoted by f(·) and g(·) can serve as a unified mathematical form describing various types of dynamic systems. This form can represent systems with algebraic constraints or higher-order derivative terms, and can also naturally degenerate into ordinary differential equations without these terms, such as dx / dt = f(x, t) + p(y, t) and dy / dt = g(y, t) + q(x, t). For the case of using a linear time-invariant model, this generalized form can also be equivalently written as a state-space structure, such as dx / dt = A_x·x + P·y and dy / dt = A_y·y + Q·x. Furthermore, the semi-discrete equations obtained through spatial discretization can also be represented in the time domain as similar ODE or DAE structures, and therefore also fall under this generalized category. Based on this unified expression, this invention can perform convolutional solutions and construct equivalent equations for dynamic systems from different sources and with different modeling methods within the same solution framework, thereby achieving the method's universal applicability to various types of systems.

[0027] After establishing the aforementioned unified mathematical representation, the method of this invention can be executed sequentially in a computer according to a predetermined process to complete the numerical solution of the dynamic system. For example... Figure 1 As shown, the dynamic system to be solved is first structurally partitioned (step 101), dividing the system into two interacting subsystems and determining their port relationships. Then, the linear dynamic characteristics of the first subsystem are extracted in the computer (step 102), including obtaining its free response under zero-input conditions and its time-domain transfer characteristics to port excitations. Based on these linear characteristics, an equivalent time-domain expression to replace the cross-subsystem interaction can be constructed in the equations of the second subsystem (step 103), forming a numerical solution structure that can be directly time-progressed. Afterwards, numerical integration is performed on the equivalent time-domain equation, generating the system's time-domain response results (step 104). All of the above steps can be automatically completed by a computer program, and data interaction is achieved through a unified data structure and interface, thus forming a clear, reusable, and universally applicable solution process. The following combines... Figure 1 The steps of the method of the present invention will be further explained below.

[0028] (Step 101) System Partitioning In this step, the overall system is divided into a first subsystem and a second subsystem based on its structural characteristics, and the interaction between the two is established through port quantities. This division method allows the influence across subsystems to be handled in a unified form, thus providing clear computational boundaries for subsequent steps such as impulse response calculation and convolution decomposition.

[0029] Within the aforementioned unified framework, the coupling relationship between the two subsystems can be summarized as follows: f(x, dx / dt, d²x / dt², …, t) = p(y, t) g(y, dy / dt, d²y / dt², …, t) = q(x, t) Here, x and y represent the state variables of the two subsystems, f(·) and g(·) are generalized operators describing their respective dynamic characteristics, and p(·) and q(·) are port functions used to abstract the port inputs applied by the other subsystem.

[0030] In a preferred partitioning method, the first subsystem maintains linear constant-coefficient characteristics and does not contain nonlinear operators that vary with state or time, so that its impulse response is mathematically unique and numerically reusable. The second subsystem can contain nonlinear, time-varying, or discrete event characteristics. For example, in mechanical or structural systems, nonlinear materials, friction, contact, or gaps can be classified into the second subsystem; in electrical systems, nonlinear devices such as diodes, transistors, and saturated inductors can be classified into the second subsystem; in hydraulic systems, valve flow equations, accumulator gas compression processes, and friction and cavitation effects can also preferably be classified into the second subsystem, thereby preserving the nonlinear behavior of the overall system while ensuring the linear time-invariant characteristics of the first subsystem.

[0031] In some implementations, the equations of the original dynamic system do not directly represent the port input structure adopted in this invention. For example, coupling terms caused by higher-order derivatives, inertial coupling, damping coupling, or other cross-domain operators may not appear in the form of first-order terms. To enable the system to be incorporated into the unified solution framework of this invention, the computer can perform a formal transformation on the original equations, explicitly organizing the cross-subsystem actions into port input terms. This transformation can be achieved through explicit processing, state expansion, or equation rearrangement. For example, the first derivative of displacement, velocity, current, and flux linkage derivative can be introduced as new state variables, the higher-order equations can be expanded into a first-order state-space form, and the cross-subsystem actions can be uniformly merged into the input terms corresponding to the port functions. After the above processing, the first subsystem receives the action from the second subsystem only through the port function p(y, t), and the second subsystem receives the action from the first subsystem only through q(x, t), thereby transforming the originally complex coupling relationship into a clear and stable first-order port coupling structure within the framework of this invention.

[0032] In a preferred embodiment, by appropriately dividing the subsystem, the port input of the second subsystem is mathematically represented as a linear linear term, i.e.: g(y, dy / dt, d²y / dt², …, t) = Q · x Here, Q is the port matrix, which uses a linear relationship with constant coefficients to describe the port effect exerted by the first subsystem on the second subsystem. This linear structure allows the impulse response of the first subsystem to be directly incorporated into the equations of the second subsystem through time-domain convolution, thus providing a clear and stable mathematical foundation for subsequent convolution solutions and the construction of equivalent equations.

[0033] The interaction across subsystems is transmitted through port quantities, which can consist of a pair of related quantities reflecting energy or signal exchange, such as force and displacement, voltage and current, or pressure and flow rate. The interface relationships established through port quantities can represent the coupling modes of different models under a unified structure. Whether derived from a system-level model or a discretized continuous medium model, they can all be abstracted into port function forms. The specific structure of the port function depends on the selection of interface variables and can be linear, nonlinear, or time-varying. In a preferred embodiment, the port input of the second subsystem uses a linear first-order term Q·x to obtain a clear and directly processable coupling form.

[0034] Through the above division and port definition, dynamic systems from different fields such as electrical, mechanical, hydraulic, thermal, electromagnetic, and acoustic can be included in a unified solution framework and models can be constructed in a modular manner, providing a consistent and stable foundation for subsequent impulse response calculation, convolution solution, and equivalent equation construction.

[0035] (Step 102) Calculation of characteristics of linear subsystem In this step, the computer calculates the discrete free response sequence of the first subsystem under zero input conditions and the discrete impulse response sequence under unit impulse input at the port, thereby obtaining the independent solution x0[n] and the impulse response h[n]. The time domain is discretized into a series of equally spaced time nodes n·Δt, where Δt is the simulation time step and n is a non-negative integer used to identify the value taken by the computer at the nth time node. The independent solution x0[n] represents the free response of the first subsystem under the initial conditions when the external input is zero at time n·Δt; the impulse response h[n] represents the discrete response of the first subsystem to unit impulse excitation at time n·Δt, used to characterize the dynamic transfer characteristics of the port input on the discrete time grid. The above sequences are all stored in the computer in the form of discrete arrays, which can be directly used for subsequent equivalent equation construction and time step integration.

[0036] The two types of responses mentioned above, as pre-calculated linear characteristic data of the first subsystem, will be directly introduced in the subsequent construction of equivalent equations as independent solution terms and convolution terms, thereby avoiding repeated processing of the linear part and improving the efficiency and stability of time-domain simulation.

[0037] In the implementation of finding the independent solution x0[n], the computer can select different numerical strategies according to the model structure of the first subsystem.

[0038] In one implementation, the computer can use an analytical method to obtain independent solutions. When the first subsystem can be represented in the first-order state-space form dx / dt = A·x + B·y, under the condition that y(t) = 0, its continuous domain free response x0(t) can be obtained first, and the independent solution sequence x0[n] can be obtained by sampling at the discrete time node n·Δt. For example, its analytical expression can be written as x0(t) = exp(A·t) · x0(0), and the corresponding discrete sequence can be obtained by taking values ​​at each n·Δt. For higher-order equations or models containing algebraic constraints, the computer can construct its free response through state extension, equivalent transformation of equations, or generalized inverse operators, and similarly form the independent solution sequence x0[n] on the discrete-time grid.

[0039] In another implementation, the computer can use time-domain numerical integration to obtain independent solutions. Given initial conditions, integrating the linear equation step-by-step yields a discrete-time series of independent solutions. Integration methods can include Euler's method, trapezoidal rule, Runge-Kutta method, or generalized-α method, which are suitable for non-diagonalizable matrix structures or large-scale sparse systems.

[0040] In another implementation, the computer can use a frequency domain method to obtain independent solutions. When the first subsystem exhibits periodic oscillations, weakly damped decay, or dynamic behavior that can be described by frequency components under zero-input conditions, frequency analysis can be performed on its linear operators to obtain the main frequency components constituting the independent solutions. The frequency domain analysis can be achieved through Fourier analysis, harmonic balancing, or characteristic frequency extraction. Based on the obtained frequency components, the computer can construct discrete-time independent solutions in two ways: on the one hand, the harmonic superposition value of the independent solutions can be directly calculated at the discrete-time node n·Δt based on the selected frequency components, thereby forming an independent solution sequence x0[n]; on the other hand, an inverse Fourier transform can be performed on the determined frequency components to generate the corresponding time-domain sequence, and x0[n] can be obtained by sampling on the discrete-time grid.

[0041] In a preferred embodiment, the linear subsystem characteristic calculation module can automatically select the strategy for obtaining independent solutions based on the model structure and dynamic characteristics of the first subsystem. When the first subsystem has a linear constant coefficient structure and its free response can be directly obtained through analytical expressions, the program can preferentially use analytical methods to calculate the free response in the continuous domain and sample at discrete time nodes n·Δt to form an independent solution sequence x0[n]. For models whose free responses can converge within a short time range, the program can use numerical integration to perform short-time integration on the linear equations under given initial conditions and generate x0[n] by discrete sampling. For models exhibiting periodic oscillations or low-damping decay characteristics, the program can use frequency domain methods to obtain the main frequency components through Fourier analysis, harmonic balance, or characteristic frequency extraction, and generate time-domain sequences by directly calculating the harmonic superposition values ​​at discrete nodes or performing inverse Fourier transform to construct independent solutions x0[n]. Through the above strategy selection mechanism, the independent solutions of the linear subsystem can be obtained in an efficient and model-adaptive manner during the program initialization phase, providing unified and reusable discrete sequence data for subsequent equivalent equation construction.

[0042] In the method of this invention, the impulse response of the first subsystem is used to characterize its time evolution characteristics when subjected to a unit impulse input at a port. For the single-port case, the impulse response can be represented as a sequence h[n] that varies with discrete time step n. For the multi-port case, since a unit impulse at each input port will produce a corresponding response at all output ports, the values ​​of the impulse response at each time step constitute a matrix between output ports and input ports, which can be written as h[n][i,j], where i represents the output port and j represents the input port. As time progresses, these matrices are arranged in a time sequence to form a three-dimensional response structure with time dimension, output port dimension, and input port dimension, thereby completely representing the dynamic transmission relationship of the multi-port linear subsystem. The specific storage method of this structure is not limited and can be represented as a three-dimensional array, a matrix sequence, or other equivalent data organization form. As long as it can reflect the response relationship between ports at discrete time steps and support the calculation of independent solution terms, historical action terms, and immediate action terms in subsequent convolution solutions, it should be considered to fall within the scope of the method of this invention.

[0043] To accommodate different system architectures and modeling requirements, impulse responses can be obtained in various ways. Three typical implementation methods are given below.

[0044] In one implementation, the computer uses a time-domain method to obtain the impulse response sequence h[n] of the first subsystem. The core idea of ​​this method is to apply a perturbation in the form of a unit pulse at the port and record the response of each port at discrete time nodes, thereby obtaining the dynamic transmission characteristics of the system to instantaneous excitation.

[0045] In the single-port case, the computer applies a unit impulse input to the port at time step n=0. In the discrete solution framework, this impulse can be directly represented as a discrete unit impulse sequence that takes the value 1 at n=0 and 0 at n>0. If the first subsystem is solved in the form of a continuous differential equation and through a numerical integrator, the program can apply a narrow pulse signal with area normalization within a very short time interval of t≈0 to approximate a unit impulse in a continuous sense. Subsequently, the computer records the response value of the port at each discrete time node n·Δt and generates a one-dimensional impulse response sequence h[n] accordingly.

[0046] In the multi-port scenario, the program executes the above process for each port separately. For the j-th port, the program applies a unit impulse perturbation to its port variable and records the response values ​​of all ports on the same time grid n·Δt. Thus, at each time step k, the system's response to the j-th port can be represented as a port mapping matrix h[n][i, j], where i is the index of the responding port and j is the index of the port to which the impulse was applied. As time progresses, {h[n][i, j]} forms a set of arrays with a uniform structure in the time dimension, which can completely describe the dynamic transmission relationship of the multi-port linear subsystem.

[0047] In numerical implementation, the program automatically selects either a discrete unit pulse or a narrow pulse approximation based on the solution method to ensure that the input disturbance has a unit area without introducing additional high-frequency errors. For models using continuous integrators, the program can continue to advance several steps after the pulse ends to capture the free decay process, ensuring that the impulse response sequence accurately reflects the dynamic characteristics of the system over a sufficiently long time range.

[0048] In another implementation, the computer uses a frequency domain method to obtain the impulse response sequence h[n] of the first subsystem. This method obtains the port frequency response of the system at a set of discrete frequency points and performs an inverse numerical transformation on the frequency domain data to obtain the corresponding time-domain impulse response. The program first establishes frequency sampling points within a set frequency range and obtains the complex frequency response between ports at each sampling frequency. The frequency response can come from the calculation results of the system's frequency domain equations, data from actual measurements, or amplitude and phase information provided by external files; in some implementations, the frequency response can be simplified or reduced in order before use, but its source still corresponds to the three types of approaches mentioned above.

[0049] Before performing the inverse transform, the program can perform sampling consistency checks and format adaptation on the frequency response data to ensure that the frequency point spacing, data arrangement, and complex number format meet the requirements of the inverse transform algorithm. After completing the above preprocessing, the frequency response data can be used for numerical inverse transform.

[0050] During the inverse transform, the program selects the appropriate data organization method based on the number of ports. For a single-port system, the frequency response varies only with frequency and can be represented as a one-dimensional array H(ω). The computer directly performs a discrete Fourier inverse transform on this sequence to obtain the time-domain impulse response h[n] on a unified time grid. For a system with multiple coupled ports, the frequency response varies with three dimensions: frequency, input port, and output port. It can be organized into a three-dimensional structure H(ω)[i, j], where j is the applied port and i is the response port. In the numerical implementation, the computer performs an element-wise IFFT along the frequency dimension, that is, it independently performs an inverse transform on the frequency sequence H(ω)[i, j] of each pair of ports (i, j) and organizes the results into h[n][i, j] according to the time index. The resulting multi-input multi-output impulse response array can be directly used for subsequent convolution solutions and equivalent equation construction.

[0051] In practical applications, since the frequency response is typically sampled at only a finite number of frequency points, such as several equally spaced points within the range of 0 Hz to 10 kHz, the resulting data belongs to the band-limited spectrum. When performing an inverse transform directly on this type of data, the time-domain results may exhibit fluctuations and ringing phenomena; this error is known as the spectral truncation effect or the Gibbs phenomenon. Furthermore, when the frequency response lacks certain spectral components, has an uneven sampling distribution, or exhibits asymmetry in the complex plane, the resulting time-domain response may show non-zero values ​​in the negative time interval, i.e., it exhibits a non-causal response. A non-causal response means that the system produces an output before the excitation, which is physically impossible and can lead to numerical instability or response distortion in time-domain simulations.

[0052] To overcome the aforementioned problems, this invention introduces a causality correction mechanism during the inverse frequency transform process. This mechanism, based on the analytical properties and causal constraints of linear time-invariant systems, restores or compensates for missing spectral components in the frequency domain, ensuring that the time-domain response after the inverse transform physically satisfies the causality condition. Through this process, the computer can effectively eliminate non-causal components caused by finite sampling, incomplete spectra, or uneven frequency distribution, obtaining a stable and physically consistent time-domain impulse response.

[0053] The correction process can be achieved using various mathematical relationships, including: the Hilbert transform relationship between the real and imaginary parts of the frequency response, the Hilbert transform relationship between the logarithmic magnitude and phase, spectral symmetry, or minimum phase reconstruction. These methods can be used individually or in combination. Their theoretical basis includes the Kramers-Kronig relation for causal systems and the principle of spectral symmetry in analytic signal theory.

[0054] The impulse response after causality correction accurately reflects the transient characteristics of the first subsystem to the port input in the time domain, and can serve as the basic data input for constructing subsequent equivalent equations. In different implementations, the computer can automatically select an appropriate causality correction method based on the type and completeness of the acquired frequency response data to ensure that the time-domain response obtained by the inverse transformation satisfies the causality condition, i.e., the system does not produce a response before the excitation. All of the above-mentioned correction methods fall within the protection scope of this invention.

[0055] In one correction method (correction method 1), the computer recovers the complete frequency response by reconstructing or correcting the imaginary part from the real part. When the imaginary part of the frequency response is missing, interfered with by noise, or does not satisfy the Hilbert duality relation, the program performs a Hilbert transform on the real part to obtain or correct the corresponding imaginary part, specifically... Im{H(ω)} = -Hilbert[Re{H(ω)}], Here, Hilbert represents the Hilbert transform operator. The complete complex frequency response can be constructed using this relationship. H(ω) = Re{H(ω)} + j·Im{H(ω)}.

[0056] The computer then performs a discrete inverse Fourier transform on the obtained H(ω) to obtain the causal time-domain impulse response h[n]. This method ensures that the real and imaginary parts satisfy the Kramers–Kronig relation through frequency domain reciprocity, thereby recovering the phase characteristics consistent with the linear time-invariant system.

[0057] In another correction method (correction method 2), the computer directly processes the time-domain sequence after the inverse real-part transformation using a weighted truncation method. The specific steps are as follows: First, Re{H(ω)} is extracted from the frequency response, and a discrete inverse transform is performed on its real part to obtain the initial time-domain sequence x[n]. Subsequently, the computer weights and amplifies the positive time portion and sets the negative time portion to zero to achieve energy compensation and causal constraints. This operation can be represented as follows: h[n] = 2·x[n], when n>0; h[0] = x[0]; h[n] = 0 when n < 0.

[0058] By multiplying the positive time component by 2, the system maintains a constant overall power after removing the negative time energy, thus ensuring amplitude consistency. The processed h[n] is the impulse response that satisfies causality.

[0059] In another correction method (correction method 3), the computer employs a symmetric point superposition method. Let x[n] be the original time-domain sequence obtained from the frequency-domain samples through the inverse discrete Fourier transform. This method eliminates negative time energy and strengthens the positive time response by superimposing and folding the positive and negative time symmetric points in the time-domain sequence, thereby ensuring the causality of the system's impulse response. The calculation process is as follows: h[n] = x[n] + x[Nn], when 1 ≤ n ≤ N / 2-1; h[0] = x[0]; h[n] = 0 when n ≥ N / 2.

[0060] This method folds the energy originally distributed in the negative time to the positive time by superimposing the values ​​of the positive and negative time symmetry points, thereby keeping the total energy and DC component of the system unchanged, while ensuring that the output is zero in the negative time period.

[0061] In another correction method (correction method 4), the computer generates the complete frequency response by reconstructing or correcting the real part from the imaginary part. When the real part of the frequency response is missing, affected by noise, or does not satisfy the Hilbert duality relation, the program performs a Hilbert transform on its imaginary part to obtain or correct the corresponding real part, specifically... Re{H(ω)} = Hilbert[Im{H(ω)}], Here, Hilbert represents the Hilbert transform operator. The complete complex frequency response can be reconstructed from this. H(ω) = Re{H(ω)} + j·Im{H(ω)}.

[0062] Subsequently, the computer performs a discrete inverse Fourier transform on the obtained H(ω) to obtain the time-domain impulse response h[n] that conforms to causality.

[0063] In the numerical implementation, the program can also correct the value of the impulse response at time zero to eliminate the offset caused by spectral truncation. This correction can be calculated using the following formula: h[0] = (1 / N) · ∑_k H(ω_k), Where N is the number of frequency sampling points, and ∑_k represents the summation over all frequency points.

[0064] This method is suitable for situations where the imaginary part information is relatively complete but the real part data is missing or has high noise. It can recover the amplitude-phase consistent spectrum and generate a stable, physically reasonable causal response. Based on the Kramers–Kronig relation and the principle of spectral reciprocity, this method ensures that the recovered frequency response satisfies analytical consistency between the real and imaginary parts, thereby eliminating non-causal components and maintaining the physical consistency of the system.

[0065] In another correction method (correction method 5), the computer employs a weighted truncation method based on the imaginary part. This method takes the imaginary part Im{H(ω)} of the frequency response as input, first performing an inverse discrete Fourier transform on it to obtain the initial time-domain sequence x[n]. Subsequently, the program achieves causal correction through weighted and truncation operations. The specific calculation formula is as follows: h[n] = 2·x[n], when n > 0; h[n] = 0 when n ≤ 0; h[0] = (1 / N) · ∑_k H(ω_k), Where N is the number of frequency sampling points, and ∑_k represents the summation over all frequency points.

[0066] This processing method eliminates the negative time component while maintaining the consistency of amplitude and energy, and recovers the DC component by averaging the zeros. This method is suitable for situations where the imaginary part is the primary input or where the imaginary part information is relatively reliable, and can effectively reconstruct a causal time-domain impulse response even when frequency domain information is incomplete.

[0067] In another correction method (correction method 6), the computer employs the symmetric point difference method. Let x[n] be the original time-domain sequence obtained from the frequency-domain samples through the inverse discrete Fourier transform. This method calculates the difference between the positive and negative time symmetric points in the time-domain sequence to remove even-symmetric components and retain odd-symmetric components, thereby ensuring the unidirectionality of the system's impulse response. The calculation process is as follows: h[n] = x[n] - x[Nn], when 1 ≤ n ≤ N / 2-1; h[0] = x[0]; h[n] = 0 when n ≥ N / 2.

[0068] This method complements correction method 3 and is suitable for situations where the spectrum contains strong even-symmetric components. It can effectively suppress spurious pre-response caused by symmetric truncation.

[0069] In another correction method (correction method 7), the computer reconstructs the phase using logarithmic amplitude to achieve minimum phase causal correction. Let |H(ω)| represent the amplitude of the frequency response. The program first calculates its logarithmic amplitude ln|H(ω)|, and then performs a Hilbert transform on the result to generate the corresponding phase: φ(ω) = - Hilbert[ln|H(ω)|].

[0070] Then the amplitude and phase are combined as H(ω) = |H(ω)| · exp(j·φ(ω)), An inverse transform is then performed on the system to obtain the causal minimum-phase time-domain response h[n]. This method is particularly effective when the system has minimum-phase characteristics (such as no zeros in the right half-plane or upper half-band), and the resulting response is stable and convergent in the time domain.

[0071] In another correction method (correction method 8), the computer reconstructs the logarithmic magnitude by phase reconstruction to recover the complete frequency response. Let φ(ω) represent the phase of the frequency response, and the program performs a Hilbert transform on the phase to generate the corresponding logarithmic magnitude: ln|H(ω)| = Hilbert[φ(ω)], And from this, the amplitude can be obtained: |H(ω)| = exp(ln|H(ω)|).

[0072] The computer then combines the reconstructed amplitude with the original phase to form a complex frequency response: H(ω) = |H(ω)| · exp(j·φ(ω)).

[0073] By performing an inverse discrete Fourier transform on the frequency response, the time-domain impulse response h[n] can be obtained.

[0074] To maintain amplitude normality, the program can correct the value of the impulse response at time zero in the final step: h[0] = (1 / N) · ∑_k H(ω_k), Where N is the number of frequency sampling points, and ∑_k represents the summation over all frequency points.

[0075] Based on the analytical characteristics of minimum phase systems and the principle of amplitude complementarity, this method can recover a stable and physically consistent time-domain impulse response when the phase information is reliable but the amplitude data is incomplete.

[0076] From a mathematical perspective, the eight correction methods mentioned above can be classified into two main categories: correction methods 1 to 6 are based on the analytical relationship between the real and imaginary parts of the frequency response; correction methods 7 and 8 are based on the analytical relationship between the logarithmic magnitude and the phase.

[0077] These two types of relationships correspond to two equivalent representations of causal systems in the complex frequency domain. The former reflects the even and odd symmetry constraints of the frequency response in real and imaginary spaces, while the latter reflects the Hilbert transform coupling law between amplitude and phase.

[0078] In correction methods 1 to 6, the real and imaginary parts of the system follow the Kramers–Kronig relation: If the system is linear, time-invariant, and causal, then its frequency response H(ω) can be written as: H(ω) = Re{H(ω)} + j·Im{H(ω)}, And the two parts satisfy the following mutual relationship: Im{H(ω)} = - Hilbert[Re{H(ω)}], Re{H(ω)} = Hilbert[Im{H(ω)}].

[0079] This means that the spectral information of any one part can completely determine the other part.

[0080] Therefore, correction methods 1 to 3 are essentially three equivalent implementation paths, achieving the same physical constraints through frequency domain reconstruction, time domain weighting, or time domain symmetric superposition, respectively.

[0081] Similarly, correction methods 4 to 6 are also three equivalent forms, except that the spectrum reconstruction is performed starting from the imaginary part.

[0082] In theory, methods 1–3 and 4–6 are mirror images of each other in the spectral space and can be derived from each other through simple phase reversal and Fourier duality.

[0083] These methods are universal and applicable to all linear time-invariant systems, regardless of whether the data comes from numerical calculations, experimental measurements, or external model interfaces.

[0084] Correction methods 7 and 8 are based on the analytical relationship between logarithmic magnitude and phase, and are usually referred to as the "minimum phase reconstruction" principle.

[0085] For any stable linear system without zeros in the right half-plane, its logarithmic magnitude ln|H(ω)| and phase φ(ω) satisfy the Hilbert duality: φ(ω) = -Hilbert[ln|H(ω)|], ln|H(ω)| = Hilbert[φ(ω)].

[0086] This means that as long as any part of the amplitude or phase information is known, the other part can be reconstructed through the Hilbert transform, thereby recovering the complete minimum phase frequency response.

[0087] Therefore, correction method 7 and correction method 8 can be regarded as two inverse implementation paths. Unlike the previous six methods, these two methods not only achieve causal correction, but also ensure the stability and minimum phase characteristics of the system, which is particularly suitable for situations in measurement or simulation where only one-sided data (amplitude or phase) can be obtained.

[0088] In summary, the correction method based on the relationship between the real and imaginary parts of the frequency response and the correction method based on the amplitude-phase relationship constitute two independent and complementary causal correction systems in mathematics. The former restores the even-odd symmetry of the spectrum through analytical relations, thereby ensuring the causal response of the linear system; the latter reconstructs the minimum-phase system through the consistency of amplitude and phase, thereby ensuring the stability of the time-domain response and energy conservation.

[0089] In practical implementation, this invention allows for adjustments to the Hilbert transform sign, phase sign, or complex exponent direction in the above mathematical formulas based on the phase definition, coordinate direction, or signal convention of the actual system, to adapt to different software frameworks or simulation environments without affecting the essential principle of this invention.

[0090] The computer can automatically select an appropriate correction strategy based on the type, completeness, and noise level of the input frequency response data, or combine multiple methods in a multi-stage calculation process to restore the physical consistency and numerical accuracy of the system to the greatest extent possible while satisfying causal constraints.

[0091] For a system with multiple coupled ports, causality correction can be applied to each port pair sequence of the frequency response array H(ω)[i, j]. The computer performs correction operations on H(ω)[i, j] in the frequency dimension to ensure that the amplitude and phase relationship meets the causality requirement; at the same time, the program maintains the consistency of the spectral structure of all port pairs during the correction process to ensure the overall coordination of the coupling relationship between the ports. When the corrected frequency response is subsequently inversely transformed into the time domain sequence h[n][i, j], it can be guaranteed that the obtained multi-port impulse response meets the causality requirement in the time domain.

[0092] In another implementation, the computer can use mode decomposition to obtain the discrete impulse response sequence h[n][i, j] of the first subsystem. This method does not rely on time-domain integration or frequency-domain inverse transform, but instead utilizes the characteristic structure of the linear system to represent the dynamic behavior of the system as several modal components, and reconstructs the impulse response of each port pair in the time domain by modal superposition.

[0093] In one implementation, for a linear subsystem that can be converted into a first-order state-space form dx / dt = A·x + B·u, the computer can perform eigenvalue decomposition or generalized eigenvalue decomposition on the system matrix A to obtain eigenvalues ​​and their left and right eigenvectors, thereby constructing the system's mode set. When a unit impulse input is applied to the j-th port, the time-domain response of each mode can be calculated analytically and evaluated at each discrete time node n·Δt. Subsequently, the program superimposes the contributions of all modes at each n·Δt to obtain h[n][i, j]. For systems represented by second-order or higher-order semi-discrete equations, the program can introduce additional state variables to convert them into a first-order form before performing the above mode extraction process.

[0094] In another implementation, when the system matrix A cannot be directly obtained or the linear operator exists as a black box, the computer can extract the dominant mode of the system based on the frequency response H(ω)[i, j]. Pole-residual expansion, subspace projection, vector fitting, or similar mode identification methods can be used to identify the dominant poles and residuals of the system from the amplitude and phase information of the frequency response. Based on the obtained modal parameters, the program can generate the corresponding time-domain modal response at each discrete time point and combine them into h[n][i, j] by port pair.

[0095] In a further implementation, to improve the accuracy of the reconstructed impulse response in the high-frequency region or to compensate for dynamic components that have not entered the dominant mode set, the computer can add a residual vector or an equivalent correction vector. This residual part can come from static correction, additional modes, or high-frequency asymptotic analysis, and its time-domain contribution is superimposed with the dominant mode, so that the obtained h[n][i, j] can still maintain good amplitude accuracy, phase consistency, and numerical stability with a finite number of modes.

[0096] Through the aforementioned modal decomposition, modal parameter identification, and residual compensation mechanisms, the computer can obtain the impulse response array h[n][i,j] that meets engineering accuracy requirements at a low computational cost without performing large-scale repeated matrix solving. This array can then be directly applied to convolution decomposition and equivalent equation construction, thereby enabling the reuse of linear dynamic characteristics.

[0097] In summary, the impulse response can be obtained through the time-domain method, the inverse frequency-domain transform method, or the mode decomposition method. All three methods are used to obtain the time-domain characteristics of a linear subsystem to a unit impulse excitation, but they differ in their implementation paths and applicable scenarios. The time-domain method directly calculates the impulse response through numerical integration, with a clear process and compatibility with existing time-domain solvers. It is suitable for systems with small model sizes or those requiring consistency with the existing simulation framework. The inverse frequency-domain transform method, based on the system's frequency response, obtains the time-domain characteristics through numerical inverse transformation. It facilitates parallel processing and can be combined with causal correction, making it suitable for complex structures, scenarios where only frequency-domain data is provided, or scenarios requiring spectral consistency. The mode decomposition method utilizes the characteristic structure of the system matrix or the dominant modes identified by the frequency response to reconstruct the impulse response analytically or semi-analytically. It is suitable for systems of moderate size or whose dynamic characteristics can be accurately characterized by finite modes. Computers can flexibly select or combine these methods according to system size, matrix properties, available data types, and performance requirements to obtain time-domain impulse responses that meet engineering accuracy requirements while ensuring numerical stability.

[0098] In some implementations, the computer can truncate or window the obtained impulse response sequence h[n] to improve its stability and spectral consistency during numerical convolution. Since impulse responses typically have a finite but not strictly truncated time length, directly using the original response within a finite time window may lead to high-frequency ringing or spectral leakage, thus affecting the accuracy of the convolution results. To reduce such truncation effects, the program can apply a smooth weighting function to h[n] within a selected time window, allowing the response to naturally decay to near zero at the end of the window. Suitable windowing methods include exponentially decaying windows, Hanning windows, Kaiser windows, or other weighted window functions with continuous boundary characteristics. By smoothly transitioning the sequence boundaries, the continuity in the time domain can be improved, non-physical oscillations can be suppressed, and the numerical stability of subsequent convolution solutions can be improved while maintaining overall energy consistency.

[0099] During simulation, the impulse response sequence and independent solution sequence may be pre-generated with a sampling step size different from that of the numerical integrator. To ensure that the convolution solution is correctly aligned on the time grid, the computer can resample or interpolate the above sequences before performing the convolution operation to make them consistent with the time step size used in the simulation. In the implementation of variable step integration, the program can dynamically update the resampling results according to the new time step size before the start of each time step, so as to maintain the energy consistency of the convolution terms at different time scales and maintain the numerical stability of the overall solution process.

[0100] Thus, this step completes the calculation of the linear subsystem characteristics, providing the theoretical and numerical foundation for the entire method. The independent solution reflects the system's free response characteristics, while the impulse response describes the transmission law of external excitation or coupled input. Together, they constitute the input terms of the subsequent equivalent equations, enabling the system coupling to be handled analytically. Therefore, after this step, the dynamic characteristics of the first subsystem have been fully characterized, and the system partitioning and linear part solution stage is completed, laying the foundation for the next stage of constructing the equivalent equations.

[0101] In some implementations, to further improve computational efficiency, the computer can perform model order reduction (MOR) on the linear equations before calculating the impulse response. This reduction can include methods such as modal truncation, static condensation, equilibrium truncation, matrix projection, and vector fitting. The simplified system matrix, while preserving the main spectral and dynamic characteristics, still satisfies the linear time-invariant condition, and therefore can be directly used in the impulse response calculation and convolution steps of this invention. This processing does not affect the universality of the method of this invention, making it applicable to both the original high-dimensional system and the linear subsystems after order reduction.

[0102] (Step 103) Construction of equivalent equations After obtaining the independent solution sequence x0[n] and impulse response sequence h[n] of the first subsystem, this step aims to explicitly introduce the linear dynamic effect of the first subsystem in the time domain into the equations of the second subsystem in the form of discrete convolution, thereby constructing an equivalent driving force that can be solved by time-stepping. Through this process, the computer no longer needs to solve the linear equations of the first subsystem during simulation; it can obtain its real-time impact on the second subsystem simply by recursive convolution.

[0103] Given the linearity, the time-domain response of the first subsystem on the discrete-time grid {n·Δt} can be obtained by superimposing the independent solutions with the impulse response. Suppose the second subsystem generates a port input p(y[k], k·Δt) at the k-th time node; then the discrete response of the first subsystem at the n-th time node can be expressed as a discrete convolution: x[n] = x0[n] + ∑(k = 0 to n) h[n - k] · p(y[k], k·Δt).

[0104] Where x0[n] is the sequence of independent solutions of the first subsystem under zero input conditions, h[n - k] is the discrete value of the impulse response at a time lag of n - k, and p(y[k], k·Δt) is the port input of the second subsystem to the first subsystem at the k-th time node.

[0105] To make convolutional structures easier to implement recursively in computers, the above formula can separate the input terms at the current time point. Explicitly writing out the k = n terms in the summation interval, we get: x[n] = x0[n] + ∑(k = 0 to n - 1) h[n - k] · p(y[k], k·Δt) + h[0]· p(y[n], n·Δt) After this separation process, the influence of the first subsystem on the second subsystem is naturally divided into three parts in the discrete-time domain: the independent solution term x0[n], the historical action term ∑(k = 0 to n-1)h[n - k] · p(y[k], k·Δt) formed by the accumulation of past port inputs, and the immediate action term h[0]· p(y[n], n·Δt) directly related to the input at the current time node. This three-term structure is numerically equivalent to the continuous domain expression x(t)=x0(t)+∫H(t-τ)p(y(τ), τ)dτ, which enables the linear coupling across subsystems to maintain causality, computability, and recursiveness in discrete time, laying a clear and stable numerical foundation for the subsequent construction of the equivalent driving force of the second subsystem.

[0106] On the discrete-time grid {n·Δt}, the generalized differential equation of the second subsystem can be expressed as: g(y[n], n·Δt) = q(x[n], n·Δt).

[0107] Here, g(·) represents the discrete form of the generalized operator of the second subsystem mentioned earlier. This operator is used to describe the ordinary differential equations or differential-algebraic equations of the system in the continuous-time domain, and in the discrete domain, it takes the state variables at the interpolation point, the discrete approximation of their time derivatives, and possible algebraic constraints as inputs. Therefore, it can uniformly cover the expression structures of ODE and DAE in discrete time. q(·) is the port function, used to describe the effect of the first subsystem on the second subsystem at the current time node.

[0108] Based on the aforementioned convolution solution, the response of the first subsystem at time node n satisfies x[n] = x0[n] + ∑(k = 0 to n-1) h[n - k] · p(y[k], k·Δt) + h[0] ·p(y[n], n·Δt).

[0109] In numerical implementation, both independent solution terms and historical action terms can be directly calculated before solving for the unknowns of the second subsystem, and thus can be considered as known equivalent inputs. However, the value of the instantaneous action term h[0]·p(y[n], n·Δt) depends on y[n], and is not yet determined at the current time point, so it cannot be treated as an external known quantity. To ensure the overall causality and solvability of the equation, the computer directly incorporates this instantaneous term into the left-hand side of the operator g(·) when constructing the equivalent equation, making it an additional linear action or equivalent operator increment within the second subsystem. Through this incorporation, the discrete equations of the second subsystem are structurally consistent with its original generalized differential equations, thus without changing the existing time integration framework, and can continue to be solved using implicit or explicit ODE / DAE integration methods.

[0110] In a preferred embodiment, the port function is presented as a linear linear term structure q(x[n], t) = Q·x[n] by dividing the subsystem. The port function of the second subsystem at the nth time node can then be written as... Q·x[n] = Q·x0[n] + ∑(k = 0 to n-1) Q·h[n - k] · p(y[k], k·Δt) + Q·h[0]·p(y[n], n·Δt).

[0111] In this form, both independent solution terms and historical action terms can be used as known inputs on the right-hand side, while the instantaneous term Q·h[0]·p(y[n], n·Δt) can be merged with the linear operator of the second subsystem itself to form the updated equivalent operator. This approach not only maintains the consistency between the equation structure and the original model, but also allows the coupling terms across subsystems to be injected into the time-step update process in an explicit and recursive manner. Since the linear dynamic action of the first subsystem has been exposed through convolution, the solution of the second subsystem at each time node no longer depends on the repeated solution of the equations of the first subsystem, thereby significantly reducing the computational burden and improving the overall simulation efficiency and numerical stability while maintaining spectral fidelity.

[0112] The independent solution sequence x0[n] and impulse response sequence h[n] of the first subsystem are usually pre-calculated at a fixed sampling step size Δt_imp, while the time integration process of the second subsystem may use a different step size Δt_sim. In variable step-size integration algorithms, Δt_sim may even change dynamically with time. If the pre-calculated sequence is directly used for convolution, the time node of the convolution operation will be inconsistent with the evolution node of the integrator, thus introducing problems such as phase shift, energy imbalance, or numerical ringing. Therefore, before the convolution calculation begins, the computer performs interpolation or resampling on the original sequence x0[n] and h[n] based on the current simulation step size Δt_sim to generate sequences x0′[n] and h′[n] that are strictly aligned with the simulation time node. This interpolation process can use linear interpolation, spline interpolation, or an equivalent resampling method performed in the frequency domain. Under variable step size conditions, the program automatically updates x0′[n] and h′[n] before the start of each time step, thereby ensuring that the time base of the convolution operation is always consistent with the integrator, so that the physical meaning and numerical stability of discrete convolution are maintained under different step sizes.

[0113] For a first subsystem with multiple input ports, its effect on the second subsystem is not a single-input, single-output relationship, but rather determined by the multi-input, multi-output relationship between the ports. To obtain the complete dynamic coupling characteristics between the ports, the computer applies a unit impulse to each input port sequentially during the initialization phase, while simultaneously recording the responses of all output ports, thus obtaining the multi-port impulse response matrix h(t). The columns of the matrix correspond to the input ports, and the rows correspond to the output ports. Its discrete form h[n] describes the dynamic transmission relationship between different ports in the time domain. During the simulation, the convolution operation is performed in parallel along the port dimension. The historical input sequence corresponding to each input port is convolved with the corresponding column of the impulse response matrix one by one, and then the responses of all ports are summed to obtain the comprehensive response of the first subsystem at each output port. Since the number of ports is usually much smaller than the system degrees of freedom, the convolution operation only occurs in the low-dimensional domain of the port space, which is naturally suitable for parallel acceleration by CPU or GPU. In this way, the coupling effect of the first subsystem in the multi-port case can still be clearly expressed by the convolution decomposition structure and can be stably transmitted to the second subsystem in the discrete time domain.

[0114] After constructing the convolution decomposition structure described above, we can further analyze the computational efficiency advantages of this invention from the perspective of algorithm complexity. In existing technologies, for dynamic systems that simultaneously contain linear and nonlinear couplings, direct numerical integration methods often require solving the system of linear equations (such as LU decomposition or other matrix decompositions) on the global matrix of the system at each time step. If the total degrees of freedom of the system are N, then a single-step solution typically has a computational complexity of O(N³), and its computational load increases cubically with the system size. For large-scale engineering models, such as high-dimensional finite element systems, circuit network models, or other large-scale dynamic simulations, repeatedly performing O(N³) matrix solutions leads to extremely high computational overhead, thereby limiting real-time performance and simulation speed.

[0115] This invention expresses the dynamic behavior of the linear component as a convolutional form by pre-calculating the independent solutions and impulse responses of the first subsystem. This eliminates the need for solving linear equations in subsequent time steps, requiring only convolution operations at the port dimension. The computational complexity of this convolution depends primarily on the number of ports P and the effective length L of the impulse response, with a single-step complexity approximately O(P²·L). Since P is typically much smaller than N, and the dynamic changes of the linear subsystem are already manifested in a fixed impulse response form, the overall complexity O(P²·L) of this invention is significantly lower than the traditional method's O(N³) when the system size is large. This achieves an order-of-magnitude or even exponential efficiency improvement while maintaining numerical accuracy and spectral consistency.

[0116] In summary, after this step, all linear dynamic actions of the first subsystem have been fully mapped to the discrete equivalent equations of the second subsystem through a combination of independent solutions, historical convolution terms, and instantaneous terms. The second subsystem can then be solved independently within the unified discrete equation framework during subsequent time-step integration, no longer relying on repeated solutions to the linear operators of the first subsystem. This equivalent equation maintains the structural characteristics of the original generalized differential equations while significantly reducing the computational burden, thereby improving overall simulation efficiency and numerical stability while preserving physical consistency. After constructing the equivalent equations, the computer enters the time integration stage to solve the discrete equivalent equations of the second subsystem and generate the final time-domain response results.

[0117] (Step 104) Time-domain numerical solution and result output After obtaining the discrete equivalent equation of the second subsystem, the computer enters the time-domain numerical solution stage. This step involves performing time-step integration on the equivalent equation to progressively update the state variables y(t) of the second subsystem and generate the system's time-domain response. Since the linear dynamic action of the first subsystem has been explicitly expressed in step 103 through convolution, the equations of the second subsystem retain the structure of generalized differential equations and can be directly solved recursively using existing ODE / DAE numerical integration methods, without needing to repeatedly solve the linear part of the first subsystem in the time-step loop.

[0118] When performing integration, the program can choose between explicit or implicit numerical integration methods based on system characteristics. Explicit methods have lower computational cost and are suitable for non-rigid or low-to-medium frequency systems; implicit methods offer better stability and are suitable for situations involving high-frequency components, rigid coefficients, or algebraic constraints. The computer can automatically select an integration strategy based on model characteristics during the initialization phase, or it can employ a specific integration algorithm based on user configuration.

[0119] Before time progression begins, the program sets necessary initial conditions based on the generalized differential equations of the second subsystem, including y(0), dy / dt(0), d²y / dt²(0), or other internal variables required by the integrator. These initial values ​​can be derived from system definitions, external input files, steady-state solutions, or user-specified values. The method of setting these initial values ​​falls under conventional time-domain integration requirements and does not limit the scope of this invention.

[0120] This method is compatible with both fixed-step and variable-step integrators. For the variable-step case, before each time step, the program invokes the resampling mechanism described in step 103 based on the current step size Δt_sim to perform interpolation or resampling on the discrete data of the pre-calculated independent solutions x0[n] and impulse responses h[n], ensuring that their time points are consistent with the current integration step size. Through this mechanism, convolutional terms and instantaneous terms can maintain time alignment and energy consistency under different step sizes, thereby avoiding numerical deviations caused by step size changes.

[0121] In some implementations, when it is necessary to recover the continuous state variables x(t) of the first subsystem, the program can reconstruct them in two ways. First, based on the convolutional structure constructed in step 103, the solution results y[n] of the second subsystem at each time node are substituted into the discrete convolution expression. x[n] = x0[n] + ∑(k = 0 to n-1) h[n - k] · p(y[k], k·Δt) + h[0] ·p(y[n], n·Δt), The response sequence of the first subsystem on the discrete-time grid can then be obtained. If a continuous-time form is required, the sequence can be reconstructed using linear interpolation, spline interpolation, or other numerical methods. Secondly, after the second subsystem completes time progression, since y(t) has been fully obtained over the entire time domain, the original generalized differential equation of the first subsystem can be obtained. f(x(t), dx / dt(t), d²x / dt²(t), …, t) = p(y(t), t) The problem degenerates into an initial value problem with known input terms. The program can perform a conventional numerical integration on the equation during the post-processing stage to obtain continuous x(t). This integration process does not participate in the main simulation flow and does not increase the computational burden of the solution stage of the second subsystem. The two recovery paths are mathematically equivalent.

[0122] For a second subsystem containing nonlinear terms or algebraic constraints, the program can solve the local nonlinear equations at each time step using implicit iteration or the Newton-Raphson method. Since the linear effects of the first subsystem are already explicit through the convolution structure, the iterative process only involves the variables within the second subsystem, thereby reducing the number of linear equation solutions and improving overall computational efficiency. For large-scale models, both convolution operations and local integration steps can be executed in parallel in a multi-threaded or multi-core environment.

[0123] At the end of the simulation, the program stores or exports the time-domain data of state variables, port variables, and output variables, which may include time series, energy and power indicators, or other result formats required by the user. Through the above-described time-domain numerical solution and output mechanism, the time-domain evolution process of the system can be completely reconstructed, and together with the convolution solution step, it constitutes an efficient, stable, and universally applicable numerical solution framework for dynamic systems.

[0124] like Figure 2 As shown, in one embodiment, the method of the present invention is implemented by a program running on a computer. This program constitutes a dynamic system numerical solution system 200, which performs functions such as system partitioning, linear subsystem characteristic analysis, equivalent equation construction, and time-domain solution through internal modules, thereby realizing dynamic system simulation based on generalized differential equations. System 200 includes a system partitioning module 201, a linear subsystem characteristic calculation module 202, an equivalent equation construction module 203, and a time-domain numerical solution and result output module 204.

[0125] The system partitioning module 201 is used to divide the dynamic system to be solved into a first subsystem and a second subsystem, and establish the port function relationship between the two. During the program initialization phase, this module reads the coefficient matrix or operator parameters of the system equations and determines the scope and interface conditions of the subsystems according to a preset strategy.

[0126] The linear subsystem characteristic calculation module 202 is used to calculate the linear characteristics of the first subsystem, including independent solutions and impulse responses, without considering coupled inputs. This module can select time-domain methods, inverse frequency-domain transform methods, mode decomposition methods, or other equivalent approaches based on available data, and perform causal corrections or windowing as necessary to obtain impulse responses that meet numerical stability requirements.

[0127] The equivalent equation construction module 203 is used to establish equivalent interaction relationships between subsystems based on the impulse response. This module generates independent solution terms, immediate action terms, and historical action terms, and substitutes them into the generalized differential equation of the second subsystem to form a directly solvable equivalent equation. Specifically, when calculating the historical action terms, the computer uses discrete convolution to weight and accumulate the port input quantities of the second subsystem at each historical time step with the impulse response sequences of the corresponding time differences, thereby obtaining the cumulative action across time. Through this method, the linear action of the first subsystem is explicitly expressed in the independent solution terms, historical action terms, and immediate action terms, maintaining the recursive structure of the dynamic equations of the second subsystem, enabling direct time advancement using numerical integration methods in subsequent steps.

[0128] The time-domain numerical solution and result output module 204 is used to perform time-step integration on the equivalent equation and output the system's time-domain response. This module can call integrators such as the Euler method, trapezoidal method, Runge-Kutta method, or generalized-α method, update state variables and port variables during the integration process, and display, record, or export the results.

[0129] The modules in System 200 can be implemented as independent software units, or integrated as functions or objects within the same program framework, and communicate through a unified data structure and interface. System 200 can run in standalone, parallel, or distributed computing environments, and can utilize multi-core processors or acceleration hardware to improve the efficiency of impulse response computation and convolution operations.

[0130] When the above program runs in the computer, the processor executes the functions of each module in sequence, enabling the computer to complete... Figure 1 The numerical solution process is shown and implemented. Figure 2 The software system structure shown enables the modeling of dynamic systems, calculation of linear characteristics, construction of equivalent equations, and time-domain solution. The above modules can be implemented through software, hardware, or a combination of both; their specific implementation does not constitute a limitation of this invention.

[0131] like Figure 3 As shown, to illustrate one possible implementation, the present invention can be implemented on a computer device 300. The computer device 300 includes a processor, a memory, and a non-transitory computer-readable storage medium storing program instructions. When the program is executed by the processor, it causes the computer to perform... Figure 1The solution process is shown below, and the solution is completed. Figure 2 The module functions are shown.

[0132] The computer device may be a general-purpose computer, a server, an embedded controller, or a cloud computing node. The hardware structure is not limited to a specific form. As long as it can execute the program instructions and implement the method steps of the present invention, it is considered to fall within the scope of the present invention. Figure 3 The description of the computer device herein is merely illustrative, intended to demonstrate possible forms of implementation of the invention, and does not limit the hardware composition in any way. For those skilled in the art, any technical solution or effect substantially identical to that of the present invention achieved through hardware circuits, software programs, programmable logic, or a combination of hardware and software should be considered within the scope of protection of the present invention.

[0133] In summary, the method of this invention can be implemented in a computer, and its steps can be executed by software, hardware, or a combination of both. Once executed... Figure 1 The numerical solution process is shown and implemented. Figure 2 The software structure shown should be considered to fall within the protection scope of this invention.

Claims

1. A method for numerically solving dynamic systems described by generalized differential equations in a computer, characterized in that, include: In a computer, the dynamic system is divided into a first subsystem and a second subsystem, and the coupling relationship between the two is established through port functions; In a computer, the independent solution sequence of the first subsystem and the impulse response sequence used to characterize its effect on the port input are obtained without considering the coupled input. In the computer, independent solution terms are generated based on the independent solution sequence, and time-domain convolution operation is performed on the impulse response sequence and the port input formed by the port function of the second subsystem to generate the historical action terms and immediate action terms of the first subsystem on the second subsystem. Based on this, an equivalent equation of the second subsystem containing the above three terms is constructed, so that the equivalent equation can be solved numerically directly in the computer. In a computer, time-step integration is performed on the equivalent equation to obtain the time-domain numerical solution of the dynamic system. The generalized differential equations include ordinary differential equations and differential algebraic equations. The port function can be linear, nonlinear, or time-varying. The impulse response is used to characterize the time-domain response of the first subsystem under port input disturbance.

2. The method according to claim 1, wherein, The method is used in computers for simulation calculations of dynamic systems, which include: system-level models established through lumped parameter modeling, including circuit network models, hydraulic system models, mechanical or multibody system dynamic models, control system models, and chemical process systems; field-level models formed by discretizing using finite element, finite volume, finite difference, boundary element, or spectral element methods to form semi-discrete equations, including structural, thermal, electromagnetic, acoustic, or fluid dynamic models; and comprehensive models formed by the coupling of multiple physical domains. The method is applicable to linear systems, nonlinear systems, and dynamic systems whose parameters change over time.

3. The method according to claim 1, wherein, The impulse response is obtained by applying an approximate impulse-type test input to the coupling port of the first subsystem in a computer and recording the port output generated during the time-domain numerical solution process to construct the impulse response.

4. The method according to claim 1, wherein, The impulse response is obtained by performing an inverse Fourier transform on the frequency response characterizing the coupling port relationship of the first subsystem in a computer to generate the corresponding time-domain impulse response; wherein the frequency response data can be obtained by calculation, experimental measurement, or imported from an external file.

5. The method according to claim 4, characterized in that, In the inverse transformation of the frequency response, causal correction is performed on the frequency response itself or the impulse response generated therefrom, so that the resulting time-domain impulse response satisfies the causal condition in time, thus avoiding the generation of non-physical responses before the excitation occurs. The causal correction is performed in a computer and can be achieved by using the Hilbert transform (or equivalent spectral symmetry) based on the relationship between the real and imaginary parts of the frequency response, the Hilbert transform (also known as minimum phase reconstruction) based on the relationship between the logarithmic magnitude and phase, or a combination of both.

6. The method according to claim 1, wherein, The impulse response is obtained by performing mode decomposition on the first subsystem, obtaining mode parameters in a computer to characterize its linear dynamic characteristics, and reconstructing the impulse response in the time domain based on one or more selected dominant modes; wherein the mode parameters can be obtained by eigenvalue decomposition, generalized eigenvalue decomposition, state-space pole-residual expansion, subspace projection, vector fitting or other equivalent mode identification methods; and the mode decomposition can be combined with one or more residual vectors or equivalent correction vectors to compensate for dynamic characteristics that do not participate in mode reconstruction and improve the reconstruction accuracy of the impulse response.

7. The method according to claim 1, wherein, The impulse response is truncated or windowed in the computer to improve numerical stability, suppress high-frequency ringing, and reduce spectral leakage. The windowing process can employ an exponential decay window, Hanning window, Kaiser window, or other weighted window functions with continuous boundary characteristics to improve time-domain smoothness while maintaining the energy consistency of the impulse response and to allow the impulse response to naturally decay to near zero at the end of the time window, thereby avoiding non-physical oscillations caused by truncation errors.

8. The method according to claim 1, wherein, When the simulation time step is inconsistent with the sampling step of the impulse response, the computer resamples or interpolates the impulse response and the independent solution sequence before performing the convolution operation to make their time steps consistent. Under variable step size integration, the computer can dynamically update the resampling results according to the current time step before the start of each simulation step, so as to ensure the energy consistency and numerical stability of the convolution calculation.

9. A computer-readable storage medium having stored thereon computer program instructions that, when executed on a computer, cause the computer to perform the method as described in any one of claims 1 to 8.

10. A computer device comprising a processor and a memory, the memory storing computer program instructions that, when executed by the processor, cause the processor to perform the method as claimed in any one of claims 1 to 8.

Images