Efficient Computation of the Magnetic Polarizability Tensor Spectral Signature using POD

The eddy current model is a low frequency approximation of the Maxwell system that neglects the displacement currents, which is valid when the frequency is small and the conductivity of the body is high. A rigorous justification of the model involves the topology of the conducting body [3]. The eddy current model is described by the system

	$\displaystyle\nabla\times\bm{E}_{\alpha}$	$\displaystyle=\mathrm{i}\omega\mu\bm{H}_{\alpha},$		(1a)
	$\displaystyle\nabla\times\bm{H}_{\alpha}$	$\displaystyle=\bm{J}_{0}+\sigma\bm{E}_{\alpha}.$		(1b)

where $\bm{E}_{\alpha}$ and $\bm{H}_{\alpha}$ are the electric and magnetic interaction fields, respectively, $\bm{J}_{0}$ is an external current source, $\mathrm{i}:=\sqrt{-1}$, $\omega$ is the angular frequency, $\mu$ is the magnetic permeability and $\sigma$ is the electric conductivity. We will use the eddy current model for describing the forward and inverse problems associated with metal detection.

Following [4, 14] we define $B_{\alpha}:=\alpha B+\bm{z}$ where $B$ is a unit size object, $\alpha$ is the object size and $\bm{z}$ is the object’s translation from the origin as shown in Figure  2.

Using the derivation in [17], we state the explicit formulae for the computation of the coefficients of $\mathcal{M}$, which are particularly well suited to a FEM discretisation. The earlier explicit expressions in [14, 15, 16] are equivalent for exact fields. We use the splitting $(\mathcal{M})_{ij}:=(\mathcal{N}^{0})_{ij}+(\mathcal{R})_{ij}+\mathrm{i}(\mathcal{I})_{ij}$ obtained in  [17] with

	$\displaystyle(\mathcal{N}^{0}[\alpha B])_{ij}$	$\displaystyle:=\alpha^{3}\delta_{ij}\int_{B}(1-\mu_{r}^{-1})\mathrm{d}\bm{\xi}+\frac{\alpha^{3}}{4}\int_{B\cup B^{c}}\tilde{\mu}_{r}^{-1}\nabla\times\tilde{\bm{\theta}}_{i}^{(0)}\cdot\nabla\times\tilde{\bm{\theta}}_{j}^{(0)}\ \mathrm{d}\bm{\xi},$		(9a)
	$\displaystyle(\mathcal{R}[\alpha B,\omega])_{ij}$	$\displaystyle:=-\frac{\alpha^{3}}{4}\int_{B\cup B^{c}}\tilde{\mu}_{r}^{-1}\nabla\times\bm{\theta}_{j}^{(1)}\cdot\nabla\times\overline{\bm{\theta}_{i}^{(1)}}\ \mathrm{d}\bm{\xi},$		(9b)
	$\displaystyle(\mathcal{I}[\alpha B,\omega])_{ij}$	$\displaystyle:=\frac{\alpha^{3}}{4}\int_{B}\nu\Big{(}\bm{\theta}_{j}^{(1)}+(\tilde{\bm{\theta}}_{j}^{(0)}+\bm{e}_{j}\times\bm{\xi})\Big{)}\cdot\Big{(}\overline{\bm{\theta}_{i}^{(1)}+(\tilde{\bm{\theta}}_{i}^{(0)}+\bm{e}_{i}\times\bm{\xi})}\Big{)}\ \mathrm{d}\bm{\xi},$		(9c)

where $\mathcal{N}^{0}[\alpha B]$, $\mathcal{R}[\alpha B,\omega]$ and $\mathcal{I}[\alpha B,\omega]$ are real symmetric rank 2 tensors, which each have real eigenvalues. In the above,

$\displaystyle\tilde{\mu}_{r}(\bm{\xi}):=\left\{\begin{array}[]{ll}\mu_{r}:=\mu_{*}/\mu_{0}&\bm{\xi}\in B\\ 1&\bm{\xi}\in B^{c}\end{array}\right.,$

(12)

and $\nu:=\alpha^{2}\omega\mu_{0}\sigma_{*}$, $\delta_{ij}$ is the Kronecker delta and the overbar denotes the complex conjugate. The computation of (9) rely on the solution of the transmission problems [17]

	$\displaystyle\nabla\times\tilde{\mu}_{r}^{-1}\nabla\times\bm{\theta}_{i}^{(0)}$	$\displaystyle=\bm{0}$	$\displaystyle\textrm{in }B\cup B^{c},$		(13a)
	$\displaystyle\nabla\cdot\bm{\theta}_{i}^{(0)}$	$\displaystyle=0$	$\displaystyle\textrm{in }B\cup B^{c},$		(13b)
	$\displaystyle[{\bm{n}}\times\bm{\theta}_{i}^{(0)}]_{\Gamma}$	$\displaystyle=\bm{0}$	$\displaystyle\textrm{on }\Gamma,$		(13c)
	$\displaystyle[{\bm{n}}\times\tilde{\mu}_{r}^{-1}\nabla\times\bm{\theta}_{i}^{(0)}]_{\Gamma}$	$\displaystyle=\bm{0}$	$\displaystyle\textrm{on }\Gamma,$		(13d)
	$\displaystyle\bm{\theta}_{i}^{(0)}-{\bm{e}}_{i}\times\bm{\xi}$	$\displaystyle=\bm{O}(|\bm{\xi}|^{-1})$	$\displaystyle\textrm{as }|\bm{\xi}|\rightarrow\infty,$		(13e)

where $\Gamma:=\partial B$ and

	$\displaystyle\nabla\times{\mu}_{r}^{-1}\nabla\times\bm{\theta}_{i}^{(1)}-\mathrm{i}\nu(\bm{\theta}_{i}^{(0)}+\bm{\theta}_{i}^{(1)})$	$\displaystyle=\bm{0}$	$\displaystyle\textrm{in }B,$		(14a)
	$\displaystyle\nabla\times\nabla\times\bm{\theta}_{i}^{(1)}$	$\displaystyle=\bm{0}$	$\displaystyle\textrm{in }B^{c},$		(14b)
	$\displaystyle\nabla\cdot\bm{\theta}_{i}^{(1)}$	$\displaystyle=0$	$\displaystyle\textrm{in }B^{c},$		(14c)
	$\displaystyle[{\bm{n}}\times\bm{\theta}_{i}^{(1)}]_{\Gamma}$	$\displaystyle=\bm{0}$	$\displaystyle\textrm{on }\Gamma,$		(14d)
	$\displaystyle[{\bm{n}}\times\tilde{\mu}_{r}^{-1}\nabla\times\bm{\theta}_{i}^{(1)}]_{\Gamma}$	$\displaystyle=\bm{0}$	$\displaystyle\textrm{on }\Gamma,$		(14e)
	$\displaystyle\bm{\theta}_{i}^{(1)}$	$\displaystyle=\bm{O}(|\bm{\xi}|^{-1})$	$\displaystyle\textrm{as }|\bm{\xi}|\rightarrow\infty.$		(14f)

Note also that we choose to introduce $\tilde{\bm{\theta}}_{i}^{(0)}\vcentcolon=\bm{\theta}_{i}^{(0)}-{\bm{e}}_{i}\times\bm{\xi}$, which can be shown to satisfy the same transmission problem as (13) except with a non-zero jump condition for $[{\bm{n}}\times\tilde{\mu}_{r}^{-1}\nabla\times\tilde{\bm{\theta}}_{i}^{(0)}]_{\Gamma}$ and the decay condition $\tilde{\bm{\theta}}_{i}^{(0)}(\bm{\xi})=\bm{O}(|\bm{\xi}|^{-1})$ as $|\bm{\xi}|\rightarrow\infty$.

Following the approach advocated in [19] for magnetostatic and eddy current problems, we add a regularisation term $\varepsilon\int_{\Omega}\tilde{\bm{\theta}}_{i}^{(0)}\cdot\bm{\psi}\mathrm{d}\bm{\xi}$, where $\varepsilon$ is a small regularisation parameter, to the weak variational statement of (13), written in terms of $\tilde{\bm{\theta}}_{i}^{(0)}$, in order to circumvent the Coulomb gauge $\nabla\cdot\tilde{\bm{\theta}}_{i}^{(0)}=0$. For details of the small error induced by this approximation see [19, 36]. Then, by choosing an appropriate set of $\bm{H}(\text{curl})$ conforming finite element functions in $W^{(hp)}\subset\bm{H}(\text{curl})$, we obtain the following discrete regularised weak form for (13) : Find real solutions $\tilde{\bm{\theta}}_{i}^{(0,hp)}\in Y^{\varepsilon}\cap W^{(hp)}$ such that

	$\displaystyle\int_{\Omega}\tilde{\mu}_{r}^{-1}\nabla\times\tilde{\bm{\theta}}_{i}^{(0,hp)}\cdot\nabla\times\bm{\psi}^{(hp)}\mathrm{d}\bm{\xi}$	$\displaystyle+\varepsilon\int_{\Omega}\tilde{\bm{\theta}}_{i}^{(0,hp)}\cdot\bm{\psi}^{(hp)}\mathrm{d}\bm{\xi}$
		$\displaystyle=2\int_{B}(1-\mu_{r}^{-1})\bm{e}_{i}\cdot\nabla\times\bm{\psi}^{(hp)}\mathrm{d}\bm{\xi},$		(16)

for all $\bm{\psi}^{(hp)}\in Y^{\varepsilon}\cap W^{(hp)}$, where

$$Y^{\varepsilon}=\Big{\{}\bm{u}\in\bm{H}(\text{curl}):{\bm{n}}\times{\bm{u}}=\bm{0}\textrm{ on }\partial\Omega\Big{\}}.$$

In a similar manner, the discrete weak variational statement of (14) is: Find complex solutions ${\bm{\theta}}_{i}^{(1,hp)}\in Y^{\varepsilon}\cap W^{(hp)}$ such that

	$\displaystyle\int_{\Omega}\big{(}\mu_{r}^{-1}\nabla\times\bm{\theta}_{i}^{(1,hp)}\big{)}$	$\displaystyle\cdot\big{(}\nabla\times\overline{\bm{\psi}^{(hp)}}\big{)}\mathrm{d}\bm{\xi}-\mathrm{i}\int_{B}\nu\bm{\theta}_{i}^{(1,hp)}\cdot\overline{\bm{\psi}^{(hp)}}\mathrm{d}\bm{\xi}$
		$\displaystyle+\varepsilon\int_{\Omega\setminus B}\bm{\theta}_{i}^{(1,hp)}\cdot\overline{\bm{\psi}^{(hp)}}\mathrm{d}\bm{\xi}=\mathrm{i}\int_{B}\nu\bm{\theta}_{i}^{(0,hp)}\cdot\overline{\bm{\psi}^{(hp)}}\mathrm{d}\bm{\xi},$		(17)

for all $\bm{\psi}^{(hp)}\in Y^{\varepsilon}\cap W^{(hp)}$ where the overbar denotes the complex conjugate.

For what follows it is beneficial to restate (3.1) in the following form: Find ${\bm{\theta}}_{i}^{(1,hp)}\in Y^{\varepsilon}\cap W^{(hp)}$ such that

$$a\big{(}\bm{\theta}_{i}^{(1,hp)},\bm{\psi}^{(hp)};\bm{\omega}\big{)}=r\big{(}\bm{\psi}^{(hp)};\bm{\theta}_{i}^{(0,hp)},\bm{\omega}\big{)},$$

(18)

for all $\bm{\psi}^{(hp)}\in Y^{\varepsilon}\cap W^{(hp)}$ where

	$\displaystyle a\big{(}\bm{\theta}_{i}^{(1,hp)},\bm{\psi}^{(hp)};\bm{\omega}\big{)}\vcentcolon$	$\displaystyle=\left<\tilde{\mu}^{-1}\nabla\times\bm{\theta}_{i}^{(1,hp)},\nabla\times{\bm{\psi}^{(hp)}}\right>_{L^{2}(\Omega)}$
		$\displaystyle-\mathrm{i}\left<\nu\bm{\theta}_{i}^{(1,hp)},{\bm{\psi}^{(hp)}}\right>_{L^{2}(B)}$
		$\displaystyle+\varepsilon\left<\bm{\theta}_{i}^{(1,hp)},\bm{\psi}^{(hp)}\right>_{L^{2}(\Omega\setminus B)},$		(19a)
	$\displaystyle r\big{(}\bm{\psi}^{(hp)};\bm{\theta}_{i}^{(0,hp)},\bm{\omega}\big{)}\vcentcolon$	$\displaystyle=\mathrm{i}\left<\nu\bm{\theta}_{i}^{(0,hp)},{\bm{\psi}^{(hp)}}\right>_{L^{2}(B)},$		(19b)

$\left<\bm{u},\bm{v}\right>_{L^{2}(\Omega)}:=\int_{\Omega}{\bm{u}}\cdot\overline{\bm{v}}\mathrm{d}\bm{\xi}$ denotes the $L^{2}$ inner product over $\Omega$ and $\bm{\omega}$ indicates the list of the problem parameters $\{\omega,\sigma_{*},\mu_{r},\alpha\}$ that one might wish to vary. Note that $r\big{(}\cdot;\cdot,\cdot\big{)}$ is a function of $\mu_{r}$ as $\bm{\theta}_{i}^{(0,hp)}$ depends on $\mu_{r}$.

For the implementation of the full order model, we use NGSolve [1, 29, 28, 36] along with the hierarchic set of $\bm{H}(\text{curl})$ conforming basis functions proposed by Schöberl and Zaglmayr [30], which are available in this software. In the following, for simplicity, we focus on the treatment of $\bm{\theta}_{i}^{(1,hp)}$ and drop the index $i$ as each direction can be computed in a similar way (as can $\tilde{\bm{\theta}}_{i}^{(0,hp)}$). We denote these basis functions with $\bm{N}^{(k)}(\bm{\xi})\in W^{(hp)}$ leading to the expression of the solution function along with the weighting functions

	$\displaystyle\bm{\theta}^{(1,hp)}(\bm{\xi},\bm{\omega})$	$\displaystyle\vcentcolon=\sum_{k=1}^{N_{d}}\bm{N}^{(k)}(\bm{\xi})\mathrm{q}_{k}(\bm{\omega}),$		(20a)
	$\displaystyle\bm{\psi}^{(hp)}(\bm{\xi},\bm{\omega})$	$\displaystyle\vcentcolon=\sum_{k=1}^{N_{d}}\bm{N}^{(k)}(\bm{\xi})\mathrm{l}_{k}(\bm{\omega}),$		(20b)

where $N_{d}$ is the number of degrees of freedom. Here, and in the following, the bold italic font denotes a vector field and the bold non-italic Roman font represents a matrix (upper case) or column vector (lower case). With this distinction, we rewrite (20) in matrix form as

	$\displaystyle\bm{\theta}^{(1,hp)}(\bm{\xi},\bm{\omega})$	$\displaystyle=\textbf{N}(\bm{\xi})\textbf{q}(\bm{\omega}),$		(21a)
	$\displaystyle\bm{\psi}^{(hp)}(\bm{\xi},\bm{\omega})$	$\displaystyle=\textbf{N}(\bm{\xi})\textbf{l}(\bm{\omega}),$		(21b)

where $\textbf{N}(\bm{\xi})$ is the matrix constructed with the basis vectors $\bm{N}^{k}(\bm{\xi})$ as its columns, i.e.

$$\textbf{N}(\bm{\xi})\vcentcolon=\big{[}\bm{N}^{(1)}(\bm{\xi}),\bm{N}^{(2)}(\bm{\xi}),...,\bm{N}^{(N_{d})}(\bm{\xi})\big{]}.$$

With this, we may also rewrite (18) as follows

$$\sum_{i=1}^{N_{d}}\sum_{j=1}^{N_{d}}\overline{l_{i}(\bm{\omega})}a\big{(}\bm{N}^{(j)}(\bm{\xi}),\bm{N}^{(i)}(\bm{\xi});\bm{\omega}\big{)}q_{j}(\bm{\omega})=\sum_{i=1}^{N_{d}}\overline{l_{i}(\bm{\omega})}r\big{(}\bm{N}^{(i)}(\bm{\xi});\bm{\theta}^{(0,hp)},\bm{\omega}\big{)},$$

(22)

and, with a suitable choice of $l_{i}(\bm{\omega})$, we may rewrite (22) as the linear system of equations

$$\textbf{A}(\bm{\omega})\textbf{q}(\bm{\omega})=\textbf{r}(\bm{\theta}^{(0,hp)},\bm{\omega}),$$

(23)

where the coefficients of $\bf{A}(\bm{\omega})$ and ${\bf r}({\bm{\theta}}^{(0,hp)},\bm{\omega})$ are defined to be

	$\displaystyle(\textbf{A}(\bm{\omega}))_{ij}$	$\displaystyle\vcentcolon=a\big{(}\bm{N}^{(j)}(\bm{\xi}),\bm{N}^{(i)}(\bm{\xi});\bm{\omega}\big{)},$		(24a)
	$\displaystyle(\textbf{r}(\bm{\theta}^{(0,hp)},\bm{\omega}))_{i}$	$\displaystyle\vcentcolon=r\big{(}\bm{N}^{(i)}(\bm{\xi});\bm{\theta}^{(0,hp)},\bm{\omega}\big{)}.$		(24b)

NGSolve offers efficient approaches for the computational solution to (23) using preconditioned iterative solvers [36, 19], which we exploit. Following the solution of (23), we can obtain $\bm{\theta}^{(1,hp)}(\bm{\xi},\bm{\omega})$ using (21) and, by repeating the process for $i=1,2,3$, we get $\bm{\theta}_{i}^{(1,hp)}(\bm{\xi},\bm{\omega})$. Then $(\mathcal{M}[\alpha B,\omega,\mu_{r},\sigma_{*}])_{ij}$, for the full order model, is found by using (9).

Following the solution of (23) for $\mathbf{q}(\bm{\omega})$ for different values of the set of parameters, $\bm{\omega}$, we construct a matrix $\mathbf{D}\in\mathbb{C}^{N_{d}\times N}$ with the vector of solution coefficients as its columns in the form

$$\mathbf{D}\vcentcolon=\big{[}\mathbf{q}(\bm{\omega}_{1}),\mathbf{q}(\bm{\omega}_{2}),...,\mathbf{q}(\bm{\omega}_{N})\big{]},$$

(25)

where $N\ll N_{d}$ denote the number of snapshots. Application of a singular value decomposition (SVD) e.g. [7, 10] gives

$$\mathbf{D}=\mathbf{U\Sigma V}^{*},$$

(26)

where $\mathbf{U}\in\mathbb{C}^{N_{d}\times N_{d}}$ and $\mathbf{V}^{*}\in\mathbb{C}^{N\times N}$ are unitary matrices and $\mathbf{\Sigma}\in\mathbb{R}^{N_{d}\times N}$ is a diagonal matrix enlarged by zeros so that it becomes rectangular. In the above, $\mathbf{V}^{*}=\overline{\mathbf{V}}^{T}$ is the hermitian of $\mathbf{V}$.

The diagonal entries $(\mathbf{\Sigma})_{ii}=\sigma_{i}$ ¹¹1Note that $\sigma_{*}$ is used for conductivity and $\sigma_{i}$ for a singular value, however, it should be clear from the application as to which definition applies are the singular values of $\mathbf{D}$ and they are arranged as $\sigma_{1}>\sigma_{2}>...>\sigma_{N}$. Based on the sparse representation of the solutions to (14) as function of $\nu$, and hence $\omega$, (and hence also the sparse representation of the MPT) found in [17], we expect these to decay rapidly towards zero, which motivates the introduction of the truncated singular value decomposition (TSVD) e.g. [7, 10]

$$\mathbf{D}\approx\mathbf{D}^{M}=\mathbf{U}^{M}\mathbf{\Sigma}^{M}(\mathbf{V}^{M})^{*},$$

(27)

where $\mathbf{U}^{M}\in\mathbb{C}^{N_{d}\times M}$ are the first $M$ columns of $\mathbf{U}$, $\mathbf{\Sigma}^{M}\in{\mathbb{R}}^{M\times M}$ is a diagonal matrix containing the first $M$ singular values and $(\mathbf{V}^{M})^{*}\in{\mathbb{C}}^{M\times N}$ are the first $M$ rows of $\mathbf{V}^{*}$. The computation of (27) constitutes the off-line stage of the POD. Using (27) we can recover an approximate representation for each of our solution snapshots as follows

$$\mathbf{q}(\bm{\omega}_{j})\approx\mathbf{U}^{M}\mathbf{\Sigma}^{M}((\mathbf{V}^{M})^{*})_{j},$$

(28)

where $((\mathbf{V}^{M})^{*})_{j}$ refers to the $j$th column of $(\mathbf{V}^{M})^{*}$.

In the online stage of PODP, $\mathbf{q}^{PODP}({\bm{\omega}})\approx\mathbf{q}({\bm{\omega}})$ is obtained by taking a linear combination of the columns of ${\bf U}^{M}$ where the coefficients of this projection are contained in the vector ${\bf p}^{M}$. We choose to also approximate $\mathbf{l}({\bm{\omega}})$ in a similar way so that

	$\displaystyle\bm{\theta}^{(1,hp)}(\bm{\xi},\bm{\omega})\approx(\bm{\theta}^{(1,hp)})^{\text{PODP}}({\bm{\xi}},\bm{\omega)}:=$	$\displaystyle\textbf{N}(\bm{\xi})\mathbf{q}^{PODP}({\bm{\omega}})=\textbf{N}(\bm{\xi})\mathbf{U}^{M}\textbf{p}^{M}(\bm{\omega})\in Y^{(PODP)},$		(29a)
	$\displaystyle\bm{\psi}^{(hp)}(\bm{\xi},\bm{\omega})\approx(\bm{\psi}^{(1,hp)})^{\text{PODP}}({\bm{\xi}},\bm{\omega)}:=$	$\displaystyle\textbf{N}(\bm{\xi})\textbf{l}^{PODP}(\bm{\omega})=\textbf{N}(\bm{\xi})\mathbf{U}^{M}\textbf{o}^{M}(\bm{\omega})\in Y^{(PODP)},$		(29b)

where $\in Y^{(PODP)}\subset Y^{\varepsilon}\cap W^{(hp)}$. Substituting these lower dimensional representations in to (22) we obtain the following

	$\displaystyle\sum_{i=1}^{M}\sum_{j=1}^{M}\overline{o^{M}_{i}(\bm{\omega})}a\big{(}\bm{N}^{(j)}(\bm{\xi})(\mathbf{U}^{M})_{j}$	$\displaystyle,\bm{N}^{(i)}(\bm{\xi})(\mathbf{U}^{M})_{i};\bm{\omega}\big{)}p^{M}_{j}(\bm{\omega})$
		$\displaystyle=\sum_{i=1}^{M}\overline{o^{M}_{i}(\bm{\omega})}r\big{(}\bm{N}^{(i)}(\bm{\xi})(\mathbf{U}^{M})_{i};\bm{\theta}^{(0,hp)},\bm{\omega}\big{)},$
	$\displaystyle({\mathbf{o}^{M}}(\bm{\omega}))^{*}((\mathbf{U}^{M})^{*}\mathbf{A}(\bm{\omega})\mathbf{U}^{M})\mathbf{p}^{M}(\bm{\omega})$	$\displaystyle=({\mathbf{o}^{M}}(\bm{\omega}))^{*}(\mathbf{U}^{M})^{*}\mathbf{r}(\bm{\theta}^{(0,hp)},\bm{\omega}).$		(30)

Then, if we choose $\mathbf{o}^{M}(\bm{\omega})$ appropriately, we obtain the linear system

$$\mathbf{A}^{M}(\bm{\omega})\mathbf{p}^{M}(\bm{\omega})=\mathbf{r}^{M}(\bm{\theta}^{(0,hp)},\bm{\omega}),$$

(31)

which is of size $M\times M$ where $\mathbf{A}^{M}(\bm{\omega})\vcentcolon=(\mathbf{U}^{M})^{*}\mathbf{A}(\bm{\omega})\mathbf{U}^{M}$ and $\mathbf{r}^{M}(\bm{\theta}^{(0,hp)},\bm{\omega})\vcentcolon=(\mathbf{U}^{M})^{*}\mathbf{r}(\bm{\theta}^{(0,hp)},\bm{\omega})$. Note, since $M<N\ll N_{d}$, this is significantly smaller than (23) and, therefore, substantially computationally cheaper to solve. After solving this reduced system, and obtaining $\mathbf{p}^{M}(\bm{\omega})$, we obtain an approximate solution for $\bm{\theta}^{(1,hp)}(\bm{\xi},\bm{\omega})$ using (29).

Focusing on the particular case where $\bm{\omega}=\omega$, from (19) we observe that we can express $\mathbf{A}$ and $\mathbf{r}$ as the simple sums

	$\displaystyle\mathbf{A}(\omega)=$	$\displaystyle\mathbf{A}^{(0)}+\omega\mathbf{A}^{(1)},$
	$\displaystyle\mathbf{r}(\bm{\theta}^{(0,hp)},\omega)=$	$\displaystyle\omega\mathbf{r}^{(1)}(\bm{\theta}^{(0,hp)}),$

where the definitions of $\mathbf{A}^{(0)}$, $\mathbf{A}^{(1)}$ and $\mathbf{r}^{(1)}(\bm{\theta}^{(0,hp)})$ are obvious from (24),(19) and the definition of $\nu$. Then, by computing and storing $(\mathbf{U}^{M})^{*}\mathbf{A}^{(0)}\mathbf{U}^{M}$, $(\mathbf{U}^{M})^{*}\mathbf{A}^{(1)}\mathbf{U}^{M}$ , $(\mathbf{U}^{M})^{*}\mathbf{r}^{(1)}(\bm{\theta}^{(0,hp)})$, which are independent of $\omega$, it follows that $\mathbf{A}^{M}({\omega})$ and $\mathbf{r}^{M}(\bm{\theta}^{(0,hp)},{\omega})$ can be efficiently calculated for each new $\omega$ from the stored data. In a similar manner, by precomputing appropriate data, the MPT coefficients in (9) can also be rapidly evaluated for each new $\omega$ using the PODP solutions. This leads to further considerable computational savings. We emphasise that the PODP is only applied to obtain ROM solutions for $\bm{\theta}^{(1)}(\bm{\xi},{\omega})$ and not to $\bm{\theta}^{(0)}(\bm{\xi})$, which does not depend on $\omega$.

We follow the approach described in [11], which enables us to derive and compute certificates of accuracy on the MPT coefficients obtained with PODP, with respect to those obtained with full order model, as a function of $\omega$. To do this, we set $\bm{\epsilon}_{i}(\omega)=\bm{\theta}_{i}^{(1,hp)}(\omega)-(\bm{\theta}_{i}^{(1,hp)})^{\text{PODP}}(\omega)\in Y^{(hp)}$, where we have reintroduced the subscript $i$, as we need to distinguish between the cases $i=1,2,3$. Although $\bm{\epsilon}_{i}$ also depends on $\bm{\xi}$, we have chosen here, and in the following, to only emphasise its dependence on $\omega$. We have also introduced $Y^{(hp)}=Y^{\varepsilon}\cap W^{(hp)}$ for simplicity of notation, and note that this error satisfies

$\displaystyle a(\bm{\epsilon}_{i}(\omega),\bm{\psi};{\omega})=r(\bm{\psi};\bm{\theta}_{i}^{(0,hp)},\omega)\qquad\forall\bm{\psi}\in Y^{(hp)},$

(32)

which is called the error equation [11] and

$\displaystyle a(\bm{\epsilon}_{i}(\omega),\bm{\psi};\omega)=0\qquad\forall\bm{\psi}\in Y^{(PODP)},$

(33)

which is called Galerkin orthogonality [11]. The Riesz representation [11] of $r(\cdot;\bm{\theta}_{i}^{(0,hp)},\omega)$ denoted by $\hat{\bm{r}}_{i}(\omega)\in Y^{(hp)}$ is such that

$\displaystyle(\hat{\bm{r}}_{i}(\omega),\bm{\psi})_{Y^{(hp)}}=r(\bm{\psi};\bm{\theta}_{i}^{(0,hp)},\omega)\qquad\forall\bm{\psi}\in Y^{(hp)},$

(34)

so that

$\displaystyle a(\bm{\epsilon}_{i}(\omega),\bm{\psi};{\omega})=(\hat{\bm{r}}_{i}(\omega),\bm{\psi})_{Y^{(hp)}}\qquad\forall\bm{\psi}\in Y^{(hp)}.$

(35)

Then, by using the alternative set of formulae for the tensor coefficients [17]

	$\displaystyle(\mathcal{R}[\alpha B,\omega])_{ij}$	$\displaystyle=-\frac{\alpha^{3}}{4}\int_{B}\nu\text{Im}(\bm{\theta}_{j}^{(1,hp)})\cdot\bm{\theta}_{i}^{(0,hp)}\mathrm{d}\bm{\xi}=-\frac{\alpha^{3}}{4}\left<\nu\text{Im}(\bm{\theta}_{j}^{(1,hp)}),\bm{\theta}_{i}^{(0,hp)}\right>_{L^{2}(B)},$		(36a)
	$\displaystyle(\mathcal{I}[\alpha B,\omega])_{ij}$	$\displaystyle=\frac{\alpha^{3}}{4}\left(\int_{B}\nu\text{Re}(\bm{\theta}_{j}^{(1,hp)})\cdot\bm{\theta}_{i}^{(0,hp)}\mathrm{d}\bm{\xi}+\int_{B}\nu\bm{\theta}_{j}^{(0,hp)}\cdot\bm{\theta}_{i}^{(0,hp)}\mathrm{d}\bm{\xi}\right)$
		$\displaystyle=\frac{\alpha^{3}}{4}\left(\left<\nu\text{Re}(\bm{\theta}_{j}^{(1,hp)}),\bm{\theta}_{i}^{(0,hp)}\right>_{L^{2}(B)}+\left<\nu\bm{\theta}_{j}^{(0,hp)},\bm{\theta}_{i}^{(0,hp)}\right>_{L^{2}(B)}\right),$		(36b)

written in terms of the full order solutions, we obtain the certificates for the tensor entries computed using PODP stated in the lemma below. Note that the formulae stated in (9) are used for the actual POD computation of $(\mathcal{R}^{PODP}[\alpha B,\omega])_{ij}$ and $(\mathcal{I}^{PODP}[\alpha B,\omega])_{ij}$, but the form in (36) is useful for obtaining certificates. Also, as $(\mathcal{N}[\alpha B])_{ij}$ is independent of $\omega$ we have $(\mathcal{N}^{0,PODP}[\alpha B])_{ij}=(\mathcal{N}^{0}[\alpha B])_{ij}$ and we write $\mathcal{M}^{PODP}[\alpha B,\omega]=\mathcal{N}^{0,PODP}[\alpha B]+\mathcal{R}^{PODP}[\alpha B,\omega]+\mathrm{i}\mathcal{I}^{PODP}[\alpha B,\omega]$ for the MPT obtained by PODP.