CN106294283A

CN106294283A - Temporal basis functions fast algorithm based on Taylor series expansion

Info

Publication number: CN106294283A
Application number: CN201510266484.9A
Authority: CN
Inventors: 陈如山; 李威; 丁大志; 樊振宏; 程光尚
Original assignee: Nanjing University of Science and Technology
Current assignee: Nanjing University of Science and Technology
Priority date: 2015-05-22
Filing date: 2015-05-22
Publication date: 2017-01-04
Anticipated expiration: 2035-05-22
Also published as: CN106294283B

Abstract

The invention discloses a fast algorithm of time domain integral equation based on Taylor series expansion. Traditional time-domain integral equation methods cannot solve large-scale electromagnetic problems due to the limitation of computation time and memory consumption. The development of fast algorithms is the only way to solve practical engineering problems. The present invention utilizes the Taylor series expansion of the Green's function in free space to reconstruct the form of aggregation, transfer, and configuration to realize the accelerated calculation of matrix-vector multiplication. The present invention requires less computational memory and computational time for solving metal complex target scattering problems . Moreover, programming is relatively simple and easy to implement, and has strong practical engineering application value.

Description

A Fast Algorithm for Time Domain Integral Equations Based on Taylor Series Expansion

技术领域 technical field

本发明涉及金属目标瞬态电磁散射特性数值计算技术，属于目标电磁散射特性的快速计算技术领域，具体是一种基于泰勒级数展开的时域积分方程快速算法。 The invention relates to a numerical calculation technology of transient electromagnetic scattering characteristics of a metal target, belongs to the technical field of rapid calculation of electromagnetic scattering characteristics of a target, and specifically relates to a fast algorithm of time-domain integral equations based on Taylor series expansion.

背景技术 Background technique

随着计算电磁学领域研究的深入，传统的频域方法已经不能满足需要。以计算机硬件技术的发展为契机，人们逐步具有了直接在时域对具有宽频带特性的瞬变电磁场的计算分析能力，从而实现了对物理量和物理现象更深刻、更直观的理解。同频域方法相比，在时域求解目标的电磁特性不仅可以直观的揭示目标与电 With the deepening of research in the field of computational electromagnetics, the traditional frequency domain methods can no longer meet the needs. Taking the development of computer hardware technology as an opportunity, people gradually have the ability to calculate and analyze transient electromagnetic fields with broadband characteristics directly in the time domain, thus realizing a deeper and more intuitive understanding of physical quantities and physical phenomena. Compared with the frequency domain method, solving the electromagnetic characteristics of the target in the time domain can not only intuitively reveal the relationship between the target and the electromagnetic

磁波相互作用的机理，而且通过少量的计算就可以获得目标的宽频带信息，这在宽带电磁问题、瞬态电磁问题分析中具有明显的优势。 The mechanism of magnetic wave interaction, and the broadband information of the target can be obtained through a small amount of calculation, which has obvious advantages in the analysis of broadband electromagnetic problems and transient electromagnetic problems.

同基于微分方程的时域方法(FDTD、FETD等)相比，基于积分方程的时域方法在求解的未知数数量上具有明显的优势，这是因为积分方程利用格林函数建立源和场的关系，求解区域在边界上，离散后未知数的数量与边界面积成正比。其次，积分方程方法自动满足辐射边界条件，不需要强加吸收边界，而吸收边界是基于微分方程的方法所必需的。 Compared with time-domain methods based on differential equations (FDTD, FETD, etc.), time-domain methods based on integral equations have obvious advantages in the number of unknowns to be solved, because the integral equation uses Green's function to establish the relationship between source and field, The solution area is on the boundary, and the number of unknowns after discretization is proportional to the boundary area. Second, the integral equation method automatically satisfies the radiation boundary condition without imposing an absorbing boundary, which is required for differential equation-based methods.

由于计算机水平和硬件的限制，传统的时域积分方程方法无法解决大规模的电磁问题。而发展快速算法是解决实际工程问题的必由之路。时域积分方程快速算法的研究是在近20年开展起来的，其中最著名的两种：一是时域平面波算法(PWTD)，二是时域自适应积分方程(TD-AIM)。对于三维目标模型，TD-AIM的计算量和存储量分别为和可见其处理三维目标的效率低下。PWTD算法其原理类似于FMM算法，但是因为时间变量的存在，使得PWTD算法更加复杂，实现起来更加困难。并且该算法的适应性差，当用以求解介质或者有耗媒质等复杂问题时，还需要对算法进行特殊修改。 Due to the limitations of computer level and hardware, the traditional time-domain integral equation method cannot solve large-scale electromagnetic problems. The development of fast algorithms is the only way to solve practical engineering problems. The research on fast algorithms for time-domain integral equations has been carried out in the past 20 years, and the two most famous ones are: one is the time-domain plane wave algorithm (PWTD), and the other is the time-domain adaptive integral equation (TD-AIM). For the 3D target model, the calculation amount and storage amount of TD-AIM are respectively and It can be seen that its efficiency in dealing with three-dimensional objects is low. The principle of the PWTD algorithm is similar to the FMM algorithm, but due to the existence of time variables, the PWTD algorithm is more complicated and more difficult to implement. Moreover, the adaptability of the algorithm is poor. When it is used to solve complex problems such as media or lossy media, special modifications to the algorithm are required.

发明内容 Contents of the invention

本发明的目的在于提供一种高效、稳定的分析金属目标的时域电磁散射特性的快速方法。由于对远场部分采用泰勒级数展开重构成聚合、转移、配置的形式进而实现矩阵矢量乘的加速计算，本发明对于求解金属复杂目标散射问题需要更少的计算内存以及计算时间。而且编程相对简单易于实现。 The purpose of the present invention is to provide an efficient and stable fast method for analyzing the time-domain electromagnetic scattering characteristics of metal targets. Since the far-field part is reconstructed into the forms of aggregation, transfer, and configuration by using Taylor series expansion to realize accelerated calculation of matrix-vector multiplication, the present invention requires less computing memory and computing time for solving metal complex target scattering problems. And programming is relatively simple and easy to implement.

实现本发明目的的技术解决方案为：一种基于泰勒级数展开的时域积分方程快速算法，步骤如下： The technical solution that realizes the object of the present invention is: a kind of time-domain integral equation fast algorithm based on Taylor series expansion, the steps are as follows:

第一步，建立时域电磁场积分方程。金属目标在入射电磁波的作用下，在金属散射体的表面上会产生面感应面电流，感应电流的不断变化，会向外辐射出电磁场，由感应电流产生的电磁场就是散射场，利用入射电磁场和散射电磁场在金属表面满足的边界条件建立时域电磁场积分方程。 The first step is to establish the time-domain electromagnetic field integral equation. Under the action of the incident electromagnetic wave, the metal target will generate a surface induction surface current on the surface of the metal scatterer. The continuous change of the induction current will radiate the electromagnetic field outward. The electromagnetic field generated by the induction current is the scattering field. Using the incident electromagnetic field and The boundary conditions of the scattered electromagnetic field on the metal surface are established to establish the integral equation of the electromagnetic field in the time domain.

第二步，将散射体表面上离散得到的子散射体分组。任意两个子散射体间的互耦或自耦根据它们所在组的位置关系而分成近场组对和远场组对。当它们是近场组对时，采用直接数值计算。而当它们为远场组对时，则采用泰勒级数展开成聚合-转移-配置方法计算。 The second step is to group the discrete sub-scatterers on the surface of the scatterer. The mutual coupling or autocoupling between any two sub-scatterers can be divided into near-field pair and far-field pair according to their positional relationship. When they are near-field pairs, direct numerical calculations are used. And when they are far-field pairings, Taylor series expansion is used to calculate by aggregation-transfer-configuration method.

第三步，近场阻抗矩阵计算，将金属表面电流密度用空间基函数和时间基函数展开，并在空间域上进行伽辽金测试，时间域上进行点匹配得到矩阵元素值。 The third step is the calculation of the near-field impedance matrix. The current density of the metal surface is expanded by the space basis function and the time basis function, and the Galerkin test is performed in the space domain, and the point matching is performed in the time domain to obtain the matrix element values.

第四步，远场部分采用泰勒级数展开重构成聚合、转移、配置的形式。 In the fourth step, the far-field part is reconstructed into the form of aggregation, transfer, and configuration using Taylor series expansion.

第五步，矩阵方程求解以及电磁散射参数的计算。利用时间递推的方式求解每个时刻的电流系数。每个时刻空间某处的散射场贡献来源分为近场组源的贡献和远场组源的贡献。前者直接通过近场矩阵和电流系数相乘得到，后者通过泰勒级数展开成聚合-转移-配置方法快速计算。采用迭代法求解出最终的散射电流系数。并根据互易定理求解出雷达散射截面积。 The fifth step is to solve the matrix equation and calculate the electromagnetic scattering parameters. The current coefficient at each moment is solved by time recursion. The source of scattered field contribution somewhere in space at each moment is divided into the contribution of the near-field group source and the contribution of the far-field group source. The former is directly obtained by multiplying the near-field matrix and the current coefficient, and the latter is quickly calculated by the aggregation-transfer-configuration method through Taylor series expansion. The final scattering current coefficient is solved by an iterative method. And according to the reciprocity theorem, the radar cross-sectional area is obtained.

本发明与现有技术相比，其显著优点为：(1)可以快速分析复杂电大目标时域电磁散射特性；(2)远场计算是基于泰勒级数展开的，避免了PWTD算法需要的复杂的时频域转换、谱域积分等操作，编程相对简单易于实现；(3)适用性高，可以只需作较少的改动便可应用于求解介质、有耗、色散等复杂问题的分析中。 Compared with the prior art, the present invention has the following significant advantages: (1) time-domain electromagnetic scattering characteristics of complex electrically large targets can be quickly analyzed; (2) far-field calculation is based on Taylor series expansion, which avoids the complexity required by the PWTD algorithm Operations such as time-frequency domain conversion, spectral domain integration, etc., programming is relatively simple and easy to implement; (3) High applicability, it can be applied to the analysis of solving complex problems such as medium, lossy, and dispersion with only a few changes .

附图说明 Description of drawings

图1是本发明的应用于基于泰勒级数展开的时域积分方程快速算法中的多层分组示意图。 Fig. 1 is a schematic diagram of multi-layer grouping applied in the fast algorithm of time domain integral equation based on Taylor series expansion according to the present invention.

图2是本发明的应用于基于泰勒级数展开的时域积分方程快速算法中的近远场关系示意图。 Fig. 2 is a schematic diagram of the near-far field relationship applied in the fast algorithm of the time-domain integral equation based on Taylor series expansion of the present invention.

图3是本发明的应用于基于泰勒级数展开的时域积分方程快速算法中的场源基函数位置向量分布示意图。 Fig. 3 is a schematic diagram of the position vector distribution of field source basis functions applied in the fast algorithm of time domain integral equation based on Taylor series expansion according to the present invention.

图4是本发明的应用于基于泰勒级数展开的时域积分方程快速算法中的简易导弹模型不同频点处双站RCS比较图。(a)30MHz；(b)150MHz；(c)270MHz。 Fig. 4 is a comparison diagram of two-station RCS at different frequency points of the simple missile model applied in the fast algorithm of time domain integral equation based on Taylor series expansion according to the present invention. (a) 30MHz; (b) 150MHz; (c) 270MHz.

具体实施方式 detailed description

下面结合附图对本发明作进一步详细描述。 The present invention will be described in further detail below in conjunction with the accompanying drawings.

本发明一种基于泰勒级数展开的时域积分方程快速算法，步骤如下： A kind of fast algorithm of time domain integral equation based on Taylor series expansion of the present invention, the steps are as follows:

第一步，建立时域电磁场积分方程。 The first step is to establish the time-domain electromagnetic field integral equation.

金属目标在入射电磁波{E^inc(r,t),H^inc(r,t)}的作用下，在金属散射体的表面上会产生面感应面电流J(r,t)，感应电流的不断变化，会向外辐射出电磁场，由感应电流产生的电磁场就是散射场{E^sca(r,t),H^sca(r,t)}，利用入射电磁场和散射电磁场在金属表面满足的边界条件建立时域电磁场积分方程。 Under the action of the incident electromagnetic wave {E ^inc (r, t), H ^inc (r, t)}, the metal target will generate a surface induced surface current J(r, t) on the surface of the metal scatterer, and the induced current continues changes, the electromagnetic field will be radiated outward, and the electromagnetic field generated by the induced current is the scattered field {E ^sca (r,t),H ^sca (r,t)}, which is established by using the boundary conditions satisfied by the incident electromagnetic field and the scattered electromagnetic field on the metal surface Integral Equations of Electromagnetic Fields in the Time Domain.

$\begin{matrix} α α \overset{^^}{n no} ((r r)) \times \times \overset{^^}{n no} ((r r)) \times \times {E E.}^{i i n no c c} ((r r,, t t)) + + ((11 - - α α)) η η \overset{^^}{n no} ((r r)) \times \times {H h}^{i i n no c c} ((r r,, t t)) \\ = = {αL α L}_{e e} {{J J ((r r,, t t))}} + + η η ((11 - - α α)) {L L}_{h h} {{J J ((r r,, t t))}} \end{matrix} - - - - - - ((11))$

$\begin{matrix} {L L}_{e e} {{J J ((r r,, t t))}} = = \overset{^^}{n no} ((r r)) \times \times \overset{^^}{n no} ((r r)) \times \times \frac{μ μ}{44 π π} &Integral; &Integral; {&Integral; &Integral;}_{s the s} d d S S \frac{11}{R R} \frac{\partial \partial J J (({r r}^{' '},, τ τ))}{\partial \partial t t} - - \\ n no ((r r)) \times \times \overset{^^}{n no} ((r r)) \times \times \frac{11}{44 π π ϵ ϵ} &Integral; &Integral; {&Integral; &Integral;}_{s the s} d d S S {&Integral; &Integral;}_{00}^{τ τ} \frac{{&dtri; &dtri;}^{' '} \cdot &Center Dot; J J (({r r}^{' '},, τ τ))}{R R} d d t t \end{matrix} - - - - - - ((22))$

${L L}_{h h} {{J J ((r r,, t t))}} = = \frac{11}{22} J J ((r r,, t t)) - - \overset{^^}{n no} ((r r)) \times \times \frac{11}{44 π π} &Integral; &Integral; {&Integral; &Integral;}_{s the s} d d S S &dtri; &dtri; \times \times \frac{J J (({r r}^{' '},, τ τ))}{R R} - - - - - - ((33))$

其中是单位外法向矢量，α取值为0～1之间的实数，μ和ε分别是自由空间的磁导率和介电参数。η是自由空间的波阻抗。R＝|r-r'|,c是自由空间中的光速,τ＝t-R/c是延时。 in is the unit external normal vector, α is a real number between 0 and 1, μ and ε are the permeability and dielectric parameters of free space, respectively. η is the wave impedance of free space. R=|r-r'|, c is the speed of light in free space, τ=tR/c is the time delay.

第二步，将散射体表面上离散得到的子散射体分组。 The second step is to group the discrete sub-scatterers on the surface of the scatterer.

我们借鉴MLFMM里面对空间基函数进行分组的思想，将基函数的相互作用转换为组的相互作用。假设空间中有一个刚好可以把整个目标物体包围的立方体，把这个立方体等分成8个子立方体，接着在将每个子立方体等分成8个更小的立方体，如图1(a-d)所示。依次类推，直到达到预先设置的门限值，停止划分。定义第一层为划分最细的一层，第一层中的立方体中如果包含基函数，定义为第一层组，对于其它层组，依次类推即可。 We learn from the idea of grouping spatial basis functions in MLFMM, and transform the interaction of basis functions into the interaction of groups. Assuming that there is a cube in the space that can just surround the entire target object, divide the cube into 8 sub-cubes, and then divide each sub-cube into 8 smaller cubes, as shown in Figure 1 (a-d). By analogy, until the preset threshold is reached, the division is stopped. Define the first layer as the most finely divided layer. If the cube in the first layer contains basis functions, it is defined as the first layer group. For other layer groups, it can be deduced in turn.

任意两个子散射体间的互耦或自耦根据它们所在组的位置关系而分成近场组对和远场组对。对于每一层，我们都会设定一个参数β,4＜β＜6。首先从最高层N_l的各组着手划分，如果两组中心间的距离大于规定的的门限值β，就把其称之为远场组对；然后在其子层中，组对中心间距大于规定的该层的门限值，同时这两个组都没有划分到N_l层的远场组对，那么就规定这样的组对为该层的远场组对；根据以上所述依次类推，对于每一层我们就可以得到其远场组对的信息，如图2所示。当它们是近场组对时，采用直接数值计算。而当它们为远场组对时，则采用泰勒级数展开成聚合-转移-配置方法计算。 The mutual coupling or autocoupling between any two sub-scatterers can be divided into near-field pair and far-field pair according to their positional relationship. For each layer, we will set a parameter β, 4<β<6. First, divide the groups from the highest layer N _l , if the distance between the centers of the two groups is greater than the specified threshold value β, it is called a far-field group pair; then in its sub-layers, the group-to-center distance greater than the specified threshold value of this layer, and these two groups are not divided into the far-field group pair of the _N1 layer, then it is stipulated that such a group pair is the far-field group pair of the layer; , for each layer we can get the information of its far-field pair, as shown in Figure 2. When they are near-field pairs, direct numerical calculations are used. And when they are far-field pairings, Taylor series expansion is used to calculate by aggregation-transfer-configuration method.

第三步，近场阻抗矩阵计算。 The third step is to calculate the near-field impedance matrix.

把金属表面电流电流密度J(r,t)用RWG基函数Λ_n(r)作为空间基函数，三角基函数T_j(t)作时间基函数展开： The current density J(r,t) on the metal surface is expanded using the RWG basis function Λ _n (r) as the space basis function, and the triangular basis function T _j (t) as the time basis function:

$J J ((r r,, t t)) = = {Σ Σ}_{n no = = 11}^{{N N}_{S S}} {Σ Σ}_{j j = = 11}^{{N N}_{t t}} {I I}_{n no,, j j} {Λ Λ}_{n no} ((r r)) {T T}_{j j} ((t t)) - - - - - - ((44))$

其中N_s是散射体包含的RWG基函数的个数，N_t是时间基函数的个数，I_n,j是第n个RWG基函数第j个时间步上时间基函数的系数。将式(4)代入(1)，在空间上进行Galerkin测试，在时间上进行点匹配，最终得到如下的矩阵方程： where N _s is the number of RWG basis functions included in the scatterer, N _t is the number of time basis functions, I _n,j is the coefficient of the time basis function of the nth RWG basis function at the jth time step. Substitute formula (4) into (1), perform Galerkin test in space, and point matching in time, and finally get the following matrix equation:

${Z Z}^{00} {I I}^{i i} = = {V V}^{i i} - - {Σ Σ}_{j j = = 11}^{i i - - 11} {Z Z}^{i i - - j j} {I I}^{j j} - - - - - - ((55))$

${Z Z}^{i i - - j j} = = {αZ αZ}_{E E.}^{i i - - j j} + + ((11 - - α α)) {Z Z}_{M m}^{i i - - j j} - - - - - - ((66))$

${[[{Z Z}_{E E.}^{i i - - j j}]]}_{m m n no} = = \{\begin{matrix} \frac{{μ μ}_{00}}{44 π π} {&Integral; &Integral;}_{{S S}_{m m}} {&Integral; &Integral;}_{{S S}_{n no}} {Λ Λ}_{m m} ((r r)) \cdot \cdot {Λ Λ}_{n no} (({r r}^{' '})) \frac{{\partial \partial}_{τ τ} {T T}_{j j} ((i i Δ Δ t t - - R R / / c c))}{R R} {ds ds}^{' '} d d s the s + + \\ \frac{11}{44 {πϵ πϵ}_{00}} {&Integral; &Integral;}_{{S S}_{m m}} {&Integral; &Integral;}_{{S S}_{n no}} &dtri; &dtri; \cdot &Center Dot; {Λ Λ}_{m m} ((r r)) &dtri; &dtri; \cdot &Center Dot; {Λ Λ}_{n no} (({r r}^{' '})) \frac{{\partial \partial}_{τ τ}^{- - 11} {T T}_{j j} ((i i Δ Δ t t - - R R / / c c))}{R R} {ds ds}^{' '} d d s the s \end{matrix} - - - - - - ((77))$

${[[{Z Z}_{M m}^{i i - - j j}]]}_{m m n no} = = \{\begin{matrix} \frac{11}{22} {&Integral; &Integral;}_{{S S}_{m m}} {Λ Λ}_{m m} ((r r)) \cdot &Center Dot; {Λ Λ}_{n no} (({r r}^{' '})) {T T}_{j j} ((i i Δ Δ t t)) d d s the s - - \frac{11}{44 π π} {&Integral; &Integral;}_{{S S}_{m m}} {Λ Λ}_{m m} ((r r)) \\ {{\overset{^^}{n no} ((r r)) \times \times {&Integral; &Integral;}_{{S S}_{n no}} \{\begin{matrix} {\partial \partial}_{τ τ} {T T}_{j j} ((i i Δ Δ t t - - R R / / c c)) {Λ Λ}_{m m} (({r r}^{' '})) \times \times \frac{R R}{{R R}^{22}} + + \\ {T T}_{j j} ((i i Δ Δ t t - - R R / / c c)) {Λ Λ}_{m m} (({r r}^{' '})) \times \times \frac{R R}{{R R}^{33}} \end{matrix}\}}} {ds ds}^{' '} d d s the s \end{matrix} - - - - - - ((88))$

第四步，远场部分采用泰勒级数展开重构成聚合、转移、配置的形式。 In the fourth step, the far-field part is reconstructed into the form of aggregation, transfer and configuration by using Taylor series expansion.

现在通过计算源点r'处的源信号J_n(r',t)对场点r处的辐射贡献来描述该算法的基本原理。假设源点r'所在的空间基函数为Λ_n(r')，因此r'处的源信号J_n(r',t)可以展开如下： The basic principle of the algorithm is now described by calculating the radiation contribution of the source signal J _n (r',t) at the source point r' to the radiation at the field point r. Assume that the spatial basis function where the source point r' is located is Λ _n (r'), so the source signal J _n (r',t) at r' can be expanded as follows:

${J J}_{n no} (({r r}^{' '},, t t)) = = {Σ Σ}_{j j = = 11}^{{N N}_{t t}} {I I}_{n no,, j j} {Λ Λ}_{n no} (({r r}^{' '})) {T T}_{j j} ((t t)) = = {Λ Λ}_{n no} (({r r}^{' '})) {g g}_{n no} ((t t)) - - - - - - ((99))$

实际操作中我们需要将源信号分解为持续时间更短的分段子信号，保证每段子信号的持续时间满足约束条件。源信号J_n(r',t)可以被分解为L段连续的子信号J_n,l(r',t)，每一段子信号的持续时间为T_s＝(M_t+1)Δt。源信号可以写成如下形式： In practice, we need to decompose the source signal into segmented sub-signals with shorter durations to ensure that the duration of each sub-signal satisfies the constraints. The source signal J _n (r',t) can be decomposed into L segments of continuous sub-signals J _n,l (r',t), and the duration of each segment of sub-signals is T _s =(M _t +1)Δt. The source signal can be written as follows:

${J J}_{n no} (({r r}^{' '},, t t)) = = {Σ Σ}_{l l = = 11}^{L L} {J J}_{n no,, l l} (({r r}^{' '},, t t)) = = {Σ Σ}_{l l = = 11}^{L L} {Λ Λ}_{n no} (({r r}^{' '})) {g g}_{n no,, l l} ((t t)) - - - - - - ((1010))$

那么源点r'处第l段子信号在场点r处产生的测试场为： Then the test field generated by the sub-signal of segment l at the source point r' at the field point r is:

$\begin{matrix} &lang; &lang; {Λ Λ}_{m m} ((r r)),, {n no}_{m m} \times \times {H h}_{n no,, l l} ((r r,, t t)) &rang; &rang; \\ = = \frac{11}{44 π π} {&Integral; &Integral;}_{{S S}_{m m}} {Λ Λ}_{m m} ((r r)) \cdot &Center Dot; {n no}_{m m} \times \times {{{&Integral; &Integral;}_{{S S}_{n no}} {dS wxya}^{' '} {{[[{\partial \partial}_{τ τ} {g g}_{n no,, l l} ((τ τ)) {Λ Λ}_{n no} (({r r}^{' '}))]] \times \times \frac{R R}{{cR c}^{22}} + + [[{g g}_{n no,, l l} ((τ τ)) {Λ Λ}_{n no} (({r r}^{' '}))]] \times \times \frac{R R}{{R R}^{33}}}}}} d d S S \\ \approx \approx {\frac{11}{44 π π}}_{m m} (({&Integral; &Integral;}_{{S S}_{m m}} {Λ Λ}_{m m} ((r r)) d d S S \cdot \cdot {n no}_{m m} \times \times {&Integral; &Integral;}_{{S S}_{n no}} {Λ Λ}_{n no} (({r r}^{' '})) d d {S S}^{' ' \times \times R R})) [[\frac{11}{{cR c}^{22}} {\partial \partial}_{τ τ} {g g}_{n no,, l l} ((τ τ)) + + \frac{11}{{R R}^{33}} {g g}_{n no,, l l} ((τ τ))]] \\ = = \frac{11}{44 π π} (({m m}_{m m} \cdot \cdot {n no}_{m m} \times \times {m m}_{n no} \times \times R R)) [[\frac{11}{{cR c}^{22}} {\partial \partial}_{τ τ} {g g}_{n no,, l l} ((τ τ)) + + \frac{11}{{R R}^{33}} {g g}_{n no,, l l} ((τ τ))]] \end{matrix} - - - - - - ((1111))$

$\begin{matrix} &lang; &lang; {Λ Λ}_{m m} ((r r)),, {E E.}_{n no,, l l} ((r r,, t t)) &rang; &rang; \\ = = {&Integral; &Integral;}_{{S S}_{m m}} {&Integral; &Integral;}_{{S S}_{n no}} {Λ Λ}_{m m} ((r r)) \cdot &Center Dot; {Λ Λ}_{n no} (({r r}^{' '})) \frac{{μ μ}_{00} {\partial \partial}_{τ τ} {g g}_{n no,, l l} ((τ τ))}{44 π π R R} {ds ds}^{' '} d d s the s + + {&Integral; &Integral;}_{{S S}_{m m}} {&Integral; &Integral;}_{{S S}_{n no}} &dtri; &dtri; \cdot \cdot {Λ Λ}_{m m} ((r r)) &dtri; &dtri; \cdot \cdot {Λ Λ}_{n no} (({r r}^{' '})) \frac{{\partial \partial}_{τ τ}^{- - 11} {g g}_{n no,, l l} ((τ τ))}{44 {πϵ πϵ}_{00} R R} {ds ds}^{' '} d d s the s \\ \approx \approx {&Integral; &Integral;}_{{S S}_{m m}} {Λ Λ}_{m m} ((r r)) d d s the s \cdot \cdot {&Integral; &Integral;}_{{S S}_{n no}} {Λ Λ}_{n no} (({r r}^{' '})) {ds ds}^{' '} \frac{{μ μ}_{00} {\partial \partial}_{τ τ} {g g}_{n no,, l l} ((τ τ))}{44 π π R R} + + {&Integral; &Integral;}_{{S S}_{m m}} &dtri; &dtri; \cdot &Center Dot; {Λ Λ}_{m m} ((r r)) d d s the s {&Integral; &Integral;}_{{S S}_{n no}} &dtri; &dtri; \cdot &Center Dot; {Λ Λ}_{n no} (({r r}^{' '})) {ds ds}^{' '} \frac{{\partial \partial}_{τ τ}^{- - 11} {g g}_{n no,, l l} ((τ τ))}{44 {πϵ πϵ}_{00} R R} \\ = = \frac{{μ μ}_{00}}{44 π π} {m m}_{m m} \cdot &Center Dot; {m m}_{n no} \frac{{\partial \partial}_{τ τ} {g g}_{n no,, l l} ((τ τ))}{R R} &PlusMinus; &PlusMinus; \frac{{l l}_{m m} {l l}_{n no}}{44 {πϵ πϵ}_{00}} \frac{{\partial \partial}_{τ τ}^{- - 11} {g g}_{n no,, l l} ((τ τ))}{R R} \end{matrix} - - - - - - ((1212))$

其中 $m_{m} = {&Integral;}_{S_{m}} Λ_{m} (r) d s, m_{n} = {&Integral;}_{S_{n}} Λ_{n} (r^{'}) {ds}^{'} .$ in $m_{m} = {&Integral;}_{S_{m}} Λ_{m} (r) d the s, m_{no} = {&Integral;}_{S_{no}} Λ_{no} (r^{'}) {ds}^{'} .$

现假设源点r_n与场点r_m分别位于两个组内，如图3所示，两个组分别称为源组和场组，组中心分别为r_i和r_j，场源基函数之间的矢量可以表示为： Now assume that the source point r _n and the field point r _m are located in two groups respectively, as shown in Fig. 3, the two groups are called the source group and the field group respectively, the centers of the groups are r _i and r _j respectively, and the field source basis function The vector between can be expressed as:

R＝r_mi+r_ij-r_nj＝R_m-R_n (13) R＝r _mi +r _ij -r _nj ＝R _m -R _n (13)

这里，r_ij＝r_i-r_j，r_mi＝r_m-r_i，r_nj＝r_n-r_j R_m＝r_mi+r_ij/2，R_n＝r_nj-r_ij/2 Here, r _ij =r _i -r _j , r _mi =r _m -r _i , r _nj =r _n -r _j R _m =r _mi +r _ij /2, R _n =r _nj -r _ij /2

利用泰勒级数展开可以得到如下表达式： Using Taylor series expansion, the following expression can be obtained:

$\begin{matrix} {R R}^{α α} = = {((R R \cdot \cdot R R))}^{\frac{α α}{22}} = = {[[(({r r}_{m m i i} + + {r r}_{i i j j} - - {r r}_{n no j j})) \cdot \cdot (({r r}_{m m i i} + + {r r}_{i i j j} - - {r r}_{n no j j}))]]}^{\frac{α α}{22}} \\ = = {r r}_{i i j j}^{α α} {[[11 + + ((\frac{22 {r r}_{m m i i} \cdot &Center Dot; {\overset{^^}{r r}}_{i i j j}}{{r r}_{i i j j}} + + \frac{{r r}_{m m i i}^{22}}{{r r}_{i i j j}^{22}})) + + ((\frac{22 {r r}_{n no j j} \cdot &Center Dot; {\overset{^^}{r r}}_{j j i i}}{{r r}_{i i j j}} + + \frac{{r r}_{n no j j}^{22}}{{r r}_{i i j j}^{22}})) - - \frac{22 {r r}_{m m i i} \cdot &Center Dot; {\overset{^^}{r r}}_{n no j j}}{{r r}_{i i j j}^{22}}]]}^{\frac{α α}{22}} \\ \approx \approx {R R}_{m m}^{((α α))} + + {R R}_{n no}^{((α α))} \end{matrix} - - - - - - ((1414))$

${R R}_{m m}^{((α α))} = = {r r}_{i i j j}^{α α} [[\frac{11}{22} + + α α ((\frac{{r r}_{m m i i} \cdot &Center Dot; {\overset{^^}{r r}}_{i i j j}}{{r r}_{i i j j}} + + \frac{{r r}_{m m i i}^{22} + + ((α α - - 22)) {(({r r}_{m m i i} \cdot \cdot {\overset{^^}{r r}}_{i i j j}))}^{22}}{22 {r r}_{i i j j}^{22}}))]] - - - - - - ((1515))$

${R R}_{n no}^{((α α))} = = {r r}_{i i j j}^{α α} [[\frac{11}{22} + + α α ((\frac{{r r}_{n no j j} \cdot \cdot {\overset{^^}{r r}}_{j j i i}}{{r r}_{i i j j}} + + \frac{{r r}_{n no j j}^{22} + + ((α α - - 22)) {(({r r}_{n no j j} \cdot &Center Dot; {\overset{^^}{r r}}_{j j i i}))}^{22}}{22 {r r}_{i i j j}^{22}}))]] - - - - - - ((1616))$

则式(11)(12)可以写成聚合、转移、配置的形式。 Then formulas (11) and (12) can be written in the forms of aggregation, transfer and configuration.

第五步，矩阵方程求解以及电磁散射参数的计算。 The fifth step is to solve the matrix equation and calculate the electromagnetic scattering parameters.

利用时间递推的方式求解每个时刻的电流系数。现将任意空间上的散射场按来源分为两部分:一部分由该场点所在组的近场组NFP(α)中的源产生的；另一部分由该场点所在组的远场组FFP(α)中的源产生。基于泰勒级数展开的时域积分方程快速算法的递推公式变为： The current coefficient at each moment is solved by time recursion. Now the scattered field in any space is divided into two parts according to the source: one part is generated by the source in the near-field group NFP(α) of the group where the field point is located; the other part is produced by the far-field group FFP(α) of the group where the field point is located. Source generation in α). The recursive formula of the fast algorithm for time-domain integral equations based on Taylor series expansion becomes:

${Z Z}^{00} {I I}^{i i} = = {V V}^{i i} - - \underset{{α α}^{' ' &Element; &Element; N N F f P P ((α α))}}{Σ Σ} {Σ Σ}_{j j = = 00}^{i i - - 11} {Z Z}_{i i - - j j}^{{αα αα}^{' '}} {I I}_{j j}^{{α α}^{' '}} - - \underset{{α α}^{' '} &Element; &Element; FF FF P P ((α α))}{Σ Σ} \underset{n no &Element; &Element; {α α}^{' '} ((n no))}{Σ Σ} &lang; &lang; {Λ Λ}_{m m} ((r r)),, {αE αE}_{n no,, l l} ((r r,, t t)) + + ((11 - - α α)) {n no}_{m m} \times \times {H h}_{n no,, l l} ((r r,, t t)) &rang; &rang; - - - - - - ((1717))$

近场区域产生的贡献是由经典MOT算法计算得到，主要通过矩阵元素值与电流系数相乘得到。在远场区域，也就是对应组中的远场组对，这些组之间的相互作用，过泰勒级数展开成聚合-转移-配置方法快速计算。 The contribution of the near-field region is calculated by the classical MOT algorithm, which is mainly obtained by multiplying the matrix element value and the current coefficient. In the far-field region, that is, the far-field pair in the corresponding group, the interaction between these groups is quickly calculated by the aggregation-transfer-configuration method through Taylor series expansion.

由于时刻iΔt以前的I^j在时刻iΔt都是已知的，j＝1,2,3,...,i-1，这样每个时间步求解一次式(17)的矩阵方程，就可以得到每个时刻iΔt上的Iⁱ所以可以递推求出各时刻的电流值。最后可以根据求得的瞬态电流系数计算出我们需要的电磁散射参数。 Since I ^j before time iΔt is known at time iΔt, j=1, 2, 3,..., i-1, so solving the matrix equation of formula (17) once at each time step, we can get The I ⁱ at each time iΔt can be recursively calculated to obtain the current value at each time. Finally, the electromagnetic scattering parameters we need can be calculated according to the obtained transient current coefficient.

为了验证本发明的正确性与有效性，下面分析了简易导弹模型的的电磁散射特性。 In order to verify the correctness and effectiveness of the present invention, the electromagnetic scattering characteristics of the simple missile model are analyzed below.

算例：简易导弹模型，几何尺寸为长8m，半径0.25m。激励源设置为：调制高斯脉冲，中心频率150MHz，带宽300MHz，入射波的方向θ＝180°,(弹头方向照射)，极化方式VV极化，观察角度0°≤θ≤180°,导弹模型采用0.1λ(λ＝1m)剖分得到8034个三角形，总未知量为12051。分组大小为0.3λ。3个盒子之外为远场区域，近场门限0.3λ*3。使用Δt＝333.3ps，一共计算了600个时间步。CFIE的参数α＝0.5，采用GMRES迭代方法计算，收敛精度为1e-9。图4给出了导弹模型在30MHz、150MHz以及270MHz时的双站RCS曲线并与经典MOT算法计算结果进行对比，表格1统计了本专利提出的方法与传统的时域积分方法消耗的时间和内存对比。 Calculation example: simple missile model, the geometric dimensions are 8m in length and 0.25m in radius. The excitation source is set as: modulated Gaussian pulse, center frequency 150MHz, bandwidth 300MHz, direction of incident wave θ=180°, (Irradiation in the direction of the warhead), polarization mode VV polarization, observation angle 0°≤θ≤180°, The missile model is divided into 8034 triangles by 0.1λ (λ=1m), and the total unknown quantity is 12051. The packet size is 0.3λ. Outside the 3 boxes is the far-field area, and the near-field threshold is 0.3λ*3. Using Δt = 333.3 ps, a total of 600 time steps were calculated. The parameter of CFIE is α=0.5, which is calculated by GMRES iterative method, and the convergence accuracy is 1e-9. Figure 4 shows the two-station RCS curves of the missile model at 30MHz, 150MHz and 270MHz and compares them with the calculation results of the classic MOT algorithm. Table 1 counts the time and memory consumed by the method proposed in this patent and the traditional time domain integration method Compared.

表1.计算时间和内存消耗对比 Table 1. Computational time and memory consumption comparison

从图4中可以看出本发明方法计算的双站RCS数据和MOT计算的结果十分吻合，证明了本发明方法的精确性。从表1看出，采用本发明方法相比经典MOT算法需要更好的计算内存和计算时间。进一步证明了本发明方法的有效性。 It can be seen from Fig. 4 that the bistation RCS data calculated by the method of the present invention is in good agreement with the result of the MOT calculation, which proves the accuracy of the method of the present invention. It can be seen from Table 1 that the method of the present invention requires better computing memory and computing time than the classical MOT algorithm. Further proves the effectiveness of the method of the present invention.

本发明的实现过程简单，相较于需要的复杂的时频域转换、谱域积分等操作，本发明方法只需在八叉树分组以及多层近远场划分的基础上，对远场组对之间的距离R进行泰勒展开重构成聚合、转移、投射的操作，编程相对简单易于实现。而且本发明方法适用性高，可以只需作较少的改动便可应用于求解介质、有耗、色散等复杂问题的分析中。 The implementation process of the present invention is simple. Compared with the required complex time-frequency domain conversion, spectral domain integration and other operations, the method of the present invention only needs to divide the far-field group on the basis of octree grouping and multi-layer near-far field division. The Taylor expansion of the distance R is reconstructed into aggregation, transfer, and projection operations, and the programming is relatively simple and easy to implement. Moreover, the method of the invention has high applicability, and can be applied to the analysis of solving complicated problems such as medium, lossy, dispersion, etc. only with few changes.

Claims

1. A fast algorithm for time-domain integral equations based on Taylor series expansion, characterized in that the steps are as follows:

The first step is to establish the time-domain electromagnetic field integral equation; use the incident electromagnetic field and the scattered electromagnetic field to meet the boundary conditions on the metal surface to establish the time-domain electromagnetic field integral equation;

The second step is to group the discretized sub-scatterers on the surface of the scatterer; surround the entire target object with a cube box, divide the cube into 8 sub-cubes, and then divide each sub-cube into 8 smaller cubes, By analogy, until the preset threshold value is reached, the division is stopped; the mutual coupling or autocoupling between any two sub-scatterers is divided into a near-field interaction pair and a far-field interaction pair according to their positional relationship; when they are near-field interaction When pairing, calculate the near-field impedance matrix, and when they are far-field acting pairs, use the Taylor series expansion to form the aggregation-transfer-configuration method to calculate the field generated by the source group basis function at the field group basis function;

The third step is to calculate the near-field impedance matrix, expand the metal surface current density with the space basis function and time basis function, and perform the Galerkin test in the space domain, and perform point matching in the time domain to obtain the matrix element values;

In the fourth step, the Taylor series expansion is used to reconstruct the far-field action pairs into the form of aggregation, transfer, and configuration to calculate the field generated by the source group basis function at the field group basis function;

The fifth step is to solve the matrix equation and calculate the electromagnetic scattering parameters; use the time recursion method to solve the current coefficient at each moment, use the iterative method to solve the final induced current coefficient, and finally calculate the The required electromagnetic scattering parameters.

2. the time domain integral equation fast algorithm based on Taylor series expansion according to claim 1, is characterized in that: in described step 4, on the basis of octree grouping and multi-layer near and far field division, to far Field action pairs are reconstructed by Taylor series expansion into aggregation, transfer, and projection operations to accelerate matrix-vector multiplication;

Assuming that the spatial basis function where the source point r' is located is Λ _n (r'), the source signal J _n (r',t) at r' is expanded as follows:

{J J}_{n no} (({r r}^{' '},, t t)) = = {Σ Σ}_{j j = = 11}^{{N N}_{t t}} {I I}_{n no,, j j} {Λ Λ}_{n no} (({r r}^{' '})) {T T}_{j j} ((t t)) = = {Λ Λ}_{n no} (({r r}^{' '})) {g g}_{n no} ((t t)) - - - - - - ((11))

Wherein the space basis function Λ _n (r) is the RWG basis function, the time basis function T _j (t) is the triangle basis function, and N _t is the number of the time basis function;

The source signal J _n (r',t) is decomposed into L segments of continuous sub-signals J _n,l (r',t), and the duration of each sub-signal is T _s =(M _t +1)Δt,M _t is the length of each sub-signal, and the source signal is written as follows:

{J J}_{n no} (({r r}^{' '},, t t)) = = {Σ Σ}_{l l = = 11}^{L L} {J J}_{n no,, l l} (({r r}^{' '},, t t)) = = {Σ Σ}_{l l = = 11}^{L L} {Λ Λ}_{n no} (({r r}^{' '})) {g g}_{n no,, l l} ((t t)) - - - - - - ((22))

The test electromagnetic field generated by the sub-signal of segment l at the source point r' at the field point r is:

\begin{matrix} &lang; &lang; {Λ Λ}_{m m} ((r r)),, {n no}_{m m} \times \times {H h}_{n no,, l l} ((r r,, t t)) &rang; &rang; \\ = = \frac{11}{44 π π} {&Integral; &Integral;}_{{S S}_{m m}} {Λ Λ}_{m m} ((r r)) \cdot &Center Dot; {n no}_{m m} \times \times {{{&Integral; &Integral;}_{{S S}_{n no}} {dS wxya}^{' '} {{[[{&PartialD; &PartialD;}_{τ τ} {g g}_{n no,, l l} ((τ τ)) {Λ Λ}_{n no} (({r r}^{' '}))]] \times \times \frac{R R}{{cR c}^{22}} + + [[{g g}_{n no,, l l} ((τ τ)) {Λ Λ}_{m m} (({r r}^{' '}))]] \times \times \frac{R R}{{R R}^{33}}}}}} dS wxya \\ \approx \approx \frac{11}{44 π π} (({&Integral; &Integral;}_{{S S}_{m m}} {Λ Λ}_{m m} ((r r)) dS wxya \cdot &Center Dot; {n no}_{m m} \times \times {&Integral; &Integral;}_{{S S}_{n no}} {Λ Λ}_{n no} (({r r}^{' '})) {dS wxya}^{' '} \times \times R R)) [[\frac{11}{c c {R R}^{22}} {&PartialD; &PartialD;}_{τ τ} {g g}_{n no,, l l} ((τ τ)) + + \frac{11}{{R R}^{33}} {g g}_{n no,, l l} ((τ τ))]] \\ = = \frac{11}{44 π π} (({m m}_{m m} \cdot &Center Dot; {n no}_{m m} \times \times {m m}_{n no} \times \times R R)) [[\frac{11}{c c {R R}^{22}} {&PartialD; &PartialD;}_{τ τ} {g g}_{n no,, l l} ((τ τ)) + + \frac{11}{{R R}^{33}} {g g}_{n no,, l l} ((τ τ))]] \end{matrix} - - - - - - ((33))

\begin{matrix} &lang; &lang; {Λ Λ}_{m m} ((r r)),, {E E.}_{n no,, l l} ((r r,, t t)) &rang; &rang; \\ = = {&Integral; &Integral;}_{{S S}_{m m}} {&Integral; &Integral;}_{{S S}_{n no}} {Λ Λ}_{m m} ((r r)) \cdot &Center Dot; {Λ Λ}_{n no} (({r r}^{' '})) \frac{μ μ {&PartialD; &PartialD;}_{τ τ} {g g}_{n no,, l l} ((τ τ))}{44 πR πR} {ds ds}^{' '} ds ds + + {&Integral; &Integral;}_{{S S}_{m m}} {&Integral; &Integral;}_{{S S}_{n no}} &dtri; &dtri; \cdot \cdot {Λ Λ}_{m m} ((r r)) &dtri; &dtri; \cdot \cdot {Λ Λ}_{n no} (({r r}^{' '})) \frac{{&PartialD; &PartialD;}_{τ τ}^{- - 11} {g g}_{n no,, l l} ((τ τ))}{44 πϵR πϵR} {ds ds}^{' '} ds ds \\ \approx \approx {&Integral; &Integral;}_{{S S}_{m m}} {Λ Λ}_{m m} ((r r)) ds ds \cdot \cdot {&Integral; &Integral;}_{{S S}_{n no}} {Λ Λ}_{n no} (({r r}^{' '})) {ds ds}^{' '} \frac{μ μ {&PartialD; &PartialD;}_{τ τ} {g g}_{n no,, l l} ((τ τ))}{44 πR πR} + + {&Integral; &Integral;}_{{S S}_{m m}} &dtri; &dtri; \cdot \cdot {Λ Λ}_{m m} ((r r)) ds ds {&Integral; &Integral;}_{{S S}_{n no}} &dtri; &dtri; \cdot &Center Dot; {Λ Λ}_{n no} (({r r}^{' '})) {ds ds}^{' '} \frac{{&PartialD; &PartialD;}_{τ τ}^{- - 11} {g g}_{n no,, l l} ((τ τ))}{44 πϵR πϵR} \\ = = \frac{μ μ}{44 π π} {m m}_{m m} \cdot &Center Dot; {m m}_{n no} \frac{{&PartialD; &PartialD;}_{τ τ} {g g}_{n no,, l l} ((τ τ))}{R R} &PlusMinus; &PlusMinus; \frac{{l l}_{m m} {l l}_{n no}}{44 πϵ πϵ} \frac{{&PartialD; &PartialD;}_{τ τ}^{- - 11} {g g}_{n no,, l l} ((τ τ))}{R R} \end{matrix} - - - - - - ((44))

where Λ _m (r) is the test basis function at field point r, n _m is the unit external normal vector at field point r, E _n,l (r,t),H _n,l (r,t) are The electromagnetic and magnetic field generated by the lth sub-signal at the source point r' at the field point r, μ and ε are the permeability and permittivity of the free space respectively; R=|r-r'|, c is the The speed of light, τ=tR/c is the time delay; l _m , l _n are the side lengths of the m and n sides respectively, is the integral of the basis function in its support domain;

The source point r _n and the field point r _m are respectively located in two groups, the two groups are called the source group and the field group respectively, and the centers of the groups are r _i and r _j respectively, and the vector between the field source basis functions is expressed as:

R＝r _mi +r _ij -r _nj ＝R _m -R _n (5)

Here, r _ij =r _i -r _j , r _mi =r _m -r _i , r _nj =r _n -r _j R _m =r _mi +r _ij /2, R _n =r _nj -r _ij /2

Using Taylor series expansion, the following expression is obtained:

\begin{matrix} {R R}^{α α} = = {((R R \cdot \cdot R R))}^{\frac{α α}{22}} = = {[[(({r r}_{mi mi} + + {r r}_{ij ij} - - {r r}_{nj nj})) \cdot \cdot (({r r}_{mi mi} + + {r r}_{ij ij} - - {r r}_{ni ni}))]]}^{\frac{α α}{22}} \\ = = {r r}_{ij ij}^{α α} {[[11 + + ((\frac{{22 r r}_{mi mi} \cdot \cdot {\overset{^^}{r r}}_{ij ij}}{{r r}_{ij ij}} + + \frac{{r r}_{mi mi}^{22}}{{r r}_{ij ij}^{22}})) + + ((\frac{{22 r r}_{nj nj} \cdot &Center Dot; {\overset{^^}{r r}}_{ji the ji}}{{r r}_{ij ij}} + + \frac{{r r}_{nj nj}^{22}}{{r r}_{ij ij}^{22}})) \frac{22 {r r}_{mi mi} \cdot &Center Dot; {r r}_{nj nj}}{{r r}_{ij ij}^{22}}]]}^{\frac{α α}{22}} \\ \approx \approx {R R}_{m m}^{((α α))} + + {R R}_{n no}^{((α α))} \end{matrix} - - - - - - ((66))

{R R}_{m m}^{((α α))} = = {r r}_{ij ij}^{α α} [[\frac{11}{22} + + α α ((\frac{{r r}_{mi mi} \cdot &Center Dot; {\overset{^^}{r r}}_{ij ij}}{{r r}_{ij ij}} + + \frac{{r r}_{mi mi}^{22} + + ((α α - - 22)) {(({r r}_{mi mi} \cdot &Center Dot; {\overset{^^}{r r}}_{ij ij}))}^{22}}{22 {r r}_{ij ij}^{22}}))]] - - - - - - ((77))

{R R}_{n no}^{((α α))} = = {r r}_{ij ij}^{α α} [[\frac{11}{22} + + α α ((\frac{{r r}_{nj nj} \cdot &Center Dot; {\overset{^^}{r r}}_{ji the ji}}{{r r}_{ij ij}} + + \frac{{r r}_{nj nj}^{22} + + ((α α - - 22)) {(({r r}_{nj nj} \cdot &Center Dot; {\overset{^^}{r r}}_{ji the ji}))}^{22}}{22 {r r}_{ij ij}^{22}}))]] - - - - - - ((88))

Substituting (5-8) into formula (3) and (4) will be written in the form of aggregation, transfer, and configuration to accelerate matrix-vector multiplication.