<figure><img src="image_1.png" /></figure>

HAL
open science

# Surrogate modeling of time-domain electromagnetic wave propagation via dynamic mode decomposition and radial basis function

Kun Li, Yixin Li, Liang Li, Stéphane Lanteri

► To cite this version:

Kun Li, Yixin Li, Liang Li, Stéphane Lanteri. Surrogate modeling of time-domain electromagnetic wave propagation via dynamic mode decomposition and radial basis function. Journal of Computational Physics, 2023, 491, pp.112354. 10.1016/j.jcp.2023.112354. hal-04401269

HAL Id: hal-04401269

https://inria.hal.science/hal-04401269v1

Submitted on 17 Jan 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Surrogate modeling of time-domain electromagnetic wave propagation
via dynamic mode decomposition and radial basis function

Kun Li<sup>a</sup>, Yixin Li<sup>a</sup>, Liang Li<sup>b,*</sup>, Stéphane Lanteri<sup>c</sup>

<sup>a</sup>School of Mathematics, Southwestern University of Finance and Economics, Chengdu, P.R. China

<sup>b</sup>School of Mathematical Sciences, University of Electronic Science and Technology of China, Chengdu, P.R. China

<sup>c</sup>Université Côte d'Azur, Inria, CNRS, LJAD, Sophia Antipolis, France

**Abstract**

This work introduces an 'equation-free' non-intrusive model order reduction (NIMOR) method for surrogate modeling of time-domain electromagnetic wave propagation. The nested proper orthogonal decomposition (POD) method, as a prior dimensionality reduction technique, is employed to extract the time- and parameter-independent reduced basis (RB) functions from a collection of high-fidelity (HF) solutions (or snapshots) on a properly chosen training parameter set. A dynamic mode decomposition (DMD) method, resulting in a further dimension reduction of the NIMOR method, is then used to predict the reduced-order coefficient vectors for future time instants on the previous training parameter set. The radial basis function (RBF) is employed for approximating the reduced-order coefficient vectors at a given untrained parameter in the future time instants, leading to the applicability of DMD method to parameterized problems. A main advantage of the proposed method is the use of a multi-step procedure consisting of the POD, DMD and RBF techniques to accurately and quickly recover field solutions from a few large-scale simulations. Numerical experiments for the scattering of a plane wave by a dielectric disk, by a multi-layer disk, and by a 3-D dielectric sphere nicely illustrate the performance of the NIMOR method.

**Keywords:** Surrogate modeling; non-intrusive model order reduction; proper orthogonal decomposition; dynamic mode decomposition; radial basis function

**1. Introduction**

Numerical simulations of electromagnetic wave propagation problems primarily rely on spa-
tial and temporal discretizations of the system of time-domain Maxwell equations using finite
difference or finite element type methods. For complex and realistic 3-D situations, such a pro-
cess can be computationally prohibitive, especially in view of many-query analyses (e.g., design,
control, optimization, and uncertainty quantification (UQ) [1, 2, 3, 4]). Therefore, developing
cost-effective and high-precision surrogate models are of great practical significance. Among
the different possible approaches for building the surrogate model of a given partial differential
equation (PDE) system, the data-driven based *non-intrusive model order reduction* (NIMOR)
approaches, which are generally implemented in a fully decoupled *offline-online* paradigm [5, 6],

*Corresponding author
Email addresses: rclikun@163.com, likun@swufe.edu.cn (Kun Li), yxli19900@163.com (Yixin Li), plum_liliang@uestc.edu.cn, plum.liliang@gmail.com (Liang Li), Stephane.Lanteri@inria.fr (Stéphane Lanteri)
have recently shown new promises due to their capability of handling nonlinear or/and high dimensional problems.

During the offline stage, a set of time- and parameter-independent reduced basis (RB) func-
tions are extracted from a collection of high-fidelity (HF) solutions (also known as snapshots),
which are always generated by numerical simulations or experiments at some time/parameter lo-
cations. A reduced space, which is spaned by the RB functions, can effectively capture the main
dynamics of the original full-order model (FOM). A proper orthogonal decomposition (POD)
method [7, 8] based on a singular value decomposition (SVD) is usually used for extracting
the RB functions. However, for time-domain electromagnetic wave propagation problems, the
snapshot matrix has a considerable size, leading to very a time-consuming process if the POD
method is used directly. To overcome this difficulty, we use a nested POD method [9], in which
the POD method is implemented separately for time and parameter spaces. As an approxima-
tion to the full-order solutions, the reduced-order solutions are recovered as a linear combination
of the RB functions, where only the combination coefficients, also named the reduced-order or
projection coefficients will be calculated during the online stage. In the NIMOR method, the
reduced-order coefficients are approximated over the parameter domain by only using a database
of reduced-order information [10].

More recently, NIMOR methods (whose mainly advantage is that they are 'equation-free'
and do not require access to the governing equations of original FOM) have been developed
for solving large-scale complex physical problems. For instance, Hesthaven and Ubbiali [5]
proposed a non-intrusive RB method combining the POD method with feed forward neural
networks (FNNs) for parameterized steady-state PDEs, which is extended to time-dependent
nonlinear problems [6], in particular to unsteady flows. Chen and et al. [11] proposed a
physics-reinforced neural network (PRNN) for nonlinear problems, which is a combination of
projection driven neural network (PDNN) with physics-informed neural network (PINN) [12].
Fu and et al. [13] developed a physics-data combined machine learning (PDCML) method for
parametric reduced-order modelling of nonlinear dynamical systems in small-data regimes. In
addition, different regression and interpolation methods, such as Gaussian process regression
(GPR) [14, 15, 16], radial basis function (RBF) [17, 18] and cubic spline interpolation (CSI)
[9, 19], have also been adopted to approximate the mapping from the time/parameter to the
reduced-order coefficients.

An alternative approach to construct the non-intrusive reduced-order model (ROM) is the *dynamic mode decomposition* (DMD) method, which is also an 'equation-free' method. The DMD method was developed by Schmid [20] and has gained popularity in last years as ROM for dynamical systems [21]. The accuracy of the DMD method has been proved numerically [22] and theoretically [23]. Moreover, many variants of the DMD method have also been proposed to improve their efficiency and applicability [24, 25, 26]. For instance, Clainchey and et al. proposed a higher order DMD (HODMD) method [24] that uses time-lagged snapshots and can be seen as superimposed DMD in a sliding window, which is also applied successfully to identify and extrapolate flow patterns [27]. Despite the DMD method is wildly used for reduced-order modelling of time-dependent problems, it is significantly more challenging for parameterized problems [28, 29, 30, 31]. For solving this limitation, Syadi, Schmid and et al. proposed a DMD-based parameter-dependent ROM framework for bifurcation analysis [32], where the time series HF solutions for different parameter values are "stacked" to form an augmented snapshot matrix, and the stacked DMD modes, which are calculated by the augmented snapshot matrix, are interpolated to form a new DMD mode for a new parameter. Further, Huhn and et al. developed two computationally more efficient parametric DMD approaches for diffusion-dominated nonlinear dynamical problems, including the reduced eigen-pair interpolation (rEPI) and the reduced Koopman operator interpolation (rKOI). In addition, the DMD methods are
also employed to build the prediction models of snapshots or latent variables at each training parameter value, and then some interpolation techniques are adopted to recover the snapshots or latent variables for some new parameter values [28, 33].

In this study, we attempt to investigate the predictive ability of a higher-order extension of the DMD method [24], i.e. HODMD, for time-domain electromagnetic wave propagation from both interpolation and extrapolation perspectives. In order to extend the DMD method to parametrized problems, we propose the RBF method to interpolate the reduced-order coefficient vectors for some new parameters. The proposed NIMOR method is purely data-driven method, where the offline and online stages are fully decoupled. During the offline stage, the HF solutions as snapshots are collected, the RB bases are extracted, and the DMD models of reduced-order coefficient vectors for all training parameter values are established. In order to predict a RB solution at a new time and parameter value on the online stage, the reduced-order coefficient vector at the new time is first predicted via the surrogate models formed by the DMD models, then the new reduced-order coefficient vector in the parameter space is interpolated via the RBF method, the RB solution is finally recovered in reduced space formed by the RB bases. The main innovation of this work is the use of a multi-step procedure consisting of the nested POD, DMD and RBF methods to accurately recover field solutions from a small number of large-scale electromagnetic simulations.

The remainder of this paper is organized as follows. We briefly introduce the time-domain Maxwell's equations and the training data set in Section 2. Section 3 describes the NIMOR method including the nested POD, DMD, and RBF methods. Section 4 presents and discusses some numerical results, including the scattering of a plane wave by a dielectric disk, by a multi-layer disk, and by a 3-D dielectric sphere. Section 5 concludes this study.

## 2. Full-order solutions

### 2.1. Time-domain Maxwell's equations

We consider the electromagnetic scattering problems governed by the normalized time-domain Maxwell's equations with the Silver-Müller absorbing boundary condition (ABC) and the initial conditions

$$ 
\left\{
\begin{array}{ll}
\mu_r \frac{\partial \mathbf{H}(\mathbf{x}, t)}{\partial t} + \nabla \times \mathbf{E}(\mathbf{x}, t) = 0, & \forall (\mathbf{x}, t) \in \Omega \times \mathcal{T}, \\
\varepsilon_r \frac{\partial \mathbf{E}(\mathbf{x}, t)}{\partial t} - \nabla \times \mathbf{H}(\mathbf{x}, t) = 0, & \forall (\mathbf{x}, t) \in \Omega \times \mathcal{T}, \\
\mathcal{L}(\mathbf{E}(\mathbf{x}, t), \mathbf{H}(\mathbf{x}, t)) = \mathcal{L}(\mathbf{E}^{\text{inc}}(\mathbf{x}, t), \mathbf{H}^{\text{inc}}(\mathbf{x}, t)), & \forall (\mathbf{x}, t) \in \partial\Omega \times \mathcal{T}, \\
\mathbf{E}(\mathbf{x}, 0) = \mathbf{E}_0(\mathbf{x}), \quad \mathbf{H}(\mathbf{x}, 0) = \mathbf{H}_0(\mathbf{x}), & \forall \mathbf{x} \in \Omega,
\end{array}
\right.
\quad (1) $$

where $\Omega \subset \mathbb{R}^2$ or $\mathbb{R}^3$ is the spatial domain, $\mathcal{T} = [0, T_f]$ is the time interval, and $\partial\Omega$ is the boundary of $\Omega$; $\mathbf{E}(\mathbf{x}, t)$ and $\mathbf{H}(\mathbf{x}, t)$ respectively denote the electric and magnetic fields; $\varepsilon_r$ and $\mu_r$ are respectively the relative permittivity and permeability; the Silver-Müller ABC operator is $\mathcal{L}(\mathbf{E}(\mathbf{x}, t), \mathbf{H}(\mathbf{x}, t)) = \mathbf{n} \times \mathbf{E}(\mathbf{x}, t) + Z\mathbf{n} \times (\mathbf{n} \times \mathbf{H}(\mathbf{x}, t))$, $\mathbf{n}$ denotes the unit normal vector pointing outward to $\partial\Omega$, $\mathbf{E}^{\text{inc}}(\mathbf{x}, t)$ and $\mathbf{H}^{\text{inc}}(\mathbf{x}, t)$ stand for the incident fields, and $Z = \sqrt{\mu_r/\varepsilon_r}$; $\mathbf{E}_0$ and $\mathbf{H}_0$ stand for some given functions. In this study, we consider $\theta = (\varepsilon_{r,1}, \varepsilon_{r,2}, \dots, \varepsilon_{r,p}) \in \mathcal{P} \subset \mathbb{R}^p$ as the problem's parameters with $\varepsilon_{r,i}$ ($i = 1, 2, \dots, p$) being the relative permittivity in the $i$-th domain of $\Omega$, $\mathcal{P}$ being the parameter domain, and $p$ being the number of parameters.
## 2.2. The data set of the high-fidelity solutions

Let $\mathcal{P}_h^{tr} = \{\theta^{(1)}, \theta^{(2)}, \dots, \theta^{(N_p)}\} \subset \mathcal{P}$ be the training parameter set with uniform or uneven sampling. For each parameter $\theta^{(i)} \in \mathcal{P}_h^{tr}$, the corresponding high-fidelity (HF) solutions or snapshots for (1) can be obtained by using classical numerical methods or physical experiments in a time partition $\mathcal{T}_h = \{t^{(0)}, t^{(1)}, \dots, t^{(m)}\}$. Then, we can formulate the following time snapshot matrices $\mathcal{S}_u^i$ ($i = 1, 2, \dots, N_p$) in a time sampling $\mathcal{T}_h^{tr} = \{t^{(n_1)}, t^{(n_2)}, \dots, t^{(n_{N_t})}\} \subseteq \mathcal{T}_h$ for each parameter sampling point

$$ \mathcal{S}_u^i = \begin{pmatrix} \mathbf{u}_{h,1}^{(n_1)}(\theta^{(i)}) - \hat{\mathbf{u}}_{h,1} & \mathbf{u}_{h,1}^{(n_2)}(\theta^{(i)}) - \hat{\mathbf{u}}_{h,1} & \cdots & \mathbf{u}_{h,1}^{(n_{N_t})}(\theta^{(i)}) - \hat{\mathbf{u}}_{h,1} \\ \mathbf{u}_{h,2}^{(n_1)}(\theta^{(i)}) - \hat{\mathbf{u}}_{h,2} & \mathbf{E}_{h,2}^{(n_2)}(\theta^{(i)}) - \hat{\mathbf{u}}_{h,2} & \cdots & \mathbf{u}_{h,2}^{(n_{N_t})}(\theta^{(i)}) - \hat{\mathbf{u}}_{h,2} \\ \vdots & \vdots & \ddots & \vdots \\ \mathbf{u}_{h,N_h}^{(n_1)}(\theta^{(i)}) - \hat{\mathbf{u}}_{h,N_h} & \mathbf{u}_{h,N_h}^{(n_2)}(\theta^{(i)}) - \hat{\mathbf{u}}_{h,N_h} & \cdots & \mathbf{u}_{h,N_h}^{(n_{N_t})}(\theta^{(i)}) - \hat{\mathbf{u}}_{h,N_h} \end{pmatrix} \in \mathcal{R}^{N_h \times N_t}, \quad (2) $$

and then the snapshot matrix for all parameter in $\mathcal{P}_h^{tr}$ is

$$ \tilde{\mathcal{S}}_{\mathbf{u}} = (\mathcal{S}_{\mathbf{u}}^{1} | \mathcal{S}_{\mathbf{u}}^{2} | \cdots | \mathcal{S}_{\mathbf{u}}^{N_p}) \in \mathcal{R}^{N_h \times N_{tp}}, \mathbf{u} \in \{\mathbf{E}, \mathbf{H}\}, \quad (3) $$

where $\mathbf{u}_{h,k}^{(n_j)}(\theta^{(i)}) = \mathbf{u}_{h,k}(t^{(n_j)}, \theta^{(i)})$ is the HF solution, $\hat{\mathbf{u}}_h \in \mathcal{R}^{N_h}$ is the mean value, $N_h$ is the number of degrees of freedom (DoF) in space, and $N_{tp} = N_t \cdot N_p$ is the number of all snapshots.

# 3. Parameterized non-intrusive model order reduction method

In this section, we consider a possible application of the POD, DMD and RBF methods in the parameterized non-intrusive reduced-order modeling.

## 3.1. Nested proper orthogonal decomposition method

Based on the above data preparation, we perform a low-rank approximation by constructing a low-dimension reduced basis (RB) subspace $\mathcal{V}_{rb,u}$ with rank $L_u \ll \min\{N_h, N_{tp}\}$ ($\mathbf{u} \in \{\mathbf{E}, \mathbf{H}\}$), which has the lowest dimension possible, while still maintaining a good approximation of the original dataset. The proper orthogonal decomposition (POD) method with snapshots is a popular approach to form the RB subspace. However, for the parameterized electromagnetic scattering problems, the snapshot matrix $\tilde{\mathcal{S}}_u$ ($\mathbf{u} \in \{\mathbf{E}, \mathbf{H}\}$) has a considerable size because time is added as an input parameter, so a direct approach of a singular value decomposition (SVD) is very time consuming. Therefore, a nested POD method [9] is employed to extract the time- and parameter-independent RB functions, where the POD method is implemented separately for time and parameter spaces. The detailed pseudo-code implementation of the nested POD method is shown in Algorithm 1. Figure 1 shows the framework of the classic POD and nested POD. For a study of the corresponding projection errors we refer to [9].

Let $\Psi_u = [\psi_{u,1}, \psi_{u,2}, \dots, \psi_{u,L_u}] \in \mathcal{R}^{N_h \times L_u}$ ($\mathbf{u} \in \{\mathbf{E}, \mathbf{H}\}$) be the RB matrix. According to the projection theory, the RB solution $\mathbf{u}_h^r(t, \theta)$, an approximation to the full-order solution $\mathbf{u}_h(t, \theta)$, can be recovered as

$$ \mathbf{u}_h(t, \theta) \approx \mathbf{u}_h^r(t, \theta) = \Psi_u \alpha_u(t, \theta) + \hat{\mathbf{u}}_h = \sum_{i=1}^{L_u} \alpha_{u,i}(t, \theta) \phi_{u,i} + \hat{\mathbf{u}}_h, \quad \mathbf{u} \in \{\mathbf{E}, \mathbf{H}\}, \quad (4) $$

where $\alpha_u(t, \theta) = [\alpha_{u,1}(t, \theta), \alpha_{u,2}(t, \theta), \dots, \alpha_{u,L_u}(t, \theta)]^T \in \mathcal{R}^{L_u}$ denotes the reduced-order or projection coefficient vector. In particular, the following reduced-order coefficient vectors can be calculated approximately in the training time and parameter sets

$$ \alpha_{\mathbf{u}}(t^{(n_i)}, \theta^{(j)}) = \Psi_{\mathbf{u}}^{T}(\mathbf{u}_h(t^{(n_i)}, \theta^{(j)}) - \hat{\mathbf{u}}_h), \quad (t^{(n_i)}, \theta^{(j)}) \in \mathcal{T}_h^{tr} \times \mathcal{P}_h^{tr}, \mathbf{u} \in \{\mathbf{E}, \mathbf{H}\}, \quad (5) $$

which can be used as the training data.
<caption>**Algorithm 1 Nested POD method**</caption>

**Input:** Time trajectory matrices $\mathbb{S}_u^i$ ($i = 1, 2, \dots, N_p$, $u \in \{\mathbf{E}, \mathbf{H}\}$), and truncation tolerances $\rho_t$ and $\rho_\theta$.
**Output:** POD basis matrix $\Psi_u$ ($u \in \{\mathbf{E}, \mathbf{H}\}$)

1: **for** $i = 1$ to $N_p$ **do**

2: $\Psi_u^i = \text{POD}(\mathbb{S}_u^i, \rho_t)$ with the size $l_u^i \ll \min\{N_h, N_t\}$

3: **end for**

4: $\hat{\mathbb{S}}_u = [\Psi_u^1, \Psi_u^2, \dots, \Psi_u^{N_p}]$

5: $\Psi_u = \text{POD}(\hat{\mathbb{S}}_u, \rho_\theta)$

6: **function** $\Psi = \text{POD}(\mathbb{S}, \rho)$ **do**

7: $[\mathbb{U}, \Sigma, \mathbb{V}] = \text{svd}(\mathbb{S})$

8: $k = \text{argmin}\{\mathcal{E}(k) : \mathcal{E}(k) \ge 1 - \rho\}$ with $\mathcal{E}(k) = \sum_{i=1}^k \sigma_i^2 / \sum_{i=1}^r \sigma_i^2$ being the relative information content, where $r$ is the rank of $\mathbb{S}$ and $\sigma_i$ is the singular value of $\mathbb{S}$

9: $\Psi = \mathbb{U}(:, 1:k)$

10: **end function**

<figure><img src="image_3.png" /><figcaption>Figure 1: Classic POD and nested POD.</figcaption></figure>
**Remark 1.** In Algorithm 1, the effects of $(\rho_t, \rho_\theta)$ on the average projection and NIMOR errors have been investigated to provide guidance on the choice of truncation tolerances in our previous work [9]. So we directly select the truncation tolerances and no longer test the numerical convergences in this study.

## 3.2. Dynamic mode decomposition method

The dynamic mode decomposition (DMD) method with data-driven framework is proposed for model reduction of time-dependent problems. Based on the database of the reduced-order coefficient vectors over a time domain for each parameter, the DMD is defined as an approximate eigendecomposition of the operator that relates one time step to the next. For each training parameter $\theta^{(j)} \in \mathcal{P}_h^{tr}$, we seek an operator (or a Koopman matrix) $\mathcal{K}_u^{(j)} \in \mathbb{R}^{L_u \times L_u}$ such that

$$ \alpha_u(t^{(n_{i+1})}, \theta^{(j)}) = \mathcal{K}_u^{(j)} \alpha_u(t^{(n_i)}, \theta^{(j)}), \quad i = 1, 2, \dots, N_t - 1, \mathbf{u} \in \{\mathbf{E}, \mathbf{H}\}, \quad (6) $$

where the operator $\mathcal{K}_u^{(j)}$ can be considered as a finite approximation of the Koopman operator [34], which is infinite-dimensional. In order to establish the operator $\mathcal{K}_u^{(j)}$ satisfying (6) to an acceptable approximation, we consider the time sampling $\mathcal{T}_h^{tr}$ with equal interval $\Delta t$. In this study, the snapshots are sampled at the first 70% of a period with equal interval $\Delta t$ for each training parameter point. In order to compute the operator $\mathcal{K}_u^{(j)}$, we arrange two following sequences (whose columns are the reduced-order coefficient vectors)

$$ \mathcal{X}_{\alpha_u}^{(j)} = [\alpha_u(t^{(n_1)}, \theta^{(j)}), \alpha_u(t^{(n_2)}, \theta^{(j)}), \dots, \alpha_u(t^{(n_{N_t-1})}, \theta^{(j)})] \in \mathbb{R}^{L_u \times (N_t-1)}, \quad (7) $$

and

$$ \mathcal{Y}_{\alpha_u}^{(j)} = [\alpha_u(t^{(n_2)}, \theta^{(j)}), \alpha_u(t^{(n_3)}, \theta^{(j)}), \dots, \alpha_u(t^{(n_{N_t})}, \theta^{(j)})] \in \mathbb{R}^{L_u \times (N_t-1)}. \quad (8) $$

One can assume that $\mathcal{X}_{\alpha_u}^{(j)}$ and $\mathcal{Y}_{\alpha_u}^{(j)}$ are connected by the following linear dynamics

$$ \mathcal{Y}_{\alpha_u}^{(j)} \approx \mathcal{K}_u^{(j)} \mathcal{X}_{\alpha_u}^{(j)}, \quad j = 1, 2, \dots, N_p, \mathbf{u} \in \{\mathbf{E}, \mathbf{H}\}. \quad (9) $$

Let $\dagger$ be the Moore-Penrose pseudoinverse, one can find the matrix $\mathcal{K}_u^{(j)} = \mathcal{Y}_{\alpha_u}^{(j)}(\mathcal{X}_{\alpha_u}^{(j)})^{\dagger}$ is a best-fit linear operator relating $\mathcal{X}_{\alpha_u}$ to $\mathcal{Y}_{\alpha_u}$ [35]. However, this calculation is expensive and often difficult to handle [36]. Therefore we adopt another method, in which the DMD can be implemented by the Algorithm 2 for the approximation of the eigenvalues and eigenvectors of $\mathcal{K}_u^{(j)}$ based on the data sets $\mathcal{X}_{\alpha_u}^{(j)}$ and $\mathcal{Y}_{\alpha_u}^{(j)}$.

Performing Algorithm 2, the eigenvalues and the corresponding eigenvectors of $\mathcal{K}_u^{(j)}$ can be approximated by $\Lambda_u^{(j)}$ and the columns of $\Phi_u^{(j)}$. Then, one can predict the solution for a new time $t^*$ in the future as follows

$$ \alpha_u(t^*, \theta^{(j)}) \approx \alpha_u^{\text{DMD}}(t^*, \theta^{(j)}) \\ = \Phi_u^{(j)}(\Lambda_u^{(j)})(t^*-t^{(n_1)})/\Delta t (\Phi_u^{(j)})^{\dagger} \alpha_u(t^{(n_1)}, \theta^{(j)}), \quad j = 1, 2, \dots, N_p, \mathbf{u} \in \{\mathbf{E}, \mathbf{H}\}. \quad (11) $$

In particular, if $L_u \approx N_t$, the quality of the results obtained by the DMD method may decay [24]; the situation could be worse if $L_u < N_t$ [36]. So, to avoid this situation, we consider a higher order DMD (termed as HODMD or DMD-s) method [24], which is able to extract more significative modes from low-dimensional data. In the HODMD method, (6) is modified to the following form

$$ \alpha_u(t^{(n_{i+s})}, \theta^{(j)}) = \mathcal{K}_{u,1}^{(j)} \alpha_u(t^{(n_i)}, \theta^{(j)}) + \mathcal{K}_{u,2}^{(j)} \alpha_u(t^{(n_{i+1})}, \theta^{(j)}) + \dots + \mathcal{K}_{u,s}^{(j)} \alpha_u(t^{(n_{i+s-1})}, \theta^{(j)}), \quad (12) $$
**Algorithm 2** DMD method

**Input:** Two matrices $\mathcal{X}_{\alpha_u}^{(j)}$, and $\mathcal{Y}_{\alpha_u}^{(j)}$, $j = 1, 2, \dots, \mathcal{N}_p$, $\mathbf{u} \in \{\mathbf{E}, \mathbf{H}\}$

**Output:** DMD modes $\Phi_{\mathbf{u}}^{(j)}$ and eigenvalue matrices $\Lambda_{\mathbf{u}}^{(j)}$, $j = 1, 2, \dots, \mathcal{N}_p$, $\mathbf{u} \in \{\mathbf{E}, \mathbf{H}\}$

1. Perform the SVD of $\mathcal{X}_{\alpha_u}^{(j)} \approx \mathcal{U}_{\mathbf{u},r_j}^{(j)} \Sigma_{\mathbf{u},r_j}^{(j)} (\mathcal{V}_{\mathbf{u},r_j}^{(j)})^T$ with $r_j$ being the rank of $\mathcal{X}_{\alpha_u}^{(j)}$

2. Calculate the low-rank projection $\tilde{\mathcal{K}}_{\mathbf{u}}^{(j)}$ (or called the reduced Koopman matrix) of $\mathcal{K}_{\mathbf{u}}^{(j)}$ onto the subspace spanned by $\mathcal{U}_{\mathbf{u},r_j}^{(j)}$:

$$ \tilde{\mathcal{K}}_{\mathbf{u}}^{(j)} = (\mathcal{U}_{\mathbf{u},r_j}^{(j)})^T \mathcal{K}_{\mathbf{u}}^{(j)} \mathcal{U}_{\mathbf{u},r_j}^{(j)} = (\mathcal{U}_{\mathbf{u},r_j}^{(j)})^T \mathcal{Y}_{\alpha_u}^{(j)} (\mathcal{X}_{\alpha_u}^{(j)})^\dagger \mathcal{U}_{\mathbf{u},r_j}^{(j)} \\ = (\mathcal{U}_{\mathbf{u},r_j}^{(j)})^T \mathcal{Y}_{\alpha_u}^{(j)} \mathcal{V}_{\mathbf{u},r_j}^{(j)} (\Sigma_{\mathbf{u},r_j}^{(j)})^{-1} (\mathcal{U}_{\mathbf{u},r_j}^{(j)})^T \mathcal{U}_{\mathbf{u},r_j}^{(j)} \\ = (\mathcal{U}_{\mathbf{u},r_j}^{(j)})^T \mathcal{Y}_{\alpha_u}^{(j)} \mathcal{V}_{\mathbf{u},r_j}^{(j)} (\Sigma_{\mathbf{u},r_j}^{(j)})^{-1} \quad (10) $$

3. Perform the eigendecomposition of $\tilde{\mathcal{K}}_{\mathbf{u}}^{(j)} \mathcal{W}_{\mathbf{u}}^{(j)} = \mathcal{W}_{\mathbf{u}}^{(j)} \Lambda_{\mathbf{u}}^{(j)}$ with $\Lambda_{\mathbf{u}}^{(j)}$ being the eigenvalues matrix and $\mathcal{W}_{\mathbf{u}}^{(j)}$ being the eigenvectors matrix

4. Obtain the DMD modes $\Phi_{\mathbf{u}}^{(j)} = \mathcal{Y}_{\alpha_u}^{(j)} \mathcal{V}_{\mathbf{u},r_j}^{(j)} (\Sigma_{\mathbf{u},r_j}^{(j)})^{-1} \mathcal{W}_{\mathbf{u}}^{(j)}$ corresponding to the particular DMD eigenvalues in $\Lambda_{\mathbf{u}}^{(j)}$

where $s \ge 1$ is tunable. Moreover, (12) can be rewritten into the following compact form

$$ \hat{\alpha}_{\mathbf{u}}(t^{(n_i+1)}, \theta^{(j)}) = \hat{\mathcal{K}}_{\mathbf{u}}^{(j)} \hat{\alpha}_{\mathbf{u}}(t^{(n_i)}, \theta^{(j)}), \quad i = 1, 2, \dots, \mathcal{N}_t - s, j = 1, 2, \dots, \mathcal{N}_p, \mathbf{u} \in \{\mathbf{E}, \mathbf{H}\}, \quad (13) $$

where the modified vectors $\hat{\alpha}_{\mathbf{u}}(t^{(n_i)}, \theta^{(j)})$ and the modified Koopman matrix $\hat{\mathcal{K}}_{\mathbf{u}}^{(j)}$ are respectively defined as

$$ \hat{\mathcal{K}}_{\mathbf{u}}^{(j)} = \begin{pmatrix} \mathbf{0} & \mathbf{I} & \mathbf{0} & \cdots & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{I} & \cdots & \mathbf{0} & \mathbf{0} \\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \cdots & \mathbf{I} & \mathbf{0} \\ \mathcal{K}_{\mathbf{u},1}^{(j)} & \mathcal{K}_{\mathbf{u},2}^{(j)} & \mathcal{K}_{\mathbf{u},3}^{(j)} & \cdots & \mathcal{K}_{\mathbf{u},s-1}^{(j)} & \mathcal{K}_{\mathbf{u},s}^{(j)} \end{pmatrix} \in \mathcal{R}^{(sL_{\mathbf{u}}) \times (sL_{\mathbf{u}})}, \quad (14) $$

and

$$ \hat{\alpha}_{\mathbf{u}}(t^{(n_i)}, \theta^{(j)}) = [\alpha_{\mathbf{u}}(t^{(n_i)}, \theta^{(j)}), \alpha_{\mathbf{u}}(t^{(n_i+1)}, \theta^{(j)}), \dots, \alpha_{\mathbf{u}}(t^{(n_i+s-1)}, \theta^{(j)})]^T \in \mathcal{R}^{(sL_{\mathbf{u}})}, \quad (15) $$

with $\mathbf{0}$ being the $L_{\mathbf{u}} \times L_{\mathbf{u}}$ null matrix and $\mathbf{I}$ being the $L_{\mathbf{u}} \times L_{\mathbf{u}}$ unit matrix, leading to (7) and (8) are respectively modified as

$$ \hat{\mathcal{X}}_{\alpha_u}^{(j)} = [\hat{\alpha}_{\mathbf{u}}(t^{(n_1)}, \theta^{(j)}), \hat{\alpha}_{\mathbf{u}}(t^{(n_2)}, \theta^{(j)}), \dots, \hat{\alpha}_{\mathbf{u}}(t^{(n_{N_t-s})}, \theta^{(j)})] \in \mathcal{R}^{(sL_{\mathbf{u}}) \times (N_t-s)}, \quad (16) $$

and

$$ \hat{\mathcal{Y}}_{\alpha_u}^{(j)} = [\hat{\alpha}_{\mathbf{u}}(t^{(n_2)}, \theta^{(j)}), \hat{\alpha}_{\mathbf{u}}(t^{(n_3)}, \theta^{(j)}), \dots, \hat{\alpha}_{\mathbf{u}}(t^{(n_{N_t-s+1})}, \theta^{(j)})] \in \mathcal{R}^{(sL_{\mathbf{u}}) \times (N_t-s)}. \quad (17) $$

So, we can perform the Algorithm 2 to find an approximation of the operator $\hat{\mathcal{K}}_{\mathbf{u}}^{(j)}$ in the HODMD method.

**Remark 2.** In Algorithm 2, the number $r_j$ is selected in terms of the rank of the matrix $\mathcal{X}_{\alpha_u}^{(j)}$, $R_j \le \min\{L_{\mathbf{u}}, \mathcal{N}_t - 1\}$, to be such that

$$ \mathcal{E}\mathcal{E}(r_j) = \sum_{i=r_j+1}^{R_j} \sigma_{j,i}^2 / \sum_{i=1}^{R_j} \sigma_{j,i}^2 \le \rho_{\text{SVD}}, \quad (18) $$
where $\mathcal{E}\mathcal{E}(r_j)$ is the relative root mean square (RMS) error, $\rho_{SVD}$ is the error threshold, and $\sigma_{j,i}$ is the singular value of $\mathcal{X}_{\alpha u}^{(j)}$ sorted in decreasing order. Particularly, the reduced-order coefficient vectors are of low dimension so that $\rho_{SVD}$ is chosen as $1 \times 10^{-10}$ in this study. Besides, a possibly smaller number $\mathcal{N}_{DMD}$ ($\le r_j$) of DMD modes, after the four step of Algorithm 2, can be selected imposing that

$$ \mathcal{AMP}_{(\mathcal{N}_{DMD}+1)}/\mathcal{AMP}_{(1)} < \rho_{DMD}, \qquad (19) $$

where $\mathcal{AMP}$ are the DMD amplitudes sorted in decreasing order, which can be computed by $\mathcal{AMP} = (\Phi_{\mathbf{u}}^{(j)})^{\dagger}\alpha_{\mathbf{u}}(t^{(n_1)}, \theta^{(j)}).$ For other alternative method of the computation of DMD mode amplitudes, one can see [24].

**Remark 3.** In Algorithm 2, for the HODMD method, the time complexities of steps 1 to 4 are $O(sL_u(\mathcal{N}_t - s)^2)$, $O(r_j(sL_u)(\mathcal{N}_t - s) + (\mathcal{N}_t - s)r_j^2 + r_j^2)$, $O(r_j^3)$, and $O((sL_u)(\mathcal{N}_t - s)r_j + (sL_u)r_j^2)$, respectively.

### 3.3. Radial basis function interpolation

Based on (11), one can predict the reduced-order coefficient vector in future time instants for the training parameter $\theta^{(j)} \in \mathcal{P}_h^{tr}$, and then obtain the RB solution via (4). However, it is inapplicable for parameterized problems. In order to resolve this problem, we propose the radial basis function (RBF) to interpolate the reduced-order coefficient vector for a new or untrained parameter $\theta^* \in \mathcal{P} \setminus \mathcal{P}_h^{tr}$.

Let $\mathcal{F}_{\alpha_{\mathbf{u}}^{\mathrm{DMD}}}(\theta)$ denotes the interpolation function representing $\alpha_{\mathbf{u}}^{\mathrm{DMD}}(t, \theta)$, which is a linear combination of the $\mathcal{N}_p$ RBFs $\varphi(\cdot)$ and has the following form

$$ \mathcal{F}_{\alpha_{\mathbf{u}}^{\mathrm{DMD}}}(\theta) = \sum_{k=1}^{\mathcal{N}_p} \varpi_{\mathbf{u},k} \varphi(\|\theta - \theta^{(k)}\|), \quad \mathbf{u} \in \{\mathbf{E}, \mathbf{H}\}, \qquad (20) $$

where $\varpi_{\mathbf{u},k}$ is the weighting vector associated with the parameter $\theta^{(k)}$, and $\|\cdot\|$ is the norm of a vector, usually chosen as Euclidean distance. Some well-known RBFs include Gaussian, linear spline, cubic spline, thin-plate spline, multi-quadratic, inverse multistory, and so on. In this work, we choose the Gaussian function, which has a form of $\varphi(d) = e^{-(d/\sigma)^2}$ with $d$ being the radius or distance and $\sigma$ being a shape parameter. In particular, choosing an optimal value of the shape parameter $\sigma$ is very critical in the RBF method. Throughout this work, the shape parameter $\sigma$ is selected by trials and tests. The weighting vector $\varpi_{\mathbf{u}}(j, i) = [\varpi_{\mathbf{u},1}(j, i), \varpi_{\mathbf{u},2}(j, i), \dots, \varpi_{\mathbf{u},\mathcal{N}_p}(j, i)]^T$ for the j-th component of the reduced-order coefficient vector $\alpha_{\mathbf{u}}^{\mathrm{DMD}}$ at the $n_i$-th time point is then determined by solving the following linear system

$$ \mathcal{A}_{\mathrm{RBF}} \varpi_{\mathbf{u}}(j, i) = \mathcal{B}_{\mathbf{u},\mathrm{RBF}}(j, i), \quad \mathbf{u} \in \{\mathbf{E}, \mathbf{H}\}, \qquad (21) $$

where the coefficient matrix $\mathcal{A}_{\mathrm{RBF}}$ and the right hand side $\mathcal{B}_{\mathbf{u},\mathrm{RBF}}(j, i)$ are respectively defined as

$$ \mathcal{A}_{\mathrm{RBF}} = \begin{pmatrix} \varphi(\|\theta^{(1)} - \theta^{(1)}\|) & \varphi(\|\theta^{(1)} - \theta^{(2)}\|) & \dots & \varphi(\|\theta^{(1)} - \theta^{(\mathcal{N}_p)}\|) \\ \varphi(\|\theta^{(2)} - \theta^{(1)}\|) & \varphi(\|\theta^{(2)} - \theta^{(2)}\|) & \dots & \varphi(\|\theta^{(2)} - \theta^{(\mathcal{N}_p)}\|) \\ \dots & \dots & \dots & \dots \\ \varphi(\|\theta^{(\mathcal{N}_p)} - \theta^{(1)}\|) & \varphi(\|\theta^{(\mathcal{N}_p)} - \theta^{(2)}\|) & \dots & \varphi(\|\theta^{(\mathcal{N}_p)} - \theta^{(\mathcal{N}_p)}\|) \end{pmatrix} \in \mathcal{R}^{\mathcal{N}_p \times \mathcal{N}_p}, \quad (22) $$

and

$$ \mathcal{B}_{\mathbf{u},\mathrm{RBF}}(j, i) = [\alpha_{\mathbf{u},j}^{\mathrm{DMD}}(t^{(n_i)}, \theta^{(1)}), \alpha_{\mathbf{u},j}^{\mathrm{DMD}}(t^{(n_i)}, \theta^{(2)}), \dots, \alpha_{\mathbf{u},j}^{\mathrm{DMD}}(t^{(n_i)}, \theta^{(\mathcal{N}_p)})]^T \in \mathcal{R}^{\mathcal{N}_p}. \quad (23) $$
For a new parameter $\theta^* \in \mathcal{P}$, the j-th component of the reduced-order coefficient vector $\alpha_{\mathbf{u},j}^{\text{DMD}}(t^{(n_i)}, \theta^*)$ at the $n_i$-th time point can be calculated as

$$ \alpha_{\mathbf{u},j}^{\text{DMD}}(t^{(n_i)}, \theta^*) \approx \alpha_{\mathbf{u},j}^{\text{RBF}}(t^{(n_i)}, \theta^*) = \sum_{k=1}^{N_p} \varpi_{\mathbf{u},k}(j, i) \varphi(\|\theta^* - \theta^{(k)}\|), \quad \mathbf{u} \in \{\mathbf{E}, \mathbf{H}\}. \quad (24) $$

Finally, for the new time $t^*$ and parameter $\theta^*$, the RB solution $\mathbf{u}_h^r(t^*, \theta^*)$ can be quickly recovered as

$$ \begin{aligned} \underline{\mathbf{u}}_h^r(t^*, \theta^*) &= \Psi_{\mathbf{u}} \alpha_{\mathbf{u}}(t^*, \theta^*) + \hat{\underline{\mathbf{u}}}_h \\ &\approx \Psi_{\mathbf{u}} \alpha_{\mathbf{u}}^{\text{DMD}}(t^*, \theta^*) + \hat{\underline{\mathbf{u}}}_h \\ &\approx \Psi_{\mathbf{u}} \alpha_{\mathbf{u}}^{\text{RBF}}(t^*, \theta^*) + \hat{\underline{\mathbf{u}}}_h, \quad \mathbf{u} \in \{\mathbf{E}, \mathbf{H}\}. \end{aligned} \quad (25) $$

### 3.4. Workflow of the proposed NIMOR method

Based on the discussions in the first three subsections, we propose the workflow for the construction of NIMOR method in Algorithm 3, where the offline and online stages of the POD, DMD, and RBF-based NIMOR are decoupled. Figure 2 shows the overview of the proposed method. During the offline stage, we collect suitable snapshots, obtain the RB bases, and build the DMD models for all training parameter values. During the online stage, in order to predict a RB solution at a new time and parameter value, we first obtain the reduced-order coefficient vector at the new time via the DMD models, then interpolate the new reduced-order coefficient vector in the parameter space, and finally recover the RB solution via (25).

**Algorithm 3 NIMOR for time-domain electromagnetic wave propagation problems**

**Offline stage**

Generate the HF solutions $\underline{\mathbf{u}}_h(t^{(n_i)}, \theta^{(j)}), (t^{(n_i)}, \theta^{(j)}) \in \mathcal{T}_h^{tr} \times \mathcal{P}_h^{tr}, \mathbf{u} \in \{\mathbf{E}, \mathbf{H}\}$

Compute the POD basis functions $\Psi_{\mathbf{u}}$ via Algorithm 1

Obtain the training data $\alpha_{\mathbf{u}}(t^{(n_i)}, \theta^{(j)})$ via (5)

Compute the DMD modes $\Phi_{\mathbf{u}}^{(j)} (j = 1, 2, \dots, N_p)$ and the eigenvalue matrices $\Lambda_{\mathbf{u}}^{(j)}$ via Algorithm 2 with (13) (HODMD)

**Online stage**

Compute the DMD predictions $\alpha_{\mathbf{u}}^{\text{DMD}}(t^*, \theta^{(j)})$ for a new time $t^*$ via (11), $\theta^{(j)} \in \mathcal{P}_h^{tr}$

Calculate the weighting vectors $\varpi_{\mathbf{u}}(j, i)$ via (21) to (23)

Calculate the RBF interpolation $\alpha_{\mathbf{u}}^{\text{RBF}}(t^*, \theta^*)$ for a new parameter $\theta^* \in \mathcal{P}$

Evaluate the RB solution $\underline{\mathbf{u}}_h^r(t^*, \theta^*)$ via (25)

**Remark 4.** DMD method is used to establish the surrogate models of reduced-order coefficient vectors for all training parameter values, and RBF method is adopted to interpolate reduced-order coefficient vector for a new parameter value, which is more effective than direct predicting and interpolating snapshots.

## 4. Numerical experiments

In this section, some numerical experiments including the scattering of a plane wave by a dielectric disk, by a multi-layer disk and by a 3-D dielectric sphere are presented to validate the effectiveness and the accuracy of the POD, DMD and RBF-based NIMOR method. In order to
<figure><img src="image_4.png" /><figcaption>Figure 2: Overview of the proposed NIMOR method.</figcaption></figure>

evaluate the accuracy of numerical experiments, the relative error between the NIMOR solution and the HF solution for $\mathbf{u} \in \{\mathbf{E}, \mathbf{H}\}$ is defined

$$e_{\mathbf{u},\text{NIMOR}}(t, \theta) = \frac{\| \underline{\mathbf{u}}_h(t, \theta) - \underline{\mathbf{u}}_h^r(t, \theta) \|_{\mathcal{R}^{\mathcal{N}_h}}}{\| \underline{\mathbf{u}}_h(t, \theta) \|_{\mathcal{R}^{\mathcal{N}_h}}} = \frac{\| \underline{\mathbf{u}}_h(t, \theta) - (\Psi_{\mathbf{u}} \alpha_{\mathbf{u}}^{\text{RBF}}(t, \theta) + \hat{\underline{\mathbf{u}}}_h) \|_{\mathcal{R}^{\mathcal{N}_h}}}{\| \underline{\mathbf{u}}_h(t, \theta) \|_{\mathcal{R}^{\mathcal{N}_h}}}, \quad (26)$$

which will be compared with the relative projection error generated by the nested POD method

$$e_{\mathbf{u},\text{Pro}}(t, \theta) = \frac{\| (\underline{\mathbf{u}}_h(t, \theta) - \hat{\underline{\mathbf{u}}}_h) - \Psi_{\mathbf{u}} \Psi_{\mathbf{u}}^T (\underline{\mathbf{u}}_h(t, \theta) - \hat{\underline{\mathbf{u}}}_h) \|_{\mathcal{R}^{\mathcal{N}_h}}}{\| \underline{\mathbf{u}}_h(t, \theta) \|_{\mathcal{R}^{\mathcal{N}_h}}}, \quad (27)$$

where $\| \cdot \|_{\mathcal{R}^{\mathcal{N}_h}}$ denotes Euclidean norm. Moreover, the average relative errors on a testing time and parameter sampling $\mathcal{T}_h^{te} \times \mathcal{P}_h^{te}$ of size $N_{te}$ are calculated to measure the convergence of the NIMOR method, which are defined as

$$\bar{e}_{\mathbf{u},\text{NIMOR}} = \frac{\sum_{(t,\theta) \in \mathcal{T}_h^{te} \times \mathcal{P}_h^{te}} e_{\mathbf{u},\text{NIMOR}}(t, \theta)}{N_{te}}, \quad (28)$$

and

$$\bar{e}_{\mathbf{u},\text{Pro}} = \frac{\sum_{(t,\theta) \in \mathcal{T}_h^{te} \times \mathcal{P}_h^{te}} e_{\mathbf{u},\text{Pro}}(t, \theta)}{N_{te}}, \quad (29)$$

where the testing parameter set $\mathcal{P}_h^{te}$ is generated via randomized latin-hypercube-sampling (LHS) method [37], and the testing time set $\mathcal{T}_h^{te}$ is uniformly chosen within a period of the physical simulation.

The HF solutions to the scattering problem of a dielectric disk and a multi-layer disk are calculated by using a discontinuous Galerkin time-domain (DGTD) solver formulated on an unstructured simplicial mesh [38], which combines a centered flux scheme with a second order leap-frog (LF2) time scheme. The HF solutions to the scattering problem of a 3-D dielectric sphere are obtained by using the Mie series method [39], and the corresponding computer implementation has been developed by Zhu [40]. The DGTD, Mie series and NIMOR methods are implemented in MATLAB and all simulations are run on a workstation equipped with an Intel Core i7-10700F CPU running at 2.90 GHz, and with 16 GB of RAM memory. All SVDs are implemented via the MATLAB function svd.
### 4.1. Scattering of a plane wave by a dielectric disk

We first consider the electromagnetic scattering of a plane wave by a dielectric disk, which is described by the following normalized form of 2-D time-domain transverse magnetic (TM) waves

$$ egin{cases} \mu_r \frac{\partial \mathbf{H}(\mathbf{x}, t)}{\partial t} + \nabla \times E(\mathbf{x}, t) = 0, & \forall (\mathbf{x}, t) \in \Omega \times \mathcal{T}, \\ \varepsilon_r \frac{\partial E(\mathbf{x}, t)}{\partial t} - \nabla \times \mathbf{H}(\mathbf{x}, t) = 0, & \forall (\mathbf{x}, t) \in \Omega \times \mathcal{T}, \end{cases} \quad (30) $$

where $E = E_z(\mathbf{x}, t)$ and $\mathbf{H}(\mathbf{x}, t) = [H_x(\mathbf{x}, t), H_y(\mathbf{x}, t)]^T$. The differential operators in this 2-D setting are

$$ \nabla \times E(\mathbf{x}, t) = \left[ \frac{\partial E_z(\mathbf{x}, t)}{\partial y}, -\frac{\partial E_z(\mathbf{x}, t)}{\partial x} \right]^T, \quad \text{and} \quad \nabla \times \mathbf{H}(\mathbf{x}, t) = \frac{\partial H_y(\mathbf{x}, t)}{\partial x} - \frac{\partial H_x(\mathbf{x}, t)}{\partial y}. \quad (31) $$

The incident plane waves take the form $E_z^{\text{inc}} = \cos(\omega t - kx)$ and $H_y^{\text{inc}} = -\cos(\omega t - kx)$ with $\omega = 2\pi f$ being the angular frequency of the incident wave frequency $f = 300$ MHz, and $k$ being the wave number in vacuum. The computational domain is the square $\Omega = [-2.6 \text{ m}, 2.6 \text{ m}] \times [-2.6 \text{ m}, 2.6 \text{ m}]$. The radius of the disk is set to be 0.6 m, and the relative permeability is set to be 1 (non-magnetic material). The range of interest relative permittivity of the disk is $[1, 5]$, so, $\mathcal{P} = \{\theta : \theta = \varepsilon_r \in [1, 5]\} \subset \mathcal{R}$. The medium exterior to the disk is assumed to be vacuum.

During the offline stage, the DGTD method is used to obtain the HF solutions in the training parameter set $\mathcal{P}_h^{tr}$ with uniformly sampling, in which the total simulation time is set to be 50 periods (which corresponds to 50 m in normalized unit) and the time step $\Delta t$ is set to be $3.678 \times 10^{-3}$ m. For each selected training parameter, the training snapshot vectors are equally chosen from 49 m to 49.7 m (whose are the first 70% data in the last period), leading to $\mathcal{T}_h^{tr} = \{49.0024 \text{ m}, 49.006 \text{ m}, \dots, 49.6938 \text{ m}, 49.6975 \text{ m}\}$. A testing parameter set $\mathcal{P}_h^{te}$ with size 40 generated by the LHS method and a testing time set $\mathcal{T}_h^{te}$ with size 263 including all time points in the last period (i.e., $\mathcal{T}_h^{te} = \mathcal{T}_h^{tr} \cup \{49.7012 \text{ m}, \dots, 49.966 \text{ m}\}$) are selected to evaluate the proposed method. The details for the training and testing datasets are listed in Table 1. The decay of the singular values in the nested POD method is shown in Figure 3, in

<table><caption>Table 1: Settings for the training and testing datasets.</caption><thead><tr><th>Case</th><th>Method</th><th>Data set</th><th>Training set</th><th>Testing set</th></tr></thead><tbody><tr><td rowspan="3">Disk (2-D)</td><td rowspan="3">DGTD</td><td>Parameter sampling</td><td>65, uniform</td><td>40, random (LHS)</td></tr><tr><td>Time sampling</td><td>190, uniform</td><td>263, uniform</td></tr><tr><td>Size</td><td>12350</td><td>10520</td></tr><tr><td rowspan="3">Multi-layer disk (2-D)</td><td rowspan="3">DGTD</td><td>Parameter sampling</td><td>137, uneven (Smolyak, 4-D, L = 3)</td><td>81, random (LHS)</td></tr><tr><td>Time sampling</td><td>183, uniform</td><td>253, uniform</td></tr><tr><td>Size</td><td>25071</td><td>20493</td></tr><tr><td rowspan="3">Sphere (3-D)</td><td rowspan="3">Mie series</td><td>Parameter sampling</td><td>81, uniform</td><td>40, random (LHS)</td></tr><tr><td>Time sampling</td><td>70, uniform</td><td>100, uniform</td></tr><tr><td>Size</td><td>5670</td><td>4000</td></tr></tbody></table>

which one can see that the singular values decrease rapidly indicating that a small dimension RB basis can capture the information of the original data. Particularly, the effects of $(\rho_t, \rho_\theta)$ on the average projection and NIMOR errors have been investigated to provide guidance on the choice of truncation tolerances in our previous work [9]. So we directly select the truncation tolerances and no longer test the numerical convergences in this study. With the truncation tolerances $\rho_t = 1 \times 10^{-1}$ and $\rho_\theta = 1 \times 10^{-5}$ in the nested POD method, a set of $L_{H_x} = 59$, $L_{H_y} = 21$, and $L_{E_z} = 22$ POD basis functions are extracted. The numerical convergences of
<figure><img src="image_5.png" /><figcaption>Figure 3: Scattering of a plane wave by a dielectric disk: the decay of the singular values in the nested POD method. The first to third columns are the singular values of the 10th, 30th and 50th training parameter in the first step of the nested POD method. The fourth column is the singular values in the second step of the nested POD method.</figcaption></figure>

the NIMOR method with respect to the number $s$ of time lagged vectors in (12) and the error threshold $\rho_{\text{DMD}}$ in (19) are shown in Figure 4, where the shape parameter $\sigma$ is set to be 0.073 in (20). Figure 5 shows the size of the DMD modes with $s = 20$, $\rho_{\text{DMD}} = 1 \times 10^{-4}$, and $s = 50$,

<figure><img src="image_6.png" /><figcaption>Figure 4: Scattering of a plane wave by a dielectric disk: the numerical convergences of (a) $\bar{e}_{\text{E,NIMOR}}$, and (b) $\bar{e}_{\text{H,NIMOR}}$ on the testing set $\mathcal{T}_h^{\text{te}} \times \mathcal{P}_h^{\text{te}}$ with vary $s$ and $\rho_{\text{DMD}}$, where $\bar{e}_{\text{u,Pro}}$ is the average projection error of the HF solutions for $\mathbf{E}$ or $\mathbf{H}$ onto the corresponding RB subspace.</figcaption></figure>

$\rho_{\text{DMD}} = 1 \times 10^{-2}$ based on (19). As we can see, the dimension of the NIMOR method is further reduced. Combining Figures 4 and 5, we choose the number $s = 50$ and the error threshold $\rho_{\text{DMD}} = 1 \times 10^{-2}$. The corresponding average relative errors are listed in Table 2. Some of
<figure><img src="image_7.png" /><figcaption>Figure 5: Scattering of a plane wave by a dielectric disk: the size of the DMD modes with $s = 20$ and $\rho_{DMD} = 1 \times 10^{-4}$, $s = 50$ and $\rho_{DMD} = 1 \times 10^{-2}$, and the dimension of RB basis for $H_x$, $H_y$ and $E_z$.</figcaption></figure>

<table><caption>Table 2: Comparison of the average relative projection and NIMOR errors on the testing set $T_h^{te} \times P_h^{te}$.</caption><thead><tr><th>Case</th><th>ε<sub>E,Pro</sub></th><th>ε<sub>E,NIMOR</sub></th><th>ε<sub>H,Pro</sub></th><th>ε<sub>H,NIMOR</sub></th></tr></thead><tbody><tr><td>Disk</td><td>1.178 × 10<sup>-2</sup></td><td>1.348 × 10<sup>-2</sup></td><td>1.072 × 10<sup>-2</sup></td><td>1.388 × 10<sup>-2</sup></td></tr><tr><td>Multi-layer disk</td><td>7.622 × 10<sup>-3</sup></td><td>7.690 × 10<sup>-3</sup></td><td>6.727 × 10<sup>-3</sup></td><td>6.775 × 10<sup>-3</sup></td></tr><tr><td>Sphere</td><td>2.451 × 10<sup>-3</sup></td><td>5.908 × 10<sup>-3</sup></td><td>2.479 × 10<sup>-3</sup></td><td>5.269 × 10<sup>-3</sup></td></tr></tbody></table>
the exact and approximate (HODMD method) reduced-order coefficient vectors in $\mathcal{T}_h^{te}$ for all training parameters are displayed Figure 6. It is clear that the accuracy of the HODMD method is very well within the **testing** time region (where the set [49 m, 49.7 m] is the **training** time region, and the set [49.7 m, 50 m] is the prediction time region) for all training parameters. This implies the potential of the HODMD method in the prediction of the RB solutions outside of the training time region.

During the online stage, the reduced-order solutions are visualized for some testing parameters including $\theta^{(1)} = 1.215$, $\theta^{(2)} = 2.215$, $\theta^{(3)} = 3.215$, and $\theta^{(4)} = 4.215$. Figure 7 shows the time evolution fields $H_y$ and $E_z$ (HF and RB solutions) in $\mathcal{T}_h^{te}$ at a given space point for above four testing parameters. It can be **found** that the RB solutions and the HF solutions match each other very well. Furthermore, the corresponding time evolution of the relative projection and NIMOR errors for the four testing parameters are displayed in Figure 8. In order to have a look on visual effects of electromagnetic fields, **over the Fourier domain during the last oscillation period of the incident wave, we display** in Figure 9 the total field of $|H_x|$ calculated by the NIMOR method for the parameter $\theta^{(3)} = 3.215$, which coincides with the HF solution. Finally, Table 3 shows the average computing time of the NIMOR and HF solutions. Although we spend a long time preparing the training data for the reduced-order modeling in the offline stage, it is worthwhile especially when solving the electromagnetic field value under multiple parameters, because the online computing time of the NIMOR method is greatly shortened.

<table><caption>Table 3: Comparison of the average computing time of the RB and HF solutions. The unit of time cost is second.</caption><thead><tr><th>Case</th><th>Offlin stage<br>(HF solutions, Nested POD, DMD training)</th><th>Online stage<br>(DMD prediction, RBF, one run for new paramter)</th></tr></thead><tbody><tr><td>Disk</td><td>2.765 × 10<sup>+4</sup>, 2.015 × 10<sup>+1</sup>, 2.920 × 10<sup>+1</sup></td><td>6.690 × 10<sup>-1</sup>, 5.089 × 10<sup>-1</sup>, 4.218 × 10<sup>-2</sup> ≪ 1</td></tr><tr><td>Multi-layer disk</td><td>6.184 × 10<sup>+4</sup>, 1.071 × 10<sup>+2</sup>, 5.224 × 10<sup>+1</sup></td><td>7.040 × 10<sup>-1</sup>, 5.313 × 10<sup>-1</sup>, 5.680 × 10<sup>-2</sup> ≪ 1</td></tr><tr><td>Sphere</td><td>2.541 × 10<sup>+5</sup>, 2.194 × 10<sup>+2</sup>, 6.900 × 10<sup>+1</sup></td><td>7.800 × 10<sup>-1</sup>, 8.066 × 10<sup>-1</sup>, 5.098 × 10<sup>-2</sup> ≪ 1</td></tr></tbody></table>

## 4.2. Scattering of a plane wave by a multi-layer disk

Next, we consider a multivariate parameter problem, which is induced by a multi-layer disk illuminated by an incident plane wave, as shown in Figure 10. The computational domain is artificially truncated by a square with 6.4 m side length. The radius $r_i$ of each layer of the multi-layer disk and the corresponding relative permittivity range are summarized in Table 4. The medium exterior to the multi-layer disk is also assumed to be vacuum. So, the material

<table><caption>Table 4: Distribution and range of material parameters.</caption><thead><tr><th></th><th>Layer <i>i</i></th><th><i>P</i><sup>(<i>i</i>)</sup></th><th><i>μ</i><sub><i>r</i>,<i>i</i></sub></th><th><i>r</i><sub><i>i</i></sub></th></tr></thead><tbody><tr><td rowspan="4">Multi-layer disk</td><td>1</td><td>ε<sub><i>r</i>,1</sub> ∈ [5.0, 5.6]</td><td>1</td><td>0.15</td></tr><tr><td>2</td><td>ε<sub><i>r</i>,2</sub> ∈ [3.25, 3.75]</td><td>1</td><td>0.3</td></tr><tr><td>3</td><td>ε<sub><i>r</i>,3</sub> ∈ [2.0, 2.5]</td><td>1</td><td>0.45</td></tr><tr><td>4</td><td>ε<sub><i>r</i>,4</sub> ∈ [1.25, 1.75]</td><td>1</td><td>0.6</td></tr></tbody></table>

parameters can be represented as a 4-D vector $\theta = (\varepsilon_{r,1}, \varepsilon_{r,2}, \varepsilon_{r,3}, \varepsilon_{r,4}) \in \mathcal{P} = \mathcal{P}^{(1)} \times \mathcal{P}^{(2)} \times \mathcal{P}^{(3)} \times \mathcal{P}^{(4)}$. Other physical information is the same as the case in section 4.1.

During the offline stage, the HF solutions, in the training parameter set $\mathcal{T}_h^{tr}$ chosen by a Smolyak sparse grid method with approximation level $L = 3$, are generated by using the DGTD solver. The Smolyak sparse grid method is often used for the integration or interpolation of
<figure><img src="image_8.png" /><figcaption>(a)</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(c)</figcaption></figure>

Figure 6: Scattering of a plane wave by a dielectric disk: some of the exact and approximate (HODMD method) reduced-order coefficient vectors $\alpha_{u,l}(t, \theta)$ ($l = 1, 2, 4, 6$) of (a) $H_x$, (b) $H_y$ and (c) $E_z$ in $\mathcal{T}_h^{te}$ for all training parameters. The dashed line indicates the end of the training time region, which is $t = 49.7$ m.
<figure><img src="image_11.png" /><figcaption>Figure 7: Scattering of a plane wave by a dielectric disk: comparison of the time evolution of the fields $H_y$ (top) and $E_z$ (bottom) at a given space point for the four testing parameters in the testing time region $T_h^{te}$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 8: Scattering of a plane wave by a dielectric disk: comparison of the time evolution of the relative errors of **E** (left) and **H** (right) for the four testing parameters in the testing time region $T_h^{te}$.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 9: Comparison of the total fields of $|H_x|$ for disk (top), multi-layer disk (middle), and sphere (bottom) for the testing parameter $\theta^{(3)}$, including the HF solutions (left), the RB solutions (middle), and the absolute errors (right).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 10: Geometry of the multi-layer disk.</figcaption></figure>
multivariate functions, and can be used to deal with the so-called *curse of dimensionality* problem in high dimensional space, which is briefly described in the Appendix. The total simulation time is also set to be 50 m, and the time step $\Delta t$ is set to be $3.800 \times 10^{-3}$ m. Similarly, one can equally choose the training snapshots from 49 m to 49.7 m for each training parameter, and the corresponding training time set is $\mathcal{T}_{\bar{h}}^{\text{tr}} = \{49.000 \text{ m}, 49.004 \text{ m}, \dots, 49.6945 \text{ m}, 49.6984 \text{ m}\}$. A testing parameter set $\mathcal{P}_{\bar{h}}^{\text{te}}$ with size 81 generated by the LHS method and a testing time set $\mathcal{T}_{\bar{h}}^{\text{te}}$ with size 253 including all time points in the last period are chosen to evaluate the NIMOR method. Table 1 summarizes the training and testing sets. With the truncation errors $\rho_t = 1 \times 10^{-1}$ and $\rho_{\theta} = 1 \times 10^{-5}$ for the nested POD method, $L_{H_x} = 16$ POD basis functions are extracted for $H_x$, $L_{H_y} = 15$ for $H_y$, and $L_{E_z} = 15$ for $E_z$. As with the previous problem, Figure 11 shows the numerical histories of the NIMOR method with respect to the number $s$ and the error threshold $\rho_{\text{DMD}}$, in which we set the shape parameter $\sigma = 1$. It can be seen that the errors $\bar{e}_{u,\text{NIMOR}}$ are almost unchanged for varying $s$ and $\rho_{\text{DMD}}$ for this case. The main reason is that the first few terms of amplitude play a key role. Figure 12 (a) shows the DMD amplitudes for the first training parameter when $s = 5$. In addition, one can see that the NIMOR errors

<figure><img src="image_15.png" /><figcaption>Figure 11: Scattering of a plane wave by a multi-layer disk: the numerical histories of (a) $\bar{e}_{E,\text{NIMOR}}$, and (b) $\bar{e}_{H,\text{NIMOR}}$ on the testing set $\mathcal{T}_{\bar{h}}^{\text{te}} \times \mathcal{P}_{\bar{h}}^{\text{te}}$ with vary $s$ and $\rho_{\text{DMD}}$, where $\bar{e}_{u,\text{Pro}}$ is the average projection error of the HF solutions for $E$ or $H$ onto the corresponding RB subspace.</figcaption></figure>

very close to the projection errors $\bar{e}_{u,\text{Pro}}$. In this test, $s$ and $\rho_{\text{DMD}}$ are respectively set to be 5 and $1 \times 10^{-1}$, and the corresponding sizes of the DMD modes for all training parameters, which are shown in Figure 12 (b), are all reduced to $\mathcal{N}_{\text{DMD}} = 2$. The corresponding average relative errors are presented in Table 2, and some of the approximation and exact reduced-order coefficient vectors, in the testing time set $\mathcal{T}_{\bar{h}}^{\text{te}}$ (where the set [49 m, 49.7 m] is the **training** time region, and the set [49.7 m, 50 m] is the prediction time region) for all some parameters, are shown in Figure 13. We can draw conclusions similar to Figure 6 from Figure 13.

During the online stage, the RB solutions are calculated by using the NIMOR method for four testing parameters $\theta^{(1)} = \{(5.125, 3.375, 2.125, 1.375)\}$, $\theta^{(2)} = \{(5.425, 3.625, 2.425, 1.625)\}$, $\theta^{(3)} = \{(5.125, 3.625, 2.125, 1.625)\}$, and $\theta^{(4)} = \{(5.425, 3.375, 2.425, 1.375)\}$. An evidence of the validity of the NIMOR method is shown in Figure 14, which reports the comparison between the time evolution of RB solutions at a given space point and their HF solutions. The time evolution of the relative projection and NIMOR errors are displayed in Figure 15. The total field of $|H_x|$ calculated by the NIMOR method for the parameter $\theta^{(3)} = \{(5.125, 3.625, 2.125, 1.625)\}$
<figure><img src="image_16.png" /><figcaption>Figure 12: Scattering of a plane wave by a multi-layer disk: (a) the DMD amplitudes for the first training parameter when $s = 5$, and (b) the size of the DMD modes for all training parameters when $s = 5$ and $\rho_{DMD} = 1 \times 10^{-1}$.</figcaption></figure>

in the Fourier domain during the last period of wave oscillation is displayed in Figure 9. The time performance of the proposed NIMOR is presented in Table 3. Based on the numerical performance of the NIMOR method proposed here, this method is effective for multivariate parameter problem.

### 4.3. Scattering of a plane wave by a dielectric sphere

Finally, we consider the electromagnetic scattering of a plane wave by a 3-D dielectric sphere, where the plane wave is polarized in the $x$ direction propagating in the $+z$ direction. The incident field in frequency domain is expressed as

$$ \mathbf{E}^{\text{inc}} = \mathbf{A}_0 e^{-ikz} \mathbf{n}_x, \text{ and } \mathbf{H}^{\text{inc}} = \frac{E_x^{\text{inc}}}{\eta} e^{-ikz} \mathbf{n}_y, \qquad (32) $$

where the magnitude $\mathbf{A}_0$ is 1 V/m, the wave number $k = \omega\sqrt{\mu_r\varepsilon_r}$ with $\omega = 2\pi f$ being the angular frequency of the incident wave frequency $f = 300$ MHz, and $i$ is the imaginary unit; $\mathbf{n}_x = (1, 0, 0)^T$, and $\mathbf{n}_y = (0, 1, 0)^T$. In particular, one can obtain the following time evolution solutions once obtaining the frequency-domain solutions

$$ \underline{\mathbf{u}}(t, \theta) = \mathcal{R}(\underline{\mathbf{u}}(\theta)e^{i\omega t}), \quad \theta \in \mathcal{P}, \mathbf{u} \in \{\mathbf{E}, \mathbf{H}\}, \qquad (33) $$

where $\mathcal{R}(\cdot)$ denotes the real part of a vector. The computational domain is the square $\Omega^2 = [-0.2 \text{ m}, 2 \text{ m}]^3$. The radius of the sphere is set to be 1 m. The range of corresponding relative permittivity is $[2, 6]$, resulting in $\mathcal{P} = \{\theta : \theta = \varepsilon_r \in [2, 6]\} \subset \mathcal{R}$. The medium exterior to the sphere is assumed to be vacuum.

During the offline stage, the HF solutions, at the 81 uniformly distributed parameter points, are obtained by using the Mie series method with $\Delta t = 1.0101 \times 10^{-2}$ m in the first period, which are implemented by Zhu [40] for the frequency-domain solutions. In this test, the number of spatial grid points in $x$, $y$, and $z$ directions is all 45, i.e., $\mathcal{N}_h = 45^3 = 91125$. The training time set is evenly chosen as $\mathcal{T}_h^{tr} = \{0 \text{ m}, 1.0101 \times 10^{-2} \text{ m}, \cdots, 6.9697 \times 10^{-1} \text{ m}\}$. The testing parameter set $\mathcal{P}_h^{te}$ with size 40 is chosen by LHS method, and the testing time set $\mathcal{T}_h^{te}$ is set to be $\{0 \text{ m}, 1.0101 \times 10^{-2} \text{ m}, \cdots, 1 \text{ m}\}$. With the truncation errors $\rho_t = 1 \times 10^{-1}$ and $\rho_\theta = 1 \times 10^{-5}$, a
<figure><img src="image_17.png" /><figcaption>Figure 13: Scattering of a plane wave by a multi-layer disk: some of the exact (-) and approximate (o) (HODMD method) reduced-order coefficient vectors $\alpha_{u,l}(t, \theta)$ ($l = 1, 5, 10$) of (a) $H_x$, (b) $H_y$ and (c) $E_z$ in $T_h^{le}$ for some training parameters. Different colors indicate different parameters. The dashed line indicates the end of the training time region, which is $t = 49.7$ m.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Figure 14: Scattering of a plane wave by a multi-layer disk: comparison of the time evolution of the fields $H_x$ (top), $H_y$ (middle) and $E_z$ (bottom) at a given point for the four testing parameters in the testing time region $\mathcal{T}_h^{te}$.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Figure 15: Scattering of a plane wave by a multi-layer disk: comparison of the time evolution of the relative errors of **E** (left) and **H** (right) for the four testing parameters in the testing time region $\mathcal{T}_h^{te}$.</figcaption></figure>
set of $L_{H_x} = 38$, $L_{H_y} = 37$, $L_{H_z} = 34$, $L_{E_x} = 36$, $L_{E_y} = 32$, and $L_{E_z} = 31$ POD basis functions are extracted. Furthermore, $s$, $\rho_{DMD}$, and $\sigma$ are respectively set to be 5, $1 \times 10^{-1}$, and 0.039, resulting in $N_{DMD} = 2$. Table 2 shows the corresponding average relative errors, and Figure 16 and 17 present some of the approximation and exact reduced-order coefficient vectors in the testing set $\mathcal{T}_h^{te} \times \mathcal{P}_h^{te}$ for all training parameters, where the set $[0 \text{ m}, 0.7 \text{ m}]$ is the training time region, and the set $[0.7 \text{ m}, 1 \text{ m}]$ is the prediction time region.

During the online stage, in order to show the performance of NIMOR method for 3-D problem, we compare the RB solutions with the HF solutions generated by the Mie series method for four testing parameters $\theta^{(1)} = 2.145$, $\theta^{(2)} = 3.145$, $\theta^{(3)} = 4.145$, and $\theta^{(4)} = 5.145$. Figure 18 reports the comparison between the time evolution of the RB solutions at a given space point for the four testing parameters and the corresponding HF solutions. Figure 19 shows the time evolution of the corresponding relative projection and NIMOR errors. Figure 9 displays the total field of $|H_x|$ calculated by the NIMOR method for the parameter $\theta^{(3)} = 4.145$ in the Fourier domain during the first period of wave oscillation. Table 3 presents the average computing time for the offline and online stages.

## 5. Conclusion

In this paper, we proposed and studied an 'equation-free' NIMOR method for time-domain electromagnetic wave propagation. From the HF database on a given training parameter set with uniform or uneven sampling, a set of RB functions were extracted by using the nested POD method, which performs a prior dimensionality reduction. The DMD method, which is defined by finding the best-fit linear model of the Koopman operator, was then used to predict the low-dimensional reduced-order coefficient vectors in future time instants on this training parameter set, where the dimension of NIMOR method is further reduced. The RBF interpolation was finally used to estimate the reduced-order coefficient vectors at a new specified parameter value from the DMD calculations. The main contribution of this work is that we exploited the POD, DMD and RBF techniques to obtain a quick approximation for new specified parameters in future time instants. The accuracy and effectiveness of the NIMOR method were illustrated numerically by considering three classical test cases: the scattering of a plane wave by a dielectric disk, by a multi-layer disk, and by a 3-D dielectric sphere. Future research direction include the usage of other interpolation methods, e.g., manifold interpolation [29], to obtain the reduced DMD operator, or the reduction of geometric parameters.

## Acknowledgments

The authors are very grateful for the constructive advices from anonymous reviewers and the help of the editors.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgements

This research was supported by the National Natural Science Foundation of China (Grant No. 12101511) and by the Fundamental Research Funds for the Central Universities (Grant No. JBK22YJ37).
<figure><img src="image_20.png" /></figure>

(a)

<figure><img src="image_21.png" /></figure>

(b)

<figure><img src="image_22.png" /></figure>

(c)

Figure 16: Scattering of a plane wave by a dielectric sphere: some of the exact and approximate (HODMD method) reduced-order coefficient vectors $\alpha_{u,l}(t, \theta)$ ($l = 1, 2, 4, 6$) of (a) $E_x$, (b) $E_y$ and (c) $E_z$ in $\mathcal{T}_h^{te}$ for all training parameters. The dashed line indicates the end of the training time region, which is $t = 0.7$ m.
<figure><img src="image_23.png" /><figcaption>(a)</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>(c)</figcaption></figure>

Figure 17: Scattering of a plane wave by a dielectric sphere: some of the exact and approximate (HODMD method) reduced-order coefficient vectors $\alpha_{u,l}(t, \theta)$ ($l = 1, 2, 4, 6$) of (a) $H_x$, (b) $H_y$ and (c) $H_z$ in $\mathcal{T}_h^{te}$ for all training parameters. The dashed line indicates the end of the training time region, which is $t = 0.7$ m.
<figure><img src="image_26.png" /><figcaption>Figure 18: Scattering of a plane wave by a dielectric sphere: comparison of the time evolution of the fields $H_x$ (top, left), $H_y$ (middle, left), $H_z$ (bottom, left), $E_x$ (top, right), $E_y$ (middle, right), $E_z$ (bottom, right) at a given space point for the four testing parameters in the testing time region $T_h^{te}$.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>Figure 19: Scattering of a plane wave by a dielectric sphere: comparison of the time evolution of the relative errors of **E** (left) and **H** (right) for the four testing parameters in the testing time region $T_h^{te}$.</figcaption></figure>
# Appendix: Smolyak sparse grid

The sampling process of the training parameter points over the parameter space $P \in \mathbb{R}^p$ is briefly introduced. Let $Q_{i_k}$ be an interpolation operator in the kth-direction of p-dimensional space, and the interpolation operator of Smolyak algorithm in whole p-dimensional space can be then expressed as

$$Q(L, p) = \sum_{L-p+1 \le |i| \le L} (-1)^{L-|i|} \binom{p-1}{L-|i|} (Q_{i_1} \otimes Q_{i_2} \otimes \cdots \otimes Q_{i_p}), \quad (34)$$

where $\otimes$ is the tensor product operator, $i_k$ ($k = 1, 2, \cdots, p$) denotes the approximation level of sparse grid in the kth-direction, $L$ is the entire approximation level of sparse grid, and $|i| = i_1 + i_2 + \cdots + i_p$. The set of Smolyak sparse sampling points in (34) can be then derived as

$$S(L, p) = \bigcup_{L-p+1 \le |i| \le L} (S_{i_1} \otimes S_{i_2} \otimes \cdots \otimes S_{i_p}), \quad (35)$$

with

$$S_{i_k} = \{\alpha^{(1)}, \alpha^{(2)}, \cdots, \alpha^{(M_{i_k})}\},$$

where $S_{i_k}$ is the vector of sampling points with approximation level $i_k$ in the kth-direction, and $M_{i_k}$ denotes the corresponding number of sampling points. Particularly, one can choose the Clenshaw-Curtis points as the sampling points on each dimension. Examples of Smolyak grids under approximation levels $L = 1, L = 2$, and $L = 3$ are illustrated in Figure 20 for 2-D case. For comparison, a tensor-product grid of $5^2$ points is also shown. For a more detailed introduction of this technique, one can see [41].

<figure><img src="image_28.png" /><figcaption>Figure 20: 2-D Smolyak sparse grids versus a full tensor product grid.</figcaption></figure>

# References

[1] P. Benner, S. Gugercin, K. Willcox, A survey of projection-based model reduction methods for parametric dynamical systems, SIAM Review 57 (4) (2015) 483–531.

[2] J. S. Hesthaven, G. Rozza, B. Stamm, Certified reduced basis methods for parametrized partial differential equations, Springer, 2016.

[3] B. Peherstorfer, K. Willcox, M. Gunzburger, Survey of multifidelity methods in uncertainty propagation, inference, and optimization, SIAM Review 60 (3) (2018) 550–591.

[4] J. S. Hesthaven, C. Pagliantini, G. Rozza, Reduced basis methods for time-dependent problems, Acta Numerica 31 (2022) 265–345.
[5] J. S. Hesthaven, S. Ubbiali, Non-intrusive reduced order modeling of nonlinear problems using neural networks, Journal of Computational Physics 363 (2018) 55–78.

[6] Q. Wang, J. S. Hesthaven, D. Ray, Non-intrusive reduced order modeling of unsteady flows using artificial neural networks with application to a combustion problem, Journal of Computational Physics 384 (2019) 289–307.

[7] G. Berkooz, P. Holmes, J. L. Lumley, The proper orthogonal decomposition in the analysis of turbulent flows, Annual Review of Fluid Mechanics 25 (1) (1993) 539–575.

[8] Y. Liang, H. P. Lee, S. Lim, W. Lin, K. Lee, C. Wu, Proper orthogonal decomposition and its applications - Part I: Theory, Journal of Sound and Vibration 252 (3) (2002) 527–544.

[9] K. Li, T.-Z. Huang, L. Li, S. Lanteri, Non-intrusive reduced-order modeling of param- eterized electromagnetic scattering problems using cubic spline interpolation, Journal of Scientific Computing 87 (2) (2021). doi:https://doi.org/10.1007/s10915-021-01467-2.

[10] F. Casenave, A. Ern, T. Lelievre, A nonintrusive reduced basis method applied to aeroa- coustic simulations, Advances in Computational Mathematics 41 (5) (2015) 961–986.

[11] W. Chen, Q. Wang, J. S. Hesthaven, C. Zhang, Physics-informed machine learning for reduced-order modeling of nonlinear problems, Journal of Computational Physics 446 (2021) 110666. doi:https://doi.org/10.1016/j.jcp.2021.110666.

[12] M. Raissi, P. Perdikaris, G. E. Karniadakis, Physics-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations, Journal of Computational Physics 378 (2019) 686–707.

[13] J. Fu, D. Xiao, R. Fu, C. Li, C. Zhu, R. Arcucci, I. M. Navon, Physics-data combined machine learning for parametric reduced-order modelling of nonlinear dynamical systems in small-data regimes, Computer Methods in Applied Mechanics and Engineering 404 (2023) 115771. doi:https://doi.org/10.1016/j.cma.2022.115771.

[14] M. Guo, J. S. Hesthaven, Data-driven reduced order modeling for time-dependent problems, Computer Methods in Applied Mechanics and Engineering 345 (2019) 75–99.

[15] J. Yu, C. Yan, Z. Jiang, W. Yuan, S. Chen, Adaptive non-intrusive reduced order modeling for compressible flows, Journal of Computational Physics 397 (2019) 108855. doi:https://doi.org/10.1016/j.jcp.2019.07.053.

[16] K. Li, T.-Z. Huang, L. Li, Y. Zhao, S. Lanteri, A non-intrusive model order reduction approach for parameterized time-domain Maxwell’s equations, Discrete and Continuous Dynamical Systems-Series B 28 (1) (2022) 449–473.

[17] D. Xiao, P. Yang, F. Fang, J. Xiang, C. C. Pain, I. M. Navon, M. Chen, A non-intrusive reduced-order model for compressible fluid and fractured solid coupling and its application to blasting, Journal of Computational Physics 330 (2017) 221–244.

[18] D. Xiao, Error estimation of the parametric non-intrusive reduced order model using ma- chine learning, Computer Methods in Applied Mechanics and Engineering 355 (2019) 513– 534.
[19] X.-F. He, L. Li, S. Lanteri, K. Li, Model order reduction for parameterized electromagnetic problems using matrix decomposition and deep neural networks, Journal of Computational and Applied Mathematics 431 (2023) 115271. doi:https://doi.org/10.1016/j.cam.2023.115271.

[20] P. J. Schmid, Dynamic mode decomposition of numerical and experimental data, Journal of Fluid Mechanics 656 (2010) 5–28.

[21] J. N. Kutz, S. L. Brunton, B. W. Brunton, J. L. Proctor, Dynamic mode decomposition: data-driven modeling of complex systems, SIAM, 2016.

[22] D. Duke, J. Soria, D. Honnery, An error analysis of the dynamic mode decomposition, Experiments in Fluids 52 (2012) 529–542.

[23] H. Lu, D. M. Tartakovsky, Prediction accuracy of dynamic mode decomposition, SIAM Journal on Scientific Computing 42 (3) (2020) A1639–A1662.

[24] S. Le Clainche, J. M. Vega, Higher order dynamic mode decomposition, SIAM Journal on Applied Dynamical Systems 16 (2) (2017) 882–925.

[25] Z. Drmac, I. Mezic, R. Mohr, Data driven modal decompositions: analysis and enhancements, SIAM Journal on Scientific Computing 40 (4) (2018) A2253–A2285.

[26] H. Lu, D. M. Tartakovsky, Lagrangian dynamic mode decomposition for construction of reduced-order models of advection-dominated phenomena, Journal of Computational Physics 407 (2020) 109229. doi:https://doi.org/10.1016/j.jcp.2020.109229.

[27] S. Le Clainche, J. M. Vega, Higher order dynamic mode decomposition to identify and extrapolate flow patterns, Physics of Fluids 29 (8) (2017) 084102. doi:https://doi.org/10.1063/1.4997206.

[28] Z. Gao, Y. Lin, X. Sun, X. Zeng, A reduced order method for nonlinear parameterized partial differential equations using dynamic mode decomposition coupled with k-nearest-neighbors regression, Journal of Computational Physics 452 (2022) 110907. doi:https://doi.org/10.1016/j.jcp.2021.110907.

[29] M. W. Hess, A. Quaini, G. Rozza, A data-driven surrogate modeling approach for time-dependent incompressible navier-stokes equations with dynamic mode decomposition and manifold interpolation, Advances in Computational Mathematics 49 (2) (2023) 22. doi:https://doi.org/10.1007/s10444-023-10016-4.

[30] Q. A. Huhn, M. E. Tano, J. C. Ragusa, Y. Choi, Parametric dynamic mode decomposition for reduced order modeling, Journal of Computational Physics 475 (2023) 111852. doi:https://doi.org/10.1016/j.jcp.2022.111852.

[31] Y. Lin, Z. Gao, Y. Chen, X. Sun, A dynamic mode decomposition based reduced-order model for parameterized time-dependent partial differential equations, Journal of Scientific Computing 95 (3) (2023) 70. doi:https://doi.org/10.1007/s10915-023-02200-x.

[32] T. Sayadi, P. J. Schmid, F. Richecoeur, D. Durox, Parametrized data-driven decomposition for bifurcation analysis, with application to thermo-acoustically unstable systems, Physics of Fluids 27 (3) (2015) 037102. doi:https://doi.org/10.1063/1.4913868.

[33] J. Duan, J. S. Hesthaven, Non-intrusive data-driven reduced-order modeling for time-dependent parametrized problems, arXiv preprint arXiv:2303.02986 (2023).
[34] B. O. Koopman, Hamiltonian systems and transformation in Hilbert space, Proceedings of the National Academy of Sciences 17 (5) (1931) 315-318.

[35] R. Penrose, On best approximate solutions of linear matrix equations, in: Mathematical Proceedings of the Cambridge Philosophical Society, Vol. 52, Cambridge University Press, 1956, pp. 17–19.

[36] F. Andreuzzi, N. Demo, G. Rozza, A dynamic mode decomposition extension for the forecasting of parametric dynamical systems, arXiv preprint arXiv:2110.09155 (2021).

[37] L. L. Sharon, Sampling: design and analysis, Toronto, Canada: Nelson Education, 2009.

[38] J. S. Hesthaven, T. Warburton, Nodal discontinuous Galerkin methods: algorithms, analysis, and applications, Springer Science & Business Media, 2007.

[39] W. C. Chew, Waves and fields in inhomogenous media, John Wiley & Sons, 1999.

[40] G. K. Zhu, Sphere scattering, https://www.mathworks.com/matlabcentral/fileexchange/31119-sphere-scattering, 2023.

[41] K. L. Judd, L. Maliar, S. Maliar, R. Valero, Smolyak method for solving dynamic economic models: Lagrange interpolation, anisotropic grid and adaptive domain, Journal of Economic Dynamics and Control 44 (2014) 92–123.