Direct Identification of the Continuous Relaxation Time and Frequency Spectra of Viscoelastic Materials

Abstract

Relaxation time and frequency spectra are not directly available by measurement. To determine them, an ill-posed inverse problem must be solved based on relaxation stress or oscillatory shear relaxation data. Therefore, the quality of spectra models has only been assessed indirectly by examining the fit of the experiment data to the relaxation modulus or dynamic moduli models. As the measures of data fitting, the mean sum of the moduli square errors were usually used, the minimization of which was an essential step of the identification algorithms. The aim of this paper was to determine a relaxation spectrum model that best approximates the real unknown spectrum in a direct manner. It was assumed that discrete-time noise-corrupted measurements of a relaxation modulus obtained in the stress relaxation experiment are available for identification. A modified relaxation frequency spectrum was defined as a quotient of the real relaxation spectrum and relaxation frequency and expanded into a series of linearly independent exponential functions that are known to constitute a basis of the space of square-integrable functions. The spectrum model, given by a finite series of these basis functions, was assumed. An integral-square error between the real unknown modified spectrum and the spectrum model was taken as a measure of the model quality. This index was proved to be expressed in terms of the measurable relaxation modulus at uniquely defined sampling instants. Next, an empirical identification index was introduced in which the values of the real relaxation modulus are replaced by their noisy measurements. The identification consists of determining the spectrum model that minimizes this empirical index. Tikhonov regularization was applied to guarantee model smoothness and noise robustness. A simple analytical formula was derived to calculate the optimal model parameters and expressed in terms of the singular value decomposition. A complete identification algorithm was developed. The analysis of the model smoothness and model accuracy for noisy measurements was carried out. The equivalence of the direct identification of the relaxation frequency and time spectra has been demonstrated when the time spectrum is modeled by a series of functions given by the product of the relaxation frequency and its exponential function. The direct identification concept can be applied to both viscoelastic fluids and solids; however, some limitations to its applicability have been pointed out. Numerical studies have shown that the proposed identification algorithm can be successfully used to identify Gaussian-like and Kohlrausch–Williams–Watt relaxation spectra. The applicability of this approach to determining other commonly used classes of relaxation spectra was also examined.

1. Introduction

Although the first papers concerning relaxation time and frequency spectra determination come from the late 1940s of the 20th century [1,2], the recovery of the relaxation spectrum from the measurement data is still an active area of research in rheology and the identification of time-variable viscoelastic mechanical characteristics [3,4,5,6,7,8,9,10]. Relaxation time and frequency spectra, with no direct accessible measurements, are recovered from the stress relaxation or oscillatory shear data by applying appropriate identification methods intended for determination of the spectra. The relaxation spectrum identification task is the problem of numerically solving a system of Fredholm integral equations of the first kind obtained for discrete measurements of the relaxation modulus or storage and loss modulus data. These problems are well-known to be the ill-posed inverse problems, the solutions to which, if any, are very sensitive to even small changes in the experiment data leading to arbitrarily large changes in the determined relaxation spectrum. Therefore, special stable algorithms are requisite to determine noise-robust relaxation spectrum models.

Over the last 80 years, different analytical and numerical tools have been applied to identify the relaxation spectrum. Numerous classes of algorithms were developed to determine continuous and discrete relaxation spectra models. Many theoretical papers have been devoted to the methods and algorithms for relaxation spectra determination, e.g., see [3,6,9,11,12,13,14]. In addition, experimental studies conducted for various viscoelastic materials motivated relaxation spectra models and appropriate identification algorithms, for example as seen in [4,15,16,17,18]. Reviews of these methods and algorithms can be found in many papers, for example [5,19,20] and, most recently, in [10,18].

After a few models and algorithms were derived from an application of the Post–Widder differential formula [21,22,23], many more intricate methods and models have been obtained based on the usage of the least-squares identification applied both to the relaxation modulus measurements obtained in the stress relaxation test [7,10,24,25,26,27,28] and to the measurements of the storage and loss moduli resulted from the oscillatory shear experiment [3,4,5,6,8,11,12,14,15,16,17,18]. For example, in [26,27,28], different identification algorithms were derived for the optimal regularized least-squares identification of relaxation time and frequency spectra in the classes of models defined by a finite series of different basis functions. In consequence, the quality of the spectra models was estimated by the mean sum of relaxation modulus or dynamic moduli square errors used as a measure, the minimization of which, with or without regularization, was an essential step of the identification algorithm. In some papers, e.g., [4,7,29], the pure least-squares identification was applied, while for example in [3,5,11,26,27,28] the regularized least-squares were used with various rules applied for the choice of regularization parameters to ensure the stability of the scheme and model smoothness. Recently, in [10], the best smoothed spectrum model—which reproduces the relaxation modulus measurements with a small error of the relaxation modulus model by minimizing the integral square norm of the spectrum—was found; however, here, the identification criterion is related only to the spectrum model and not to the unknown real spectrum, and the model error is assessed in terms of the measurement-available relaxation modulus.

In this paper, a new approach is proposed based on direct approximation of the real unknown relaxation time spectrum by a series of appropriately selected basis functions. It was assumed that discrete-time noise-corrupted measurements of a relaxation modulus obtained in the stress relaxation experiment are available for identification. First, a modified relaxation frequency spectrum was defined as a quotient of the real relaxation spectrum and relaxation frequency. This spectrum was expanded into a series of exponential functions forming a basis of the space of square-integrable functions [30]. Such expansion is equivalent to the expansion of the relaxation time spectrum into a series of basis functions, these being the products of the relaxation frequency and the exponential function of it. The spectra models, given by the finite series of these basis functions, were assumed. An integral square error between the real unknown spectrum and the spectrum model was taken as a measure of the model quality index. The equivalence of such defined indices for the relaxation time spectrum and the modified frequency spectrum was proved, which means an equivalence between the respective spectra approximation tasks. Next, an empirical identification index was introduced by replacing the real relaxation modulus by their noise measurements. The resulting identification problem is a linear-quadratic optimization task in which Tikhonov regularization is applied to ensure its well-posedness. Simple analytical formula for determining the optimal model parameters was derived; the singular value decomposition can be used for algebraic computations. A complete identification algorithm for determining the optimal models of the relaxation spectra has been developed. Model smoothness and noise robustness were analyzed. The results of simulation studies conducted for uni- and double-mode Gaussian-like and Kohlrausch–Williams–Watts relaxation spectra are presented. Finally, based on the congruence of the boundary conditions of the real spectra and the model basis functions, a short analysis of the applicability of the proposed approach is outlined for different classes of the real spectra, and its limitations are pointed out. It is demonstrated that the concept of direct relaxation spectrum identification can be applied both for viscoelastic fluids and viscoelastic solids. In Appendix A, the proofs and derivations of some mathematical formulas and results are given.

The idea of using a series expansion of the spectrum model has been previously applied both in the time [26,27,28] and frequency [8] domains; however, in these papers, the identification indices, being minimized, were related to the models of the relaxation or dynamic moduli and not to the unknown spectrum model. Here, the use of appropriately selected basis functions of the relaxation spectrum model allowed for linking the model quality index, related directly to the unknown spectrum, with the relaxation modulus measurements. This means that the identification index being minimized, although expressed in terms of the relaxation modulus measurements, refers directly to the unknown relaxation spectrum, not to the measured relaxation modulus. This new approach is proposed and used in this paper for the first time.

2. Materials and Methods

2.1. Relaxation Spectra

It is widely assumed in rheology [31,32,33] that the linear relaxation modulus $G\left(t\right)$ (i.e., the stress per unit strain) has a relaxation spectrum representation of the form

$$G\left(t\right)={\int }_0^{\mathrm{\infty }}{\displaystyle \frac{\mathcal{H}\left(\tau \right)}{\tau }}{e}^{-t/\tau }d\tau ,$$

(1)

or equivalently by

$$G\left(t\right)={\int }_0^{\mathrm{\infty }}{\displaystyle \frac{H\left(v\right)}{v}}{e}^{-tv}dv,$$

(2)

where the relaxation time $\mathcal{H}\left(\tau \right)$ and frequency $H\left(v\right)$ spectra, related by

$$H\left(v\right)=\mathcal{H}\left({\displaystyle \frac{1}{v}}\right), \mathcal{H}\left(\tau \right)=H\left({\displaystyle \frac{1}{\tau }}\right),$$

(3)

characterize the distributions of relaxation times $\tau$ and frequencies $v$. They are generalizations of discrete Maxwell spectra [31,32] to continuous functions of $\tau$ and $v$. Although other definitions of the relaxation spectrum are used in the literature, for example, in [34,35,36], the definition introduced by Equations (1) and (2) dominates.

2.2. Models

Following [26], the modified spectrum is introduced

$$H^M\left(v\right)={\displaystyle \frac{H\left(v\right)}{v}},$$

(4)

where the upper index of $H^M\left(v\right)$ means “modified”. Model transformation defined by (4) is a bijection. Equation (2) can be rewritten as follows:

$$G\left(t\right)={\int }_0^{\mathrm{\infty }}{H}^M\left(v\right)e^{-tv}dv,$$

(5)

i.e., the modulus $G\left(t\right)$ is directly the Laplace integral of the spectrum $H^M\left(v\right)$.

Assume that $H^M\left(v\right)\in L^2\left(0,\right.\left.\mathrm{\infty }\right)$, where $L^2\left(0,\right.\left.\mathrm{\infty }\right)$ is the space of real-valued square-integrable functions on the interval $\left(0,\right.\left.\mathrm{\infty }\right)$. The respective sufficient conditions are given by Theorem 3 in [37]. The set of the linearly independent exponential functions $\left\{e^{-\alpha kv}\right\}$, $k=\mathrm{0,1},\dots$, where $\alpha >0$, i.e., the kernel of the Laplace transformation, form a basis of the space $L^2\left(0,\right.\left.\mathrm{\infty }\right)$ [30]. Thus, the modified relaxation spectrum can be expressed as:

$$H^M\left(v\right)={\sum }_{k=0}^{\mathrm{\infty }}{g}_k{h}_k\left(v\right),$$

(6)

with basis functions defined as follows

$$h_k\left(v\right)=e^{-\alpha kv},$$

(7)

where parameter $\alpha >0$ is a time-scaling factor expressed in seconds, while $g_k$ are constant model parameters.

By (4) and (6), for the real relaxation spectrum of the material we have:

$$H\left(v\right)={\sum }_{k=0}^{\mathrm{\infty }}{g}_k{h}_k\left(v\right)v.$$

(8)

The modified spectrum $H^M\left(v\right)\in L^2\left(0,\right.\left.\mathrm{\infty }\right)$, then, $H^M\left(v\right)\to 0$ as $v\to \mathrm{\infty }$ and the first basis function can be neglected. For practical reasons, it is convenient to replace the infinite summation in the above equation with a finite one of $K$ first terms, i.e., to approximate the relaxation spectrum $H^M\left(v\right)$ (4) by a model of the form

$$H_K^M\left(v\right)={\sum }_{k=1}^K{g}_k{h}_k\left(v\right),$$

(9)

where the lower index of $H_K^M\left(v\right)$ is the number of model summands. Spectrum $H^M\left(v\right)$ (4) is expressed in Pa·s, so also Pa·s is a unit of the model’s parameters $g_k$. The model of the original spectrum $H\left(v\right)$ related to (9) takes the form

$$H_K\left(v\right)=v{H}_K^M\left(v\right)={\sum }_{k=1}^K{g}_k{h}_k\left(v\right)v.$$

(10)

The related relaxation modulus model, by (5) and (9) is described by the following:

$$G_K\left(t\right)={\int }_0^{\mathrm{\infty }}{H}_K^M\left(v\right)e^{-tv}dv={\sum }_{k=1}^K{g}_k{\displaystyle \frac{1}{t+\alpha k}}={\sum }_{k=1}^K{g}_k{\varphi }_k\left(t\right),$$

(11)

where the basis functions, expressed in ${\mathrm{s}}^{-1}$, are as follows

$${\varphi }_k\left(t\right)={\displaystyle \frac{1}{t+\alpha k}}.$$

(12)

By the second equality in (3), (8), and (7), we obtain the following series representation of the relaxation time spectrum

$$\mathcal{H}\left(\tau \right)={\sum }_{k=0}^{\mathrm{\infty }}{g}_k{e}^{-{\displaystyle \frac{\alpha k}{\tau }}}{\displaystyle \frac{1}{\tau }},$$

(13)

Omitting, as above, the first component and considering the $K$ next terms of the series (13), the relaxation spectrum $\mathcal{H}\left(\tau \right)$ can be approximated by a model of the form

$${\mathcal{H}}_K\left(\tau \right)={\sum }_{k=1}^K{g}_k{\mathcal{h}}_k\left(\tau \right),$$

(14)

where the basis functions

$${\mathcal{h}}_k\left(\tau \right)=e^{-{\displaystyle \frac{\alpha k}{\tau }}}{\displaystyle \frac{1}{\tau }}$$

(15)

and the model parameters $g_k$, expressed in Pa·s, are identical to that of model (9). By (1) and (14), the relaxation modulus model is as follows:

$$G_K\left(t\right)={\int }_0^{\mathrm{\infty }}{\displaystyle \frac{{\mathcal{H}}_K\left(\tau \right)}{\tau }}{e}^{-t/\tau }d\tau ={\sum }_{k=1}^K{g}_k{\int }_0^{\mathrm{\infty }}{\displaystyle \frac{1}{{\tau }^2}}{e}^{-{\displaystyle \frac{t+\alpha k}{\tau }}}d\tau ={\sum }_{k=1}^K{g}_k{\varphi }_k\left(t\right),$$

where the basis functions are given by (12), i.e., it is identical to the model described by (11).

2.3. Properties of the Basis Functions

A few basis functions $h_k\left(v\right)$ (7) of the model $H_K^M\left(v\right)$ (9) are shown in Figure 1 for two different values of the time-scale factor $\alpha$, while in Figure 2, the basis functions $h_k\left(v\right)v$ of the model $H_K\left(v\right)$ (10) are given. In Figure 3, the basis functions ${\mathcal{h}}_k\left(\tau \right)$ (15) of the relaxation time spectrum model ${\mathcal{H}}_K\left(\tau \right)$ (14) are demonstrated.

Figure 4 shows the hyperbolic basis functions ${\varphi }_k\left(t\right)$ (12) of the relaxation modulus model $G_K\left(t\right)$ (11). Functions ${\varphi }_k\left(t\right)$ are almost constant in time and near zero (in the considered time intervals) for $k=100$ and $k=1000$. However, for smaller indices $k$, they are in good agreement with the real relaxation modulus obtained in the stress relaxation test.

The basis functions of the models $H_K^M\left(v\right)$, $H_K\left(v\right)$, and $G_K\left(t\right)$ are positive definite. Functions $h_k\left(v\right)$ (7) are monotonically decreasing, while basis functions $h_k\left(v\right)v$ of the model $H_K\left(v\right)$ (10) have global maxima equal to $1/\left(\alpha ke\right)$ for the relaxation frequencies $v=1/\left(\alpha k\right)$. In addition, basis functions ${\mathcal{h}}_k\left(\tau \right)$ (15) have unique maxima equal to $1/\left(\alpha ke\right)$ at the relaxation times $\tau =\alpha k$. Relaxation modulus basis functions ${\varphi }_k\left(t\right)$ (12) are monotonically decreasing. For large $v$, and in particular for $v\to \mathrm{\infty }$, the basis functions $h_k\left(v\right)$ and $h_k\left(v\right)v$ decrease exponentially to zero. Similarly, for $\tau \to \mathrm{\infty }$ functions ${\mathcal{h}}_k\left(\tau \right)\to 0$, faster than the exponential function. Functions ${\varphi }_k\left(t\right)$ tends to zero hyperbolically as $t\to \mathrm{\infty }$.

3. Results

In this section, the problem of optimal spectrum approximation in the class of models defined by a finite series of the introduced basis functions is formulated and solved. First, it is demonstrated that the problem of the optimal—in the sense of an integral square error—approximation of the modified relaxation frequency spectrum is equivalent to the problem of the optimal approximation of the relaxation time spectrum. It is also proved that, due to the choice of exponential basis functions, the integral square model error can be expressed in terms of the relaxation modulus in the sampling points uniquely determined by the basis functions of the spectrum model. Next, an empirical identification index is introduced in which the values of the real relaxation modulus are replaced by their noisy, in general, measurements. The optimal models of the relaxation time and frequency spectra are determined by solving the linear-quadratic identification task. However, this problem turned out to be ill-conditioned. Therefore, Tikhonov regularization is applied resulting in the stable, noise-robust and simple identification rule. Next, the equations and functions essential for the proposed identification scheme are described in terms of the singular value decomposition of the basic matrix of this identification problem. Model smoothness is estimated, error of the relaxation modulus is evaluated, noise robustness and convergence analysis is conducted. The complete identification algorithm is presented. The results of simulation studies for Gaussian-like and Kohlrausch–Williams–Watts relaxation spectra describing many real materials are presented. Finally, based on the congruence of the boundary conditions of the real spectra and the basis functions, a rough applicability analysis of the proposed approach is outlined for different classes of the real spectra.

3.1. Spectrum Approximation

Model $H_K^M\left(v\right)$ (9) approximates the modified spectrum $H^M\left(v\right)$ (4). As a measure of the model (9) accuracy the integral square index is taken

$$\mathcal{J}\left({\mathit{g}}_K\right)={\int }_0^{\mathrm{\infty }}{\left[H^M\left(v\right)-H_K^M\left(v\right)\right]}^{2}dv,$$

(16)

where ${\mathit{g}}_K={\left[\begin{array}{ccc}{g}_1& \cdots & g_K\end{array}\right]}^T$ is an $K$—element vector of the model (9) parameters; superscript “$T$” indicates transpose. In view of the additive form of the model (9), index $\mathcal{J}\left({\mathit{g}}_K\right)$ can be expressed as follows:

$$\mathcal{J}\left({\mathit{g}}_K\right)={\int }_0^{\mathrm{\infty }}{H}^M{\left(v\right)}^{2}dv-2{\sum }_{k=1}^K{g}_k{\int }_0^{\mathrm{\infty }}{H}^M\left(v\right)h_k\left(v\right)dv+{\sum }_{k=1}^K{\sum }_{m=1}^K{g}_k{g}_m{\phi }_{km},$$

(17)

where, by (7), the coefficients

$${\phi }_{km}={\int }_0^{\mathrm{\infty }}{h}_k\left(v\right)h_m\left(v\right)dv={\int }_0^{\mathrm{\infty }}{e}^{-\alpha \left(k+m\right)v}dv={\displaystyle \frac{1}{\alpha \left(k+m\right)}}.$$

(18)

Since, in view of (5) and (7), we have:

$${\int }_0^{\mathrm{\infty }}{H}^M\left(v\right)h_k\left(v\right)dv={\int }_0^{\mathrm{\infty }}{H}^M{e}^{-\alpha kv}dv=G\left(\alpha k\right),$$

the above index takes the form

$$\mathcal{J}\left({\mathit{g}}_K\right)={\int }_0^{\mathrm{\infty }}{H^M\left(v\right)}^{2}dv-2{\sum }_{k=1}^K{g}_{k}G\left(\alpha k\right)+{\sum }_{k=1}^K{\sum }_{m=1}^K{g}_k{g}_m{\phi }_{km}.$$

(19)

Model ${\mathcal{H}}_K\left(\tau \right)$ (14) approximates the real spectrum $\mathcal{H}\left(\tau \right)$. As a measure of the model (14) accuracy the square index, analogous to (16), is taken

$$\mathcal{J}\left({\mathit{g}}_K\right)={\int }_0^{\mathrm{\infty }}{\left[\mathcal{H}\left(\tau \right)-{\mathcal{H}}_K\left(\tau \right)\right]}^{2}d\tau ,$$

(20)

where ${\mathit{g}}_K$ is a vector of the model (14) parameters. In view of (14) and (15), the above index can be expressed as follows:

$$\begin{gathered} \mathcal{J}\left({\mathit{g}}_K\right)={\int }_0^{\mathrm{\infty }}{\mathcal{H}}^2\left(\tau \right)d\tau -2{\sum }_{k=1}^K{g}_k{\int }_0^{\mathrm{\infty }}\mathcal{H}\left(\tau \right){\mathcal{h}}_k\left(\tau \right)d\tau + \\ {\sum }_{k=1}^K{{\sum }_{m=1}^K{g}_{k}g}_m{\int }_0^{\mathrm{\infty }}{\mathcal{h}}_k\left(\tau \right){\mathcal{h}}_m\left(\tau \right)d\tau , \end{gathered}$$

(21)

where, by (15) and (1)

$${\int }_0^{\mathrm{\infty }}\mathcal{H}\left(\tau \right){\mathcal{h}}_k\left(\tau \right)d\tau ={\int }_0^{\mathrm{\infty }}\mathcal{H}\left(\tau \right)e^{-{\displaystyle \frac{\alpha k}{\tau }}}{\displaystyle \frac{1}{\tau }}d\tau =G\left(\alpha k\right),$$

while the integrals

$${\int }_0^{\mathrm{\infty }}{\mathcal{h}}_k\left(\tau \right){\mathcal{h}}_m\left(\tau \right)d\tau ={\int }_0^{\mathrm{\infty }}{{\displaystyle \frac{1}{{\tau }^2}}e}^{-{\displaystyle \frac{\alpha \left(k+m\right)}{\tau }}}d\tau ={\displaystyle \frac{1}{\alpha \left(k+m\right)}}={\phi }_{km},$$

(22)

are identical to that given by (18). Combining the above with (21) yields

$$\mathcal{J}\left({\mathit{g}}_K\right)={\int }_0^{\mathrm{\infty }}{\mathcal{H}}^2\left(\tau \right)d\tau -2{\sum }_{k=1}^K{g}_{k}G\left(\alpha k\right)+{\sum }_{k=1}^K{\sum }_{m=1}^K{g}_k{g}_m{\phi }_{km}.$$

For the material relaxation spectra $\mathcal{H}\left(\tau \right)$ and $H^M\left(v\right)$, by (4) and simple substitution based on (3), we have the following:

$${\int }_0^{\mathrm{\infty }}{H^M\left(v\right)}^{2}dv={\int }_0^{\mathrm{\infty }}{\displaystyle \frac{{H\left(v\right)}^2}{v^2}}dv={\int }_0^{\mathrm{\infty }}{\mathcal{H}\left({\displaystyle \frac{1}{v}}\right)}^2{\displaystyle \frac{1}{v^2}}dv={\int }_0^{\mathrm{\infty }}{\mathcal{H}\left(\tau \right)}^{2}d\tau ,$$

(23)

that is, the integral square index defined by (20) is identical to that defined by (16); therefore, the same notation was used.

3.2. Identification Problem

Suppose that a certain identification experiment (stress relaxation test [1,33,38]) performed on the specimen of the material under investigation resulted in a set of measurements of the relaxation modulus $\left\{\overline{G}\left(t_k\right)=G\left(t_k\right)+z\left(t_k\right)\right\}$ at the sampling instants $t_k=\alpha k$, $k=1,\dots ,K$; where $z\left(t_k\right)$ is additive measurement noise. Generally, identification consists of selecting within the given class of models, which ensures the best fit to the measurement results. Classically, the mean square identification index related to the measurements of the relaxation modulus $\overline{G}\left(t_k\right)$ is used; compare [7,24,26,27,28]. This means that the model quality index is not related directly to the unknown relaxation spectrum, which is inaccessible by measurement, but to the measurement-available relaxation modulus. Such an approach is typical in the context of the inverse problem.

Here, as a measure of the model (9), equivalently (14), the accuracy of the index $\mathcal{J}\left({\mathit{g}}_K\right)$ of the form (16) is applied. Note, that the first component of $\mathcal{J}\left({\mathit{g}}_K\right)$ given by the right-hand side of (19) depends on the unknown relaxation spectrum; the second term is determined by model parameters $g_k$ and the values of the measurable relaxation modulus at the time instant $t_k=\alpha k$; and the last component is affected only by the model parameters and the times $t_k$. Replacing in (19) the relaxation modulus $G\left(t_k\right)=G\left(\alpha k\right)$ by their measurements $\overline{G}\left(t_k\right)$, we obtain the following integral-empirical index

$${\overline{\mathcal{J}}}_K\left({\mathit{g}}_K\right)={\int }_0^{\mathrm{\infty }}{H^M\left(v\right)}^{2}dv-2{\sum }_{k=1}^K{g}_k\overline{G}\left(\alpha k\right)+{\sum }_{k=1}^K{\sum }_{m=1}^K{g}_k{g}_m{\phi }_{km}.$$

(24)

Let us introduce the vector-matrix notation

$${\mathit{\Phi }}_K=\left[\begin{array}{ccc}{\displaystyle \frac{1}{2}}& \cdots & {\displaystyle \frac{1}{K+1}}\\ ⋮& \ddots & ⋮\\ {\displaystyle \frac{1}{K+1}}& \cdots & {\displaystyle \frac{1}{2K}}\end{array}\right], {\overline{\mathit{G}}}_K=\left[\begin{array}{c}\overline{G}\left(t_1\right)\\ ⋮\\ \overline{G}\left(t_K\right)\end{array}\right].$$

(25)

Note that element $(k,m)$ of the matrix ${\mathit{\Phi }}_K$, i.e., the entry in the $k$-th row and $m$-th column of ${\mathit{\Phi }}_K$, is equal to $\alpha {\phi }_{km}$ and is dimensionless. Therefore, the algebraic properties of the matrix ${\mathit{\Phi }}_K$ do not depend on the time-scale factor $\alpha$. Using the above notation and bearing in mind (18), the identification index (24) can be expressed in compact form as follows:

$${\overline{\mathcal{J}}}_K\left({\mathit{g}}_K\right)={\int }_0^{\mathrm{\infty }}{H^M\left(v\right)}^{2}dv-2{\overline{\mathit{G}}}_K^T{\mathit{g}}_K+{\displaystyle \frac{1}{\alpha }}{\mathit{g}}_K^T{\mathit{\Phi }}_K{\mathit{g}}_K.$$

(26)

Thus, the optimal identification of the relaxation spectrum in the class of models (9), equivalently (14), consists of solving—with respect to the model parameter ${\mathit{g}}_K$—the linear-quadratic problem

$$\underset{{\mathit{g}}_K\in R^K}{min}{\overline{\mathcal{J}}}_K\left({\mathit{g}}_K\right),$$

(27)

the system of normal equations of which is as follows:

$${\mathit{\Phi }}_K{\mathit{g}}_K=\alpha {\overline{\mathit{G}}}_K.$$

(28)

The existence and properties of the solution to (28) depend on the properties of square symmetric matrix ${\mathit{\Phi }}_K$ (25) specified by the following result, which is the simple consequence of the independence of the basis functions $h_k\left(v\right)$ (7), as proved in Appendix A.1.

Lemma 1. 

The matrix ${\mathit{\Phi }}_K$, Equation (25), is positive definite for an arbitrary $K\ge 1$.

By the above lemma, the unique solution of the minimization task (27) is as follows:

$${\mathit{g}}_K=\alpha {{\mathit{\Phi }}_K^{-1}\overline{\mathit{G}}}_K.$$

(29)

The matrix ${\mathit{\Phi }}_K$, although of full-rank, is extremely ill-conditioned and must be used with care. Namely, the problem of the matrix ${\mathit{\Phi }}_K$ inversion is ill-conditioned, therefore small perturbations in ${\mathit{\Phi }}_K$ may produce large changes in ${\mathit{\Phi }}_K^{-1}$. The spectral condition number ([39] Equation (2.6.3))

$$\kappa \left({\mathit{\Phi }}_K\right)={‖{\mathit{\Phi }}_K^{-1}‖}_2·{‖{\mathit{\Phi }}_K‖}_2,$$

(30)

where ${‖·‖}_2$ is the spectral norm of matrices, and equal, in fact, to the ratio of the largest singular value of ${\mathit{\Phi }}_K$ to the smallest, measures the sensitivity of the answer to small perturbation of the data. From the first row values in

Table 1. Spectral condition numbers $\kappa \left({\mathit{\Phi }}_K\right)$ and $\kappa \left({\mathit{\Phi }}_K+\alpha \lambda {\mathbb{I}}_{K,K}\right)$ defined according to (30).

$\mathit{\alpha }\mathit{\lambda }$	$\mathit{K}$ = 2	$\mathit{K}$ = 3	$\mathit{K}$ = 4	$\mathit{K}$ = 5	$\mathit{K}$ = 8	$\mathit{K}$ = 10	$\mathit{K}$ = 15	$\mathit{K}$ = 20	$\mathit{K}$ = 100	$\mathit{K}$ = 1000	$\mathit{K}$ = 10,000
0	38.474	1.35 × 103	4.59 × 104	1.54 × 106	5.64 × 1010	6.23 × 1013	2.61 × 1017	6.45 × 1018	1.20 × 1019	9.59 × 1020	1.0 × 1022
10	1.071	1.087	1.098	1.106	1.122	1.129	1.14129	1.1496	1.18809	1.225	1.248
1	1.699	1.874	1.977	2.056	2.215	2.288	2.412	2.495	2.880	3.247	3.479
0.1	6.983	9.688	10.773	11.559	13.154	13.877	15.123	15.953	19.800	23.469	25.798
0.01	25.552	83.135	98.55	106.59	122.542	129.766	142.232	150.535	189.001	225.687	248.977
1 × 10−3	36.600	532.056	958.14	1.06 × 103	1.22 × 103	1.29 × 103	1.41 × 103	1.49 × 103	1.88 × 103	2.25 × 103	2.48 × 103
1 × 10−4	38.278	1.17 × 103	8.06 × 103	1.05 × 104	1.22 × 104	1.29 × 104	1.41 × 104	1.49 × 104	1.88 × 104	2.25 × 104	2.48 × 104
1 × 10−5	38.454	1.33× 103	3.12 × 104	9.88 × 104	1.22 × 105	1.29 × 105	1.41 × 105	1.49 × 105	1.88 × 105	2.25 × 105	2.48 × 105

, where $\kappa \left({\mathit{\Phi }}_K\right)$ is given for a few values of $K$, we see that index $\kappa \left({\mathit{\Phi }}_K\right)$ exceeds the value of 10⁵ already for $K=5$, the value of 10¹⁰ as early as $K=8$ and tends to infinity with growing $K$. Thus, positive definite ${\mathit{\Phi }}_K$ is suspected to be very ill-conditioned, even for not very large $K\ge 5$, and numerical solution of (28) results in fluctuations of the parameters vector ${\mathit{g}}_K$, which are the greater, the greater is the value of the condition number $\kappa \left({\mathit{\Phi }}_K\right)$.

Summarizing, the linear-quadratic identification task (27) is ill-conditioned [40] and when the data are noisy even small changes in ${\overline{\mathit{G}}}_K$ would lead to an arbitrarily large artefact in ${\mathit{g}}_K$ given by (29). Therefore, the numerical solution of the finite-dimensional problem (27) is fraught with the same difficulties that the original continuous ill-posed problems of numerical solution of the Fredholm Equations (2) or (5).

3.3. Regularization

To deal with the ill-conditioning, we use Tikhonov regularization [41], which is classical and because of its simplicity, probably the most common method for solving ill-posed linear-quadratic problems. For the linear-quadratic task (27), Tikhonov regularization strives in minimizing a modified square functional of the form

$${\overline{\mathcal{J}}}_K\left({\mathit{g}}_K\right)+\lambda {\mathit{g}}_K^T{\mathit{g}}_K,$$

(31)

where $\lambda >0$ is a regularization parameter. The unit of $\lambda$ must be ${\mathrm{s}}^{-1}$ to ensure dimensional consistency of the above index. For ${\overline{\mathcal{J}}}_K\left({\mathit{g}}_K\right)$ given by (26), the regularized task results in the linear-quadratic optimization task

$$\underset{{\mathit{g}}_K\in R^K}{min}{\displaystyle \frac{1}{\alpha }}{\mathit{g}}_K^T{\mathit{\Phi }}_K{\mathit{g}}_K-2{\overline{\mathit{G}}}_K^T{\mathit{g}}_K+\lambda {\mathit{g}}_K^T{\mathit{g}}_K;$$

(32)

the first summand of (26), being independent on ${\mathit{g}}_K$, does not have to be taken into account here, just as it did not affect the minimization result in the original problem (27). The set of normal equations is now as follows

$$\left({\mathit{\Phi }}_K+\alpha \lambda {\mathbb{I}}_K\right){\mathit{g}}_K=\alpha {\overline{\mathit{G}}}_K,$$

(33)

where ${\mathbb{I}}_K$ is $K\times K$ identity matrix.

The existence and properties of the solution of (33) depend on the properties of the symmetric matrix $\left({\mathit{\Phi }}_K+\alpha \lambda {\mathbb{I}}_K\right)$. Based on Lemma 1, $\left({\mathit{\Phi }}_K+\alpha \lambda {\mathbb{I}}_K\right)$ is non-singular and positive definite for any $\lambda \ge 0$. In successive rows of Table 1, the spectral conditional numbers $\kappa \left({\mathit{\Phi }}_K+\alpha \lambda {\mathbb{I}}_K\right)$ are given for a few values of the dimensionless product $\alpha \lambda$, which determines the value of $\kappa \left({\mathit{\Phi }}_K+\alpha \lambda {\mathbb{I}}_K\right)$ for given $K$. The numerical studies indicate that $\kappa \left({\mathit{\Phi }}_K+\alpha \lambda {\mathbb{I}}_K\right)$ is not greater than the numerically acceptable value 10⁵ [42] for $K\le 50$, whenever the parameters product $\alpha \lambda \ge 1.7\times {10}^{-5}$, for $K\le 100$ if $\alpha \lambda \ge 1.89\times {10}^{-5}$, and for $\alpha \lambda \ge 2.16\times {10}^{-5}$, assuming measurement points $K\le 500$.

Therefore, the problem (32) is well-posed, that is the solution exists, is unique, and continuously depends on both the matrix ${\mathit{\Phi }}_K$ and the measurements ${\overline{\mathit{G}}}_K$. By (33), the optimal regularized vector is given by the following formula:

$${\overline{\mathit{g}}}_K\left(\lambda \right)=\alpha {\left({\mathit{\Phi }}_K+\alpha \lambda {\mathbb{I}}_K\right)}^{-1}{\overline{\mathit{G}}}_K.$$

(34)

Elegant and compact Formula (34) is, however, unsuitable for computational purposes, for which the singular decomposition technique [39] will be used.

3.4. Algebraic Background

Let the singular value decomposition (SVD) of the matrix ${\mathit{\Phi }}_K$ (25) take the form [39]

$${\mathit{\Phi }}_K={\mathit{U}}_K{\mathit{\Sigma }}_K{\mathit{U}}_K^T,$$

(35)

where the diagonal $K\times K$ matrix

$${\mathit{\Sigma }}_K=diag\left({\mathit{\sigma }}_1,\dots ,{\mathit{\sigma }}_K\right),$$

(36)

is composed of the non-zero singular values ${\sigma }_1\ge \dots \ge {\sigma }_k\ge \dots \ge {\sigma }_K$ of the matrix ${\mathit{\Phi }}_K$; matrix ${\mathit{U}}_K\in R^{K,K}$ is orthogonal. The SVD (35) is uniquely determined. The singular values ${\sigma }_k$ of ${\mathit{\Phi }}_K$ do not depend on the time-scale factor $\alpha$. Therefore, for given $K$, the SVD must be computed only once even if the sampling instants $t_k=\alpha k$ dependent on the parameter $\alpha$ are changed in the experiment. Taking advantage of the diagonal structure of ${\mathit{\Sigma }}_K$ and orthogonality of the matrix ${\mathit{U}}_K$, we have:

$${\left({\mathit{\Phi }}_K+\alpha \lambda {\mathbb{I}}_K\right)}^{-1}={{\mathit{U}}_K\left({\mathit{\Sigma }}_K+\alpha \lambda {\mathbb{I}}_K\right)}^{-1}{\mathit{U}}_K^T={\mathit{U}}_K{\mathit{\Omega }}_K{{\mathit{U}}_K^T}^T,$$

(37)

where diagonal $K\times K$ matrix

$${{\mathit{\Omega }}_K=\left({\mathit{\Sigma }}_K+\alpha \lambda {\mathbb{I}}_K\right)}^{-1}=diag\left({\displaystyle \frac{1}{{\sigma }_1+\alpha \lambda }},\dots ,{\displaystyle \frac{1}{{\sigma }_K+\alpha \lambda }}\right).$$

(38)

In view of (37), (34), and (38), the optimal regularized vector ${\overline{\mathit{g}}}_{\mathit{K}}\left(\lambda \right)$ is expressed as follows:

$${\overline{\mathit{g}}}_K\left(\lambda \right)=\alpha {{\mathit{U}}_K\left({\mathit{\Sigma }}_K+\alpha \lambda {\mathbb{I}}_K\right)}^{-1}{\mathit{Y}}_K=\alpha {\mathit{U}}_K{\mathit{\Omega }}_K{\mathit{Y}}_K,$$

(39)

where the $K$ dimensional vector

$${\mathit{Y}}_K={\mathit{U}}_K^T{\overline{\mathit{G}}}_K.$$

(40)

According to (9), (10), and (14), the resulting best relaxation spectra models are as follows

$${\overline{H}}_K^M\left(v\right)={\sum }_{k=1}^K{\overline{g}}_k\left(\lambda \right)h_k\left(v\right),$$

(41)

$${\overline{H}}_K\left(v\right)={\sum }_{k=1}^K{\overline{g}}_k\left(\lambda \right)h_k\left(v\right)v,$$

(42)

and

$${\overline{\mathcal{H}}}_K\left(\tau \right)={\sum }_{k=1}^K{\overline{g}}_k\left(\lambda \right){\mathcal{h}}_k\left(\tau \right),$$

(43)

where ${\overline{g}}_k\left(\lambda \right)$ are elements of the vector ${\overline{\mathit{g}}}_K\left(\lambda \right)$ (34), or equivalent (39).

3.5. Analysis

The model smoothing and its accuracy in the case of noisy measurements of the relaxation modulus will be now analyzed. Contrary to the previous papers [26,27,28]—in which the model quality index refers to the relaxation modulus but not directly to the unknown relaxation spectrum—here we can estimate the spectra errors ${‖H^M\left(v\right)-{\overline{H}}_K^M\left(v\right)‖}_2$, ${‖{\overline{H}}_K\left(v\right)-{\tilde{H}}_K\left(v\right)‖}_2$, and ${‖\mathcal{H}\left(\tau \right)-{\overline{\mathcal{H}}}_K\left(\tau \right)‖}_2$ directly, where ${‖·‖}_2$ denotes the square norm in the space $L^2\left(0,\right.\left.\infty \right)$.

3.5.1. Model Smoothness

The purpose of the regularization applied in (31) relies on stabilization of the vector of model parameters ${\overline{\mathit{g}}}_K\left(\lambda \right)$ (34). The norms ${‖{\overline{H}}_K^M\left(v\right)‖}_2$, ${‖{\overline{H}}_K\left(v\right)‖}_2$, and ${‖{\overline{\mathcal{H}}}_K\left(\tau \right)‖}_2$ are natural measures of the spectra models’ (41)–(43) smoothness. In Appendix A.2, the following result is derived.

Proposition 1. 

Let the time-scale factor $\alpha >0$ and the regularization parameter $\lambda >0$. Then, for the optimal relaxation spectra models ${\overline{H}}_K^M\left(v\right)$ (41), ${\overline{H}}_K\left(v\right)$ (42), and ${\overline{\mathcal{H}}}_K\left(\tau \right)$ (43) we have the following

$$\begin{gathered} {{\displaystyle \frac{1}{\alpha }}{\sigma }_K{\overline{\mathit{g}}}_K^T\left(\lambda \right){\overline{\mathit{g}}}_K\left(\lambda \right)\le ‖{\overline{H}}_K^M\left(v\right)‖}_2^2={‖{\overline{\mathcal{H}}}_K\left(\tau \right)‖}_2^2={\displaystyle \frac{1}{\alpha }}{\overline{\mathit{g}}}_K^T\left(\lambda \right){\mathit{\Phi }}_K{\overline{\mathit{g}}}_K\left(\lambda \right)\le \\ {\displaystyle \frac{1}{\alpha }}{\sigma }_1{\overline{\mathit{g}}}_K^T\left(\lambda \right){\overline{\mathit{g}}}_K\left(\lambda \right), \end{gathered}$$

(44)

where the vector ${\overline{\mathit{g}}}_K\left(\lambda \right)$ is given by (39) and

$${\displaystyle \frac{2}{{\alpha }^3}}{\varsigma }_K{\overline{\mathit{g}}}_K^T\left(\lambda \right){\overline{\mathit{g}}}_K\left(\lambda \right)\le {‖{\overline{H}}_K\left(v\right)‖}_2^2={\displaystyle \frac{2}{{\alpha }^3}}{\overline{\mathit{g}}}_K^T\left(\lambda \right){\mathit{\Theta }}_K{\overline{\mathit{g}}}_K\left(\lambda \right)\le {\displaystyle \frac{2}{{\alpha }^3}}{\varsigma }_1{\overline{\mathit{g}}}_K^T\left(\lambda \right){\overline{\mathit{g}}}_K\left(\lambda \right),$$

(45)

where ${\sigma }_K$ and ${\sigma }_1$ are the minimal and maximal singular values of the matrix ${\mathit{\Phi }}_K$ (25), while ${\varsigma }_K$ and ${\varsigma }_1$ are the minimal and maximal singular values of the positive definite matrix

$${\mathit{\Theta }}_K=\left[\begin{array}{ccc}{\displaystyle \frac{1}{2^3}}& \cdots & {\displaystyle \frac{1}{{\left(K+1\right)}^3}}\\ ⋮& \ddots & ⋮\\ {\displaystyle \frac{1}{{\left(K+1\right)}^3}}& \cdots & {\displaystyle \frac{1}{{\left(2K\right)}^3}}\end{array}\right].$$

(46)

The values of square roots of the smallest and largest singular values ${\sigma }_K$, ${\sigma }_1$, ${\varsigma }_K$, and ${\varsigma }_1$ for some model summands $K$ are summarized in

Table 2. The square roots of the largest ${\sigma }_1$, ${\varsigma }_1$ and minimal ${\sigma }_K$, ${\varsigma }_K$ singular value of the matrices ${\mathit{\Phi }}_K$ (25) and ${\mathit{\Theta }}_K$ (46) for $K$ model summands.

$\mathit{K}$	10	50	100	500	1000	5000	10,000
$\sqrt{{\sigma }_1}$	1.1348	1.3153	1.3711	1.4677	1.4989	1.5556	1.5747
$\sqrt{{\sigma }_K}$	1.437 × 10−7	1.919 × 10−10	3.955 × 10−10	1.923 × 10−10	4.838 × 10−11	2.323 × 10−11	1.573 × 10−11
$\sqrt{{\varsigma }_1}$	0.3737	0.3737	0.3737	0.3737	0.3737	0.3737	0.3737
$\sqrt{{\varsigma }_K}$	6.829 × 10−8	1.151 × 10−11	1.847 × 10−12	1.017 × 10−13	2.455 × 10−14	2.769 × 10−15	6.473 × 10−15

. Due to the ill-conditioning of matrices ${\mathit{\Phi }}_K$ and ${\mathit{\Theta }}_K$, the lower bounds in (44) and (45) are not too useful. Since $\sqrt{{\sigma }_1}$ and $\sqrt{{\varsigma }_1}$ grows with $K$, from the analysis of Table 2 data and the right inequalities in (44) and (45), the next result follows immediately.

Proposition 2. 

Let the time-scale factor $\alpha >0$ and the regularization parameter $\lambda >0$. If the number of relaxation modulus measurements $K\le {10}^4$, then the optimal relaxation spectra models ${\overline{H}}_K^M\left(v\right)$ (41), ${\overline{H}}_K\left(v\right)$ (42), and ${\overline{\mathcal{H}}}_K\left(\tau \right)$ (43) are such that

$${‖{\overline{H}}_K^M\left(v\right)‖}_2={‖{\overline{\mathcal{H}}}_K\left(\tau \right)‖}_2\le {\displaystyle \frac{1.5747}{\sqrt{\alpha }}}{‖{\overline{\mathit{g}}}_K\left(\lambda \right)‖}_2,$$

(47)

$${‖{\overline{H}}_K\left(v\right)‖}_2\le {\displaystyle \frac{0.5285}{\alpha \sqrt{\alpha }}}{‖{\overline{\mathit{g}}}_K\left(\lambda \right)‖}_2,$$

(48)

where the vector ${\overline{\mathit{g}}}_K\left(\lambda \right)$ (39) and ${‖·‖}_2$ denotes here the square norm in Euclidean space $R^K$.

Since, by the right equality in (39) and the orthogonality of ${\mathit{U}}_K$, we obtain:

$${\overline{\mathit{g}}}_K^T\left(\lambda \right){\overline{\mathit{g}}}_K\left(\lambda \right)={\alpha }^2{\mathit{Y}}_K^T{{\mathit{\Omega }}_K{\mathit{U}}_K^T\mathit{U}}_K{\mathit{\Omega }}_K{\mathit{Y}}_K^T={\alpha }^2{\mathit{Y}}_K^T{\mathit{\Omega }}_K^2{\mathit{Y}}_K^T,$$

bearing in mind the diagonal structure of ${\mathit{\Omega }}_K$ (38) we have the formula:

$${{‖{\overline{\mathit{g}}}_K\left(\lambda \right)‖}_2^2=\overline{\mathit{g}}}_K^T\left(\lambda \right){\overline{\mathit{g}}}_K\left(\lambda \right)={\alpha }^2{\sum }_{k=1}^K {\displaystyle \frac{y_k^2}{{\left({\sigma }_k+\alpha \lambda \right)}^2}},$$

(49)

where $y_k$ are the elements of the vector ${\mathit{Y}}_K$ (40), which directly illustrates the mechanism of the regularization. The following rule holds: the greater the regularization parameter $\lambda$ is, the more highly bounded the fluctuations of the vector ${\overline{\mathit{g}}}_K\left(\lambda \right)$ are.

Summarizing, the smoothness of the optimal vector ${\overline{\mathit{g}}}_K\left(\lambda \right)$ of the model’s parameters guarantees that the fluctuations of the resulting relaxation spectra models ${\overline{H}}_K^M\left(v\right)$, ${\overline{\mathcal{H}}}_K\left(\tau \right)$, and ${\overline{H}}_K\left(v\right)$ are also bounded. Both time-scale factor $\alpha$ and the regularization parameter $\lambda$ affect the smoothness of the spectrum models. However, it should be remembered that the inequalities (47) and (48) give only the upper bounds of the respective norms.

3.5.2. Noise Robustness and Convergence

The model of the modified spectrum that we would obtain for the same time-scale factor $\alpha$ and regularization parameter $\lambda$ on the basis of ideal (noise-free) relaxation modulus measurements:

$${\tilde{H}}_K^M\left(v\right)={\sum }_{k=1}^K{\tilde{g}}_k\left(\lambda \right)h_k\left(v\right),$$

(50)

where ${\tilde{\mathit{g}}}_K\left(\lambda \right)$ is the vector model parameters given by the following (compare (39) and (40)):

$${\tilde{\mathit{g}}}_K\left(\lambda \right)=\alpha {\mathit{U}}_K{\mathit{\Omega }}_K{\mathit{U}}_K^T{\mathit{G}}_K$$

(51)

for the noise-free relaxation modulus ${\mathit{G}}_K={\left[\begin{array}{ccc}G\left(t_1\right)& \cdots & G\left(t_K\right)\end{array}\right]}^T$, which will be considered as a reference point for the model ${\overline{H}}_K^M\left(v\right)$ (41). The respective noise-free optimal regularized models of the relaxation frequency and time spectra are as follows

$${\tilde{H}}_K\left(v\right)={\sum }_{k=1}^K{\tilde{g}}_k\left(\lambda \right)h_k\left(v\right)v,$$

(52)

and

$${\tilde{\mathcal{H}}}_K\left(\tau \right)={\sum }_{k=1}^K{\tilde{g}}_k\left(\lambda \right){\mathcal{h}}_k\left(\tau \right).$$

(53)

In Appendix A.3, the following estimations are derived.

Proposition 3. 

For an arbitrary time-scale factor $\alpha$ and arbitrary regularization parameter $\lambda$, the errors between the relaxation spectra models ${\overline{H}}_K^M\left(v\right)$ (41), ${\overline{H}}_K\left(v\right)$ (42), and ${\overline{\mathcal{H}}}_K\left(\tau \right)$ (43) and related noise-free models ${\tilde{H}}_K^M\left(v\right)$ (50), ${\tilde{H}}_K\left(v\right)$ (52), and ${\tilde{\mathcal{H}}}_K\left(\tau \right)$ (53) are estimated by the following inequalities:

$${‖{\overline{H}}_K^M\left(v\right)-{\tilde{H}}_K^M\left(v\right)‖}_2={‖{{\overline{\mathcal{H}}}_K\left(\tau \right)-\tilde{\mathcal{H}}}_K\left(\tau \right)‖}_2\le \sqrt{\alpha }\gamma {‖{\mathit{z}}_N‖}_2,$$

(54)

 where parameter

$$\gamma =\underset{1\le k\le K}{\mathit{max}} {\displaystyle \frac{\sqrt{{\sigma }_k}}{{\sigma }_k+\alpha \lambda }},$$

(55)

 and

$${‖{\overline{H}}_K\left(v\right)-{\tilde{H}}_K\left(v\right)‖}_2\le {\displaystyle \frac{\sqrt{2 {\varsigma }_1}}{\sqrt{\alpha }\left({\sigma }_K+\alpha \lambda \right)}}{‖{\mathit{z}}_N‖}_2,$$

(56)

where ${\varsigma }_1$ is the maximal singular value of ${\mathit{\Theta }}_K$ (46); ${\mathit{z}}_N={\left[\begin{array}{ccc}z\left(t_1\right)& \cdots & z\left(t_N\right)\end{array}\right]}^T$ is the vector of measurement noises.

According to inequalities (54) and (56), the accuracy of the noise-free optimal spectra approximation depends on the measurement noises, the regularization parameter, the time-scale factor, and on the singular values of the matrices ${\mathit{\Phi }}_K$ (25) and ${\mathit{\Theta }}_K$ (46), this being dependent on the number of measurements. By (54) and (56), having in mind the continuity of all the spectra considered here, we conclude that the spectra ${\overline{H}}_K^M\left(v\right)$ (41), ${\overline{H}}_K\left(v\right)$ (42), and ${\overline{\mathcal{H}}}_K\left(\tau \right)$ (43) tend to their noise-free counterparts for each $v>0$ and $\tau >0$ linearly with respect to the norm ${‖{\mathit{z}}_N‖}_2$, as ${‖{\mathit{z}}_N‖}_2\to 0$, and the faster the larger the regularization parameter $\lambda$.

3.5.3. Error of the Relaxation Modulus Model

The approximation of the material spectrum $H\left(v\right)$ by series of functions $H_K^M\left(v\right)$ (9) results in the relaxation modulus $G\left(t\right)$ approximation by the series $G_K\left(t\right)$ (11) of basis functions ${\varphi }_k\left(t\right)$ (12). Therefore, the relaxation modulus model corresponding to the relaxation spectra models (41)–(43) is described by the following equation:

$${\overline{G}}_K\left(t\right)={\sum }_{k=1}^K{\overline{g}}_k\left(\lambda \right){\varphi }_k\left(t\right).$$

(57)

The mean square error of the relaxation modulus model is as follows

$$Q_K\left({\overline{\mathit{g}}}_K\left(\lambda \right)\right)={\displaystyle \frac{1}{K}}{\sum }_{k=1}^K{\left[\overline{G}\left(t_k\right)-{\overline{G}}_K\left(t_k\right)\right]}^2.$$

(58)

In Appendix A.4, the following result is derived.

Proposition 4. 

For an arbitrary time-scale factor $\alpha$ and arbitrary regularization parameter $\lambda$, the square error of the relaxation modulus model ${\overline{G}}_K\left(t\right)$ (57) defined by (58) for the optimal model parameter ${\overline{\mathit{g}}}_K\left(\lambda \right)$ (39) is given by the following formula:

$$Q_K\left({\overline{\mathit{g}}}_K\left(\lambda \right)\right)={\displaystyle \frac{1}{K}}{\left[{\overline{\mathit{G}}}_K-{\displaystyle \frac{1}{\alpha }}{\mathit{\Phi }}_K{\overline{\mathit{g}}}_K\left(\lambda \right)\right]}^T\left[{\overline{\mathit{G}}}_K-{\displaystyle \frac{1}{\alpha }}{\mathit{\Phi }}_K{\overline{\mathit{g}}}_K\left(\lambda \right)\right]={\sum }_{k=1}^K{\displaystyle \frac{y_k^2}{{\left(\frac{{\sigma }_k}{\alpha \lambda }+1\right)}^2}},$$

(59)

where ${\sigma }_k$, $k=1,\dots ,K$, are the singular values of the matrix; ${\mathit{\Phi }}_K$ (25) and $y_k$ are the elements of the vector ${\mathit{Y}}_K$ (40).

The equality (59) yields that the accuracy of the relaxation modulus approximation depends on the following: the real relaxation modulus ${\mathit{G}}_K$ and measurement noises ${\mathit{z}}_K$ affecting the value of ${\mathit{Y}}_K={\mathit{U}}_K^T{\mathit{G}}_K+{\mathit{U}}_K^T{\mathit{z}}_K$; the scaling factor $\alpha$; the regularization parameter $\lambda$; and singular values of the matrix ${\mathit{\Phi }}_K$ (25). These, in turn, depend on the number of measurements. Note also that only the product $\alpha \lambda$, not $\alpha$ and $\lambda$ independently, affects the index $Q_K\left({\overline{\mathit{g}}}_K\left(\lambda \right)\right)$. Since the first derivative

$${\displaystyle \frac{\partial Q_K\left({\overline{\mathit{g}}}_K\left(\lambda \right)\right)}{\partial \left(\alpha \lambda \right)}}=2{\sum }_{k=1}^K{\displaystyle \frac{y_k^2{\sigma }_k\alpha \lambda }{{\left({\sigma }_k+\alpha \lambda \right)}^3}},$$

is positive for any $\alpha \lambda >0$, the error of the relaxation modulus model grows with increasing regularization parameter $\lambda$, slow for very small and very large $\alpha \lambda$.

3.6. Identification Algorithm

Allowing the above, the calculation of the relaxation spectra models involves the following steps.

• For the studied material, perform the preliminary experiment (stress relaxation test [1,33,38]) and record the measurements $\overline{G}\left(t_i\right)$, $i=1,\dots ,N$, of the relaxation modulus for pre-selected time instants (e.g., sampled with the constant period in the time interval $\left[0,T\right]$, $T<\infty$);
• Choose the time-scaling factor $\alpha$ and the number $K$ of model components comparing, for different values of $\alpha$, a few functions from the sequence $\left\{{\varphi }_k\left(t\right)\right\}$ given by (12), and creating relaxation modulus model $G_K\left(t\right)$ (11) with the experiment results $\left\{\overline{G}\left(t_i\right)\right\}$;
• Perform the experiment and record the measurements $\overline{G}\left(t_k\right)$ of the relaxation modulus at times $t_k=\alpha ·k$, $k=1,\dots ,K$;
• Compute the matrix ${\mathit{\Phi }}_K$ (25), and next, determine the SVD (35) with the singular values ${\sigma }_1,\dots ,{\sigma }_k,\dots {\sigma }_K$ of ${\mathit{\Phi }}_K$;
• Select the regularization parameter $\lambda$ such that for assumed $\alpha$ and $K$ the spectral condition number is such that

  $$\kappa \left({\mathit{\Phi }}_K+\alpha \lambda {\mathbb{I}}_{K,K}\right)={\displaystyle \frac{{\sigma }_1+\alpha \lambda }{{\sigma }_K+\alpha \lambda }}\le {10}^5;$$

  (60)
• For chosen $\lambda$, compute the regularized solution ${\overline{\mathit{g}}}_K\left(\lambda \right)$ according to (39);
• Determine the modified spectrum of relaxation frequencies ${\overline{H}}_K^M\left(v\right)$ according to (41);
• Determined the spectra of relaxation time ${\overline{\mathcal{H}}}_K\left(\tau \right)$ and frequency ${\overline{H}}_K\left(v\right)$ according to (43) and (42), respectively, as the linear combinations of the respective basis functions.

The matrix ${\mathit{\Phi }}_K$ (25) and, in particular, the singular values ${\sigma }_k$ of ${\mathit{\Phi }}_K$ depend only on the number of measurements and do not depend even on the time-scale factor $\alpha$ and on the sampling points $t_k$. Therefore, for the fixed $K$ matrix ${\mathit{\Phi }}_K$ and the SVD of ${\mathit{\Phi }}_K$—being the most space- and time-consuming task of the scheme of computational complexity $\mathcal{O}\left(N{K}^2\right)$ [39]—these must be determined only once when the identification scheme is applied for successive samples of the same material (step 4). The SVD is accessible in the form of optimized numerical procedures in most computational packets.

The condition (60) from step 5 means that the condition number does not exceed the numerically acceptable value 10⁵ [42]; data from Table 1 may be useful here.

3.7. Simulational Studies

This section presents the results of the approach with proposed numerical studies for three simulated materials whose viscoelastic properties are described by the Gauss-like and Kohlrausch–Williams–Watts (KWW) models. Gauss distributions of the relaxation spectra are examined while developing new identification methods; the best examples are as follows: ([5] (Figure 2)), ([14] (Figures 9, 11 and 17)) and ([43] (Figures 2, 3, 6–11 and 14)). The Gaussian-like distributions of the relaxation spectra were used to describe the viscoelastic properties of a lot of real materials, mainly polymers, for example, poly(methyl methacrylate) [44], polyacrylamide gels ([45] (Figure A4)), polyethylene [46] and carboxymethylcellulose (CMC) [47]. Gaussian nature has also the spectra of many biopolymers, e.g., fresh egg white-hydrocolloids [47], cold gel-like emulsions stabilized with bovine gelatin [48], xanthan gum water solution [47], some (potato, corn, wheat, and banana) native starch gels [49], and wood [24,50]. Recently, Gaussian-type relaxation spectra have also been determined for the modified asphalt binder blends ([51] (Figures 3a,c and 5a,c,e)).

The KWW model of the stretched exponential relaxation has been found by many researchers to be more appropriate than standard exponentials to describe viscoelastic processes of many materials, for example, polymer melts [52], the segmental dynamics and the glass transition behavior of poly(2-vinylpyridine) [53], relaxation of bone and bone collagen [54], and alginate films while considering glycerol concentration [55]. The KWW model, initially introduced to describe the viscoelastic relaxation processes [56,57], has also been used to model other relaxation processes occurring in materials, for example, enthalpy relaxation in Cu₄₆Zr₄₅Al₇Y₂ and Zr₅₅Cu₃₀Ni₅Al₁₀ bulk metallic glasses [58], isothermal enthalpy relaxation and density relaxation investigated for bulk Pd_42.5Cu₃₀Ni_7.5P₂₀ and Pd₄₀Ni₄₀P₂₀ metallic glasses [59], and structural relaxation of a Hf-microalloyed Co-based glassy alloy [60].

Applying the proposed identification algorithm, the best spectra models ${\overline{H}}_K^M\left(v\right)$ (41), ${\overline{H}}_K\left(v\right)$ (42), and ${\overline{\mathcal{H}}}_K\left(\tau \right)$ (43) were determined for a few numbers of measurements. The smoothness of the models is estimated by their integral square norms described in Proposition 1. For the spectra ${\overline{H}}_K^M\left(v\right)$ and ${\overline{\mathcal{H}}}_K\left(\tau \right)$, these norms are equal and uniquely characterized by the middle equality in (44), which yields the following:

$${‖{\overline{H}}_K^M\left(v\right)‖}_2={‖{\overline{\mathcal{H}}}_K\left(\tau \right)‖}_2={\displaystyle \frac{1}{\sqrt{\alpha }}}\sqrt{{\overline{\mathit{g}}}_K^T\left(\lambda \right){\mathit{\Phi }}_K{\overline{\mathit{g}}}_K\left(\lambda \right)}.$$

(61)

By (23), for arbitrary relaxation spectra $\mathcal{H}\left(\tau \right)$ and $H^M\left(v\right)$ of the real material, the analogous equality of the norms holds, i.e., ${‖H^M\left(v\right)‖}_2={‖\mathcal{H}\left(\tau \right)‖}_2$. The formula describing the square norm of ${\overline{H}}_K\left(v\right)$ (42) directly results from the equality in (45), from which the following can be obtained:

$${‖{\overline{H}}_K\left(v\right)‖}_2={\displaystyle \frac{\sqrt{2}}{\alpha \sqrt{\alpha }}}\sqrt{{\overline{\mathit{g}}}_K^T\left(\lambda \right){\mathit{\Theta }}_K{\overline{\mathit{g}}}_K\left(\lambda \right)},$$

(62)

with the matrix ${\mathit{\Theta }}_K$ defined by (46).

The errors of the relaxation modulus models are estimated using index $Q_K\left({\overline{\mathit{g}}}_K\left(\lambda \right)\right)$ (58) and expressed for the optimal models ${\overline{H}}_K^M\left(v\right)$ (41), ${\overline{H}}_K\left(v\right)$ (42), and ${\overline{\mathcal{H}}}_K\left(\tau \right)$ (43) by Formula (59) from Proposition 4.

The errors of the relaxation spectra models are measured directly by the integral $\mathcal{J}\left({\mathit{g}}_K\right)$, introduced by Equation (16) for model $H_K^M\left(v\right)$ and by Equation (20) for the relaxation time spectrum model. By (19), having in mind the notation (25), this index for the optimal model ${\overline{H}}_K^M\left(v\right)$ (41) can be expressed as follows (compare to (26)):

$$\mathcal{J}\left({\overline{\mathit{g}}}_K\left(\lambda \right)\right)={\int }_0^{\mathrm{\infty }}{H^M\left(v\right)}^{2}dv-2{\mathit{G}}_K^T{\overline{\mathit{g}}}_K\left(\lambda \right)+{\displaystyle \frac{1}{\alpha }}{\overline{\mathit{g}}}_K^T\left(\lambda \right){\mathit{\Phi }}_K{\overline{\mathit{g}}}_K\left(\lambda \right),$$

(63)

where ${\mathit{G}}_K$ is the vector of noise-free values of relaxation modulus defined below Equation (51). The error $\mathcal{J}\left({\overline{\mathit{g}}}_K\left(\lambda \right)\right)$ for the model ${\overline{\mathcal{H}}}_K\left(\tau \right)$ (43) is obviously identical due to the identity of indices (16) and (20) proved above.

The error $\mathcal{J}\left({\overline{\mathit{g}}}_K\left(\lambda \right)\right)$ related to square of the norm of real spectra ${‖H^M\left(v\right)‖}_2^2={‖\mathcal{H}\left(\tau \right)‖}_2^2$ is measured by the relative index

$${\mathcal{J}}_{rel}\left({\overline{\mathit{g}}}_K\left(\lambda \right)\right)={\displaystyle \frac{\mathcal{J}\left({\overline{\mathit{g}}}_K\left(\lambda \right)\right)}{{‖\mathcal{H}\left(\tau \right)‖}_2^2}}.$$

(64)

The “real” materials and the optimal models were simulated in Matlab R2023b, The Mathworks, Inc., Natick, MA, USA. For the singular value decomposition procedure, svd was applied.

3.8. Identification of Uni-Mode Gauss-like Spectrum

Consider material whose rheological properties are characterized by the uni-modal Gauss-like distribution [10,61]:

$$\mathcal{H}\left(\tau \right)=\vartheta e^{-{\left({\displaystyle \frac{1}{\tau }}-m\right)}^2/q}/\tau ,$$

(65)

where the parameters are as follows [10,27,61]: $\vartheta =31.52 \mathrm{k}\mathrm{P}\mathrm{a}·\mathrm{s}$, $m=0.0912 {\mathrm{s}}^{-1}$, and $q=3.25\times {10}^{-3} {\mathrm{s}}^{-2}$. The relaxation modulus is displayed below [10]:

$$G\left(t\right)={\displaystyle \frac{\sqrt{\pi q}}{2}}\vartheta e^{{\displaystyle \frac{1}{4}}{t}^{2}q-mt}erfc\left({\displaystyle \frac{\frac{1}{2}tq-m}{\sqrt{q}}}\right),$$

(66)

where the complementary error function $erfc\left(x\right)$ is defined as follows ([62] Equation (8.250.4)):

$$erfc\left(x\right)={\displaystyle \frac{2}{\sqrt{\pi } }}{\int }_x^{\infty }{e}^{{-z}^2}dz.$$

(67)

By the first equality in (3), the following spectrum of relaxation frequencies corresponds to (65)

$$H\left(v\right)=\vartheta v{e}^{-{\left(v-m\right)}^2/q},$$

(68)

whence, in view of (4), the modified spectrum is described as outlined below:

$$H^M\left(v\right)=\vartheta e^{-{\left(v-m\right)}^2/q}.$$

(69)

In Appendix A.5, the analytical Formulas (A10) and (A13) are derived describing the norms of the spectra $\mathcal{H}\left(\tau \right)$ (65), $H\left(v\right)$ (68), and $H^M\left(v\right)$ (69). These norms are as follows: ${{‖H^M\left(v\right)‖}_2=‖\mathcal{H}\left(\tau \right)‖}_2=8.422432\mathrm{k}\mathrm{P}\mathrm{a}·{\mathrm{s}}^{1/2}$ and ${‖H\left(v\right)‖}_2=0.805043\mathrm{k}\mathrm{P}\mathrm{a}·{\mathrm{s}}^{-1/2}$.

The preliminary relaxation test experiment was performed (step 1) and the measurements of the relaxation modulus $G\left(t\right)$ (66) were recorded for 200 s, selected following [10,27,61]. Then, the time scale factors $\alpha$ have been selected by comparison of the courses of experiment results $\left\{\overline{G}\left(t_i\right)\right\}$ and basis functions ${\varphi }_k\left(t\right)$ (12) for a few $k$. Next, to simulate the experiment, $K$ sampling instants $t_k=\alpha k$ were generated with the constant period $\alpha$ for $K=\mathrm{20,50,100,150,200}$ measurements. Additive measurement noises $z\left(t_k\right)$ were selected independently by random choice with uniform distribution on the interval $\left[-10, 10\right] Pa$, i.e., double stronger than noises assumed in the previous papers [10,61]. The measurements $\overline{G}\left(t_k\right)$ were recorded. For successive $K$, the matrix ${\mathit{\Phi }}_K$ (25) and SVD (35) were determined. Next, the regularization parameters $\lambda$ were selected according to the spectral condition number rule (60); their values are given in

Table 3. For the uni-mode Gauss-like spectrum, $\mathcal{H}\left(\tau \right)$ (65), and the models ${\overline{H}}_K^M\left(v\right)$ (41),${\overline{H}}_K\left(v\right)$ (42), and ${\overline{\mathcal{H}}}_K\left(\tau \right)$ (43): time-scale factors $\alpha$; numbers of model summands $K$; regularization parameters $\lambda$; the model’s smoothness indices ${‖{\overline{H}}_K^M\left(v\right)‖}_2={‖{\overline{\mathcal{H}}}_K\left(\tau \right)‖}_2$ (61) and ${‖{\overline{H}}_K\left(v\right)‖}_2$ (62); the mean square relaxation modulus approximation index $Q_K\left({\overline{\mathit{g}}}_K\left(\lambda \right)\right)$ (58); norms ${‖{\overline{\mathit{g}}}_K\left(\lambda \right)‖}_2$ (49) of the model parameter vectors; the integral square approximation indices $\mathcal{J}\left({\overline{\mathit{g}}}_K\left(\lambda \right)\right)$ (63); and relative index ${\mathcal{J}}_{rel}\left({\overline{\mathit{g}}}_K\left(\lambda \right)\right)$ (64).

$\mathit{K}$	$\mathit{\alpha }$$\left[\mathbf{s}\right]$	$\mathit{\lambda }$$\left[{\mathit{s}}^{-1}\right]$	${‖{\overline{\mathcal{H}}}_{\mathit{K}}\left(\mathit{\tau }\right)‖}_2$$[\mathbf{k}\mathbf{P}\mathbf{a}·{\mathbf{s}}^{1/2}]$	${‖{\overline{\mathit{H}}}_{\mathit{K}}\left(\mathit{v}\right)‖}_2$$\left[\mathbf{k}\mathbf{P}\mathbf{a}·{\mathbf{s}}^{-1/2}\right]$	${\mathit{Q}}_{\mathit{K}}\left({\overline{\mathit{g}}}_{\mathit{K}}\left(\mathit{\lambda }\right)\right)$$\left[{\mathbf{k}\mathbf{P}\mathbf{a}}^2\right]$	${‖{\overline{\mathit{g}}}_{\mathit{K}}\left(\mathit{\lambda }\right)‖}_2$ $\left[\mathbf{k}\mathbf{P}\mathbf{a}·\mathbf{s}\right]$	$\mathcal{J}\left({\overline{\mathit{g}}}_{\mathit{K}}\left(\mathit{\lambda }\right)\right)$$\left[{\mathbf{k}\mathbf{P}\mathbf{a}}^2·\mathbf{s}\right]$	${\mathcal{J}}_{\mathit{r}\mathit{e}\mathit{l}}\left({\overline{\mathit{g}}}_{\mathit{K}}\left(\mathit{\lambda }\right)\right)$
20	4.90	3.1 × 10−6	8.194688	1.000985	3.681218 × 10−5	8.752841 × 103	3.548595	0.0500
50	3.75	9.0 × 10−6	8.167268	0.847115	3.486716 × 10−5	4.639282 × 103	2.461727	0.0347
100	3.98	7.8 × 10−6	8.193556	0.861458	3.272773 × 10−5	7.334379 × 103	2.439776	0.0344
150	3.75	5.5 × 10−6	8.248417	0.865259	3.310334 × 10−5	1.281206 × 104	2.625005	0.0370
200	4.08	9.5 × 10−6	8.173553	0.869989	3.299329 × 10−5	8.550751 × 103	2.875342	0.0405

. The optimal model parameters ${\overline{\mathit{g}}}_K\left(\lambda \right)$ (39) and the models ${\overline{H}}_K^M\left(v\right)$ (41), ${\overline{H}}_K\left(v\right)$ (42), and ${\overline{\mathcal{H}}}_K\left(\tau \right)$ (43) were determined. The best models are depicted in Figure 5, Figure 6 and Figure 7 together with the real spectra (65), (68), and (69) marked by red lines. Small subfigures show fitting near the maximum of the real spectrum. The respective relaxation modulus models ${\overline{G}}_K\left(t\right)$ (57) are plotted in Figure 8 for $K=20$ and 200, where the measurements $\overline{G}\left(t_k\right)$ of the real modulus $G\left(t\right)$ (66) are also marked; the small subfigures confirm the excellent model fit. For $K=200$, logarithmic time scale is used. In Table 3, the norms ${‖{\overline{H}}_K^M\left(v\right)‖}_2={‖{\overline{\mathcal{H}}}_K\left(\tau \right)‖}_2$ (61), ${‖{\overline{H}}_K\left(v\right)‖}_2$ (62), and the norms ${‖{\overline{\mathit{g}}}_K\left(\lambda \right)‖}_2$—expressed in (49)—of the optimal model parameters are given. In addition, the integral square approximation index $\mathcal{J}\left({\overline{\mathit{g}}}_K\left(\lambda \right)\right)$ (63) together with the relative index ${\mathcal{J}}_{rel}\left({\overline{\mathit{g}}}_K\left(\lambda \right)\right)$ (64) and the mean square approximation index $Q_K\left({\overline{\mathit{g}}}_K\left(\lambda \right)\right)$ (58) are provided in Table 3.

The relative spectrum approximation index ${\mathcal{J}}_{rel}\left({\overline{\mathit{g}}}_K\left(\lambda \right)\right)$ (64) does not exceed 5%; additionally, the values of the relaxation modulus approximation index $Q_K\left({\overline{\mathit{g}}}_K\left(\lambda \right)\right)$ (58) indicate the excellent model fit; models ${\overline{G}}_K\left(t\right)$ (57) practically coincide with the measurement points $\overline{G}\left(t_k\right)$, see Figure 8. An inspection of Figure 5, Figure 6 and Figure 7 shows that for the number of $K\ge 50$ measurements, satisfactory approximation of the relaxation spectra was obtained while maintaining the consistency of the maxima of real spectra and their models.

3.9. Identification of Double-Mode Gauss-like Spectrum

Consider now the viscoelastic material of the relaxation spectrum described by the double-mode Gauss-like distribution considered in [10,27,28,46]:

$$\mathcal{H}\left(\tau \right)=\left[{\vartheta }_1{e}^{-{\left({\displaystyle \frac{1}{\tau }}-m_1\right)}^2/q_1}+{\vartheta }_2{e}^{-{\left({\displaystyle \frac{1}{\tau }}-m_2\right)}^2/q_2}\right]/\tau ,$$

(70)

where the parameters are as follows [27,28]: ${\vartheta }_1=467 \mathrm{P}\mathrm{a}·\mathrm{s}$, $m_1=0.0037 {\mathrm{s}}^{-1}$, $q_1=1.124261\times {10}^{-6} {\mathrm{s}}^{-2}$, ${\vartheta }_2=39 \mathrm{P}\mathrm{a}·\mathrm{s}$, $m_2=0.045 {\mathrm{s}}^{-1}$, and $q_2=1.173\times {10}^{-3} {\mathrm{s}}^{-2}$. Therefore, the corresponding spectrum of relaxation frequencies is as follows

$$H\left(v\right)={\vartheta }_{1}v{e}^{-{\left(v-m_1\right)}^2/q_1}+{\vartheta }_{2}v{e}^{-{\left(v-m_2\right)}^2/q_2},$$

(71)

and, in view of (4), the modified spectrum is described by $H^M\left(v\right)=H\left(v\right)/v$. By (66), the related real relaxation modulus is outlined below:

$$G\left(t\right)={\displaystyle \frac{\sqrt{\pi }}{2}}\left[{\vartheta }_1\sqrt{q_1} e^{{\displaystyle \frac{1}{4}}{t}^2{q}_1-m_{1}t}erfc\left({\displaystyle \frac{\frac{1}{2}t{q}_1-m_1}{\sqrt{q_1}}}\right)+{\vartheta }_2\sqrt{q_2} e^{{\displaystyle \frac{1}{4}}{t}^2{q}_2-m_{2}t}erfc\left({\displaystyle \frac{\frac{1}{2}t{q}_2-m_2}{\sqrt{q_2}}}\right)\right].$$

(72)

In Appendix A.6, the analytical Formulas (A17) and (A19) are derived to describe the square norms of the “real” spectra $H\left(v\right)$ (71), $\mathcal{H}\left(\tau \right)$ (70), and $H^M\left(v\right)$, which are as follows: ${{‖H^M\left(v\right)‖}_2=‖\mathcal{H}\left(\tau \right)‖}_2=19.257051 \mathrm{P}\mathrm{a}·{\mathrm{s}}^{1/2}$ and ${‖H\left(v\right)‖}_2=0.394490\mathrm{P}\mathrm{a}·{\mathrm{s}}^{-1/2}$.

Based on the course of the modulus $G\left(t\right)$ (72), in the preliminary experiment, $N=5000$ sampling instants were generated with the constant period in the time interval $\mathcal{T}=\left[0,1550\right]$ s, c.f., [27,28]. Following [27,28], additive measurement noises $z\left(t_i\right)$ were selected independently by random choice with uniform distribution on the interval $\left[-0.005, 0.005\right] \mathrm{P}\mathrm{a}$. The same as before, several $K$ relaxation spectra models were determined using the proposed identification algorithm. The values of selected regularization parameters $\lambda$, the norms ${‖{\overline{H}}_K^M\left(v\right)‖}_2={‖{\overline{\mathcal{H}}}_K\left(\tau \right)‖}_2$, ${‖{\overline{H}}_K\left(v\right)‖}_2$, and ${‖{\overline{\mathit{g}}}_K\left(\lambda \right)‖}_2$ and the indices $\mathcal{J}\left({\overline{\mathit{g}}}_K\left(\lambda \right)\right)$, ${\mathcal{J}}_{rel}\left({\overline{\mathit{g}}}_K\left(\lambda \right)\right)$, and $Q_K\left({\overline{\mathit{g}}}_K\left(\lambda \right)\right)$ (58) are presented in

Table 4. For the double-mode Gauss-like spectrum $\mathcal{H}\left(\tau \right)$ (70) and the models ${\overline{H}}_K^M\left(v\right)$ (41),${\overline{H}}_K\left(v\right)$ (42) and ${\overline{\mathcal{H}}}_K\left(\tau \right)$ (43): time-scale factors $\alpha$; numbers of model summands $K$; regularization parameters $\lambda$; the model’s smoothness indices ${‖{\overline{H}}_K^M\left(v\right)‖}_2={‖{\overline{\mathcal{H}}}_K\left(\tau \right)‖}_2$ and ${‖{\overline{H}}_K\left(v\right)‖}_2$; the mean square relaxation modulus approximation index $Q_K\left({\overline{\mathit{g}}}_K\left(\lambda \right)\right)$ (58); norms ${‖{\overline{\mathit{g}}}_K\left(\lambda \right)‖}_2$ (49) of the model parameter vectors; the integral square approximation indices $\mathcal{J}\left({\overline{\mathit{g}}}_K\left(\lambda \right)\right)$ (63); and relative index ${\mathcal{J}}_{rel}\left({\overline{\mathit{g}}}_K\left(\lambda \right)\right)$ (64).

$\mathit{K}$	$\mathit{\alpha }$$\left[\mathit{s}\right]$	$\mathit{\lambda }$$\left[{\mathit{s}}^{-1}\right]$	${‖{\overline{\mathcal{H}}}_{\mathit{K}}\left(\mathit{\tau }\right)‖}_2$$[\mathbf{P}\mathbf{a}·{\mathbf{s}}^{1/2}]$	${‖{\overline{\mathit{H}}}_{\mathit{K}}\left(\mathit{v}\right)‖}_2$$\left[\mathbf{P}\mathbf{a}·{\mathbf{s}}^{-1/2}\right]$	${\mathit{Q}}_{\mathit{K}}\left({\overline{\mathit{g}}}_{\mathit{K}}\left(\mathit{\lambda }\right)\right)$$\left[{\mathbf{P}\mathbf{a}}^2\right]$	${‖{\overline{\mathit{g}}}_{\mathit{K}}\left(\mathit{\lambda }\right)‖}_2$$\left[\mathbf{P}\mathbf{a}·\mathbf{s}\right]$	$\mathcal{J}\left({\overline{\mathit{g}}}_{\mathit{K}}\left(\mathit{\lambda }\right)\right)$$\left[{\mathbf{P}\mathbf{a}}^2·\mathbf{s}\right]$	${\mathcal{J}}_{\mathit{r}\mathit{e}\mathit{l}}\left({\overline{\mathit{g}}}_{\mathit{K}}\left(\mathit{\lambda }\right)\right)$
50	22.5	8 × 10−7	17.008811	0.377417	8.91256 × 10−6	2.63874 × 104	83.212916	0.224394
100	16.3	1.2 × 10−6	16.749729	0.415414	8.20297 × 10−6	2.38674 × 104	89.54064	0.241457
150	9.35	2.1 × 10−6	16.145506	0.440330	8.33036 × 10−6	1.68329 × 104	90.59845	0.244309
200	6.5	3.1 × 10−6	16.447062	0.404095	8.38748 × 10−6	1.321209 × 104	85.71129	0.2311311
300	5.2	4 × 10−6	16.650568	0.401779	8.04576 × 10−6	1.228249 × 104	87.92225	0.237093

. The optimal models ${\overline{H}}_K^M\left(v\right)$ (41), ${\overline{H}}_K\left(v\right)$ (42), and ${\overline{\mathcal{H}}}_K\left(\tau \right)$ (43) are depicted in Figure 9, Figure 10 and Figure 11 together with the real spectra plotted by red lines. The respective relaxation modulus models ${\overline{G}}_K\left(t\right)$ (57) are plotted in Figure 12 with the modulus $G\left(t\right)$ (72) measurements.

However, for double-mode spectrum, the relative spectrum approximation index ${\mathcal{J}}_{rel}\left({\overline{\mathit{g}}}_K\left(\lambda \right)\right)$ is as much as 25%, an inspection of Figure 9, Figure 10 and Figure 11 indicates a satisfactory approximation of the real spectra while maintaining the locations of both their maxima. Excellent models ${\overline{G}}_K\left(t\right)$ fit is confirmed by the values $Q_K\left({\overline{\mathit{g}}}_K\left(\lambda \right)\right)$ and Figure 12.

3.10. Identification of KWW Relaxation Spectrum

The relaxation spectrum of the KWW model of the stretched exponential relaxation is described by the following [56]:

$$G\left(t\right)=G_0{e}^{-{\left({\displaystyle \frac{t}{{\tau }_r}}\right)}^{\beta }},$$

(73)

where the stretching exponent $0<\beta <1$, ${\tau }_r$ is the relaxation time, and $G_0$ is the initial shear modulus, which has a unimodal [57] relaxation spectrum described by the infinite series [56,57]:

$$\mathcal{H}\left(\tau \right)={\displaystyle \frac{G_0}{\pi }} {\sum }_{k=1}^{\infty }{\displaystyle \frac{{\left(-1\right)}^{k+1}}{k!}}\mathit{sin}\left(\pi \beta k\right) \mathsf{\Gamma }\left(\beta k+1\right) {\left({\displaystyle \frac{\tau }{{\tau }_r}}\right)}^{\beta k},$$

(74)

where $\mathsf{\Gamma }\left(n\right)$ is Euler’s gamma function ([62] Equation (8.310.1)). However, for the stretching exponent $\beta =0.5$, spectrum $\mathcal{H}\left(\tau \right)$ has simple analytical form [57]:

$$\mathcal{H}\left(\tau \right)={\displaystyle \frac{G_0}{2\sqrt{\pi }}}\sqrt{{\displaystyle \frac{\tau }{{\tau }_r}}} e^{-{\displaystyle \frac{\tau }{4{\tau }_r }}}.$$

(75)

The stretching exponent 0.5 is assumed, for which the relaxation spectrum is given by the analytical formula because the effectiveness of the identification method can only be verified when the assumed spectrum is exactly known. The exponent $\beta =0.5$ has been reported by Plazek and Ngai [63] for poly(methylphenylsiloxane) at the glass temperature $T_g=200 \mathrm{K}$; however, the related relaxation time ${\tau }_r$ is not reported in [63]. In 1993, Böhmer et al. [64], based on the literature and also private communications, have presented data concerning stretched exponential relaxation from about 70 amorphous polymeric glass formers (supercooled liquids and disordered crystals). The exponent $\beta =0.5$ has been experimentally obtained for sorbitol, dehydroabietic acid, BBKDE, 1,4-cis-polyisoprene, and silicate flint glass ([64] (Table I)). The coefficients $\beta$ near 0.5 have been found for toluene (25%) ($\beta =0.52$), BCDE ($\beta =0.51$), polyisobutylene ($\beta =0.55$), and several other forms of glass. However, this paper also does not contain data on the related relaxation times. Recently, Chen et al. [65] applied the KWW model to describe the viscoelastic properties of the cross-linked polystyrene estimating the following parameters of the model (73) at temperature $110 °\mathrm{C}$: $G_0=0.78 \mathrm{M}\mathrm{P}\mathrm{a}$, $\beta =0.59$, and ${\tau }_r=1.08 \mathrm{s}$ ([65] (Table 4)). For the purpose of numerical tests of the algorithm, $\beta =0.59$ was replaced here by $\beta =0.5$. The discrepancy between the relaxation modulus $G\left(t\right)$ (73) for $\beta =0.5$ and $\beta =0.59$ is illustrated in Figure 13; the mean least-squares error for 500 equidistant sampling points between these modulus is equal to 1.2853 × 10⁻⁴ MPa².

The relaxation frequency spectrum corresponding to (75) is as follows

$$H\left(v\right)={\displaystyle \frac{G_0}{2\sqrt{\pi }}}\sqrt{{\displaystyle \frac{1}{{\tau }_r·v}}} e^{-{\displaystyle \frac{1}{4 {\tau }_r·v}}},$$

(76)

while the modified spectrum is described by the following:

$$H^M\left(v\right)={\displaystyle \frac{G_0}{2\sqrt{\pi }}}\sqrt{{\displaystyle \frac{1}{{\tau }_r·v^3}}} e^{-{\displaystyle \frac{1}{4 {\tau }_r·v}}}.$$

(77)

In the preliminary experiment, $N=500$ sampling instants were generated with the constant period in the time interval $\mathcal{T}=\left[0,50\right]$ s, selected based on the course of the modulus $G\left(t\right)$ (73) in Figure 13. Additive measurement noises $z\left(t_i\right)$ were selected independently by random choice with uniform distribution on the interval $\left[-0.5, 0.5\right] k\mathrm{P}\mathrm{a}$. The values of the time scale factors $\alpha$ selected for a few values of $K$ and the regularization parameters $\lambda$ selected according to the rule (60) are given in

Table 5. For the KWW spectrum $\mathcal{H}\left(\tau \right)$ (75) and the models ${\overline{H}}_K^M\left(v\right)$ (41), ${\overline{H}}_K\left(v\right)$ (42), and ${\overline{\mathcal{H}}}_K\left(\tau \right)$ (43): time-scale factors $\alpha$; numbers of model summands $K$; regularization parameters $\lambda$; the model’s smoothness indices ${‖{\overline{H}}_K^M\left(v\right)‖}_2={‖{\overline{\mathcal{H}}}_K\left(\tau \right)‖}_2$ and ${‖{\overline{H}}_K\left(v\right)‖}_2$; the mean square relaxation modulus approximation index $Q_K\left({\overline{\mathit{g}}}_{\mathit{K}}\left(\mathit{\lambda }\right)\right)$ (58); norms ${‖{\overline{\mathit{g}}}_K\left(\lambda \right)‖}_2$ (49) of the model parameter vectors; the integral square approximation indices $\mathcal{J}\left({\overline{\mathit{g}}}_K\left(\lambda \right)\right)$ (63); and relative index ${\mathcal{J}}_{rel}\left({\overline{\mathit{g}}}_K\left(\lambda \right)\right)$ (64).

$\mathit{K}$	$\mathit{\alpha }\left[\mathit{s}\right]$	$\mathit{\lambda }\left[{\mathit{s}}^{-1}\right]$	${‖{\overline{\mathcal{H}}}_{\mathit{K}}\left(\mathit{\tau }\right)‖}_2[\mathbf{M}\mathbf{P}\mathbf{a}·{\mathbf{s}}^{1/2}]$	${‖{\overline{\mathit{H}}}_{\mathit{K}}\left(\mathit{v}\right)‖}_2\left[\mathbf{M}\mathbf{P}\mathbf{a}·{\mathbf{s}}^{-1/2}\right]$	${\mathit{Q}}_{\mathit{K}}\left({\overline{\mathit{g}}}_{\mathit{K}}\left(\mathit{\lambda }\right)\right)\left[\mathbf{M}{\mathbf{P}\mathbf{a}}^2\right]$	${‖{\overline{\mathit{g}}}_{\mathit{K}}\left(\mathit{\lambda }\right)‖}_2\left[\mathbf{M}\mathbf{P}\mathbf{a}·\mathbf{s}\right]$	$\mathcal{J}\left({\overline{\mathit{g}}}_{\mathit{K}}\left(\mathit{\lambda }\right)\right)\left[{\mathbf{M}\mathbf{P}\mathbf{a}}^2·\mathbf{s}\right]$	${\mathcal{J}}_{\mathit{r}\mathit{e}\mathit{l}}\left({\overline{\mathit{g}}}_{\mathit{K}}\left(\mathit{\lambda }\right)\right)$
25	0.8	2 × 10−5	0.454723	0.267896	9.46218 × 10−8	76.901664	1.83447 × 10−3	8.77098 × 10−3
50	0.65	7 × 10−5	0.456956	0.284479	8.70335 × 10−8	29.800961	9.82909 × 10−4	4.69948 × 10−3
75	0.6	7.5 × 10−5	0.4571396	0.289756	8.52022 × 10−8	33.705039	8.50259 × 10−4	4.06526 × 10−3
100	0.65	8.5 × 10−5	0.456765	0.284243	8.18454 × 10−8	33.6572199	8.26024 × 10−4	3.94938 × 10−3
150	0.6	1 × 10−4	0.457356	0.289670	8.28962 × 10−8	35.262493	8.41199 × 10−4	4.02194 × 10−3
200	0.6	1.5 × 10−4	0.456657	0.289207	8.25479 × 10−8	27.087999	7.34565 × 10−4	3.51210 × 10−3
300	0.55	1.6 × 10−4	0.456875	0.294931	7.97708 × 10−8	30.574734	7.12159 × 10−4	3.40497 × 10−3
400	0.55	1.6 × 10−4	0.456791	0.294901	8.10268 × 10−8	35.581507	7.12571 × 10−4	3.40694 × 10−3

. The best models ${\overline{H}}_K^M\left(v\right)$ (41), ${\overline{H}}_K\left(v\right)$ (42), and ${\overline{\mathcal{H}}}_K\left(\tau \right)$ (43) are depicted in Figure 14, Figure 15 and Figure 16 together with the real spectra (red lines). The respective models ${\overline{G}}_K\left(t\right)$ (57) are plotted in Figure 17, with the modulus $G\left(t\right)$ (73) measurements. In Table 5, the norms ${‖{\overline{H}}_K^M\left(v\right)‖}_2={‖{\overline{\mathcal{H}}}_K\left(\tau \right)‖}_2$, ${‖{\overline{H}}_K\left(v\right)‖}_2$, and the norms ${‖{\overline{\mathit{g}}}_K\left(\lambda \right)‖}_2$ of the optimal parameters ${\overline{\mathit{g}}}_K\left(\lambda \right)$ (39) are given. Since from (75) and (A3) we have the following:

$${‖\mathcal{H}\left(\tau \right)‖}_2^2={\displaystyle \frac{G_0^2}{4\pi {\tau }_r}}{\int }_0^{\mathrm{\infty }}\tau e^{-{\displaystyle \frac{\tau }{2{\tau }_r }}}d\tau ={\displaystyle \frac{G_0^2{\tau }_r}{\pi }},$$

the norm ${‖\mathcal{H}\left(\tau \right)‖}_2={‖H^M\left(v\right)‖}_2=G_0\sqrt{{\tau }_r/\pi }=0.457331 \mathrm{M}\mathrm{P}\mathrm{a}·{\mathrm{s}}^{1/2}$. In turn, the norm of the relaxation frequency spectrum (76) is infinite; however, in Table 5, the model’s norms ${‖{\overline{H}}_K\left(v\right)‖}_2$ are given. The integral indices $\mathcal{J}\left({\overline{\mathit{g}}}_K\left(\lambda \right)\right)$ (63), the relative index ${\mathcal{J}}_{rel}\left({\overline{\mathit{g}}}_K\left(\lambda \right)\right)$, (64) and the mean square relaxation modulus approximation index $Q_K\left({\overline{\mathit{g}}}_K\left(\lambda \right)\right)$ (58) are also presented in Table 5.

The relative integral square index of the spectra approximation ${\mathcal{J}}_{rel}\left({\overline{\mathit{g}}}_K\left(\lambda \right)\right)$ does not exceed 0.5% for $K\ge 50$ measurements, which means a better approximation of the assumed relaxation spectrum in the whole range of time/frequency relaxation variation, i.e., from zero to infinity, than in the case of Gaussian spectra. Also related to the relaxation modulus index, $Q_K\left({\overline{\mathit{g}}}_K\left(\lambda \right)\right)$, not exceeding 10⁻⁷, confirms the perfect approximation of the relaxation modulus measurements. In the case of this unimodal spectrum, increasing the number of measurements, i.e., the components of the series that create the models, does not significantly affect the quality of these models, which, in addition to the indices in Table 5, is also confirmed by a review of Figure 14, Figure 15 and Figure 16. For $K\ge 100$, the courses of the spectra models for increasing $K$ remain practically almost identical, although a slight improvement in the fit to the real spectra can be seen in the values of the indices $\mathcal{J}\left({\overline{\mathit{g}}}_K\left(\lambda \right)\right)$ and ${\mathcal{J}}_{rel}\left({\overline{\mathit{g}}}_K\left(\lambda \right)\right)$. The relative index of the spectrum $\mathcal{H}\left(\tau \right)$ approximation for $K\ge 100$ falls below 0.41%.

3.11. Applicability of the Approach for Identification of Relaxation Spectra of Different Types

The rough condition of the approach’s successful applicability follows from the boundary properties of the optimal models ${\overline{\mathcal{H}}}_K\left(\tau \right)$ (43) and ${\overline{H}}_K\left(v\right)$ (42), yielded by the properties of the basis functions ${\mathcal{h}}_k\left(\tau \right)$ (15) and $h_k\left(v\right)v$, where $h_k\left(v\right)$ is given by (7). Since for $\tau \to 0^{+}$ and $\tau \to \mathrm{\infty }$, the basis functions ${\mathcal{h}}_k\left(\tau \right)\to 0$, the best model ${\overline{\mathcal{H}}}_K\left(\tau \right)$ (43) also tends to zero as the relaxation time $\tau$ tends to zero and to infinity, which limit the scope of applicability of this model to real relaxation time spectra that satisfy zero boundary conditions. For the relaxation frequencies $v=0$ and $v\to \mathrm{\infty }$, the basis functions $h_k\left(v\right)v\to 0$, that is the basis functions $h_k\left(v\right)v$ of the relaxation frequency model ${\overline{H}}_K\left(v\right)$ (42) also have zero boundary conditions. Therefore, in terms of the relaxation frequency, the scope of applicability of the model and method to real relaxation frequency spectra is confined to the spectra of zero boundary conditions, too.

The real relaxation time and frequency spectra and the known spectra models tend to zero as the relaxation time $\tau$ and the relaxation frequency $v$ tend to infinity. Therefore, the properties of the spectra for $\tau \to 0^{+}$ and $v\to 0^{+}$ are essential here.

The examples presented above showed that the approach proposed can be applied for Gauss-like relaxation spectra, both uni- and double mode, and for the KWW spectrum of the stretching exponent $\beta =0.5$. However, it is easy to check that for the relaxation spectrum $\mathcal{H}\left(\tau \right)$ (74) both zero boundary conditions are satisfied. Therefore, the proposed identification method can be applied to determine the spectrum of materials whose relaxation processes have KWW stretched exponential nature. This is also important that the optimal model ${\overline{\mathcal{H}}}_K\left(\tau \right)$ (43), given by a finite series, may prove to be more useful than the original KWW infinite series spectrum (74) for many applications.

A multiplicative model that combines the power law with the stretched exponential relaxation described by Equation (8) in [66]:

$$\mathcal{H}\left(\tau \right)=n_{\alpha }{G}_c{\left({\displaystyle \frac{\tau }{{\tau }_{\alpha }}}\right)}^{n_{\alpha }} e^{-{\left({\displaystyle \frac{\tau }{{\tau }_{\alpha }}}\right)}^{\beta }},$$

(78)

where ${\tau }_{\alpha }$ is the longest relaxation time, $G_c$ is the plateau modulus, the stretching parameter $0<\beta \le 1$, and the exponent $0<n_{\alpha }<1$, was applied for modeling spectrum of bitumen in the vicinity of the glass transition [66]. The unimodal spectrum (78), named by the authors as the broadened power-law spectrum model [66], satisfies both zero boundary conditions—compare to ([66] (Figure 11a))—and therefore, is also within the scope of the proposed algorithm’s applicability.

However, the well-known Baumgaertel, Schausberger, and Winter (BSW) spectrum [15,67] used to describe the viscoelasticity of polybutadiene (PBD) [68], polydisperse polymer melts [8], polymethylmethacrylate (PMMA) [68], and many other materials, is described by the following model:

$$\mathcal{H}\left(\tau \right)=\left\{{\beta }_1{\left({\displaystyle \frac{\tau }{{\tau }_c}}\right)}^{{\rho }_1}+{\beta }_2{\left({\displaystyle \frac{\tau }{{\tau }_c}}\right)}^{{\rho }_2}\right\}{e}^{-{\displaystyle \frac{\tau }{{\tau }_{max}}}},$$

with positive coefficients ${\beta }_1$, ${\beta }_2$ and relaxation times ${\tau }_c$, ${\tau }_{max}$, which tends to infinity for $\tau \to 0$ whenever at least one of the parameters ${\rho }_1$ and ${\rho }_2$ is negative. Therefore, this is the case for real material models, compare [8,15,67,68], when the optimal model ${\overline{\mathcal{H}}}_K\left(\tau \right)$ (43) cannot well-approximate this spectrum.

Likewise, the real relaxation spectra modeled by pure inverse power laws [69], for example, a combined four-interval power model with fractional exponents describing a solution-polymerized styrene butadiene rubber [70] or a power type spectrum with an exponent of −1/2 describing the cross-linking polymers at their gel point [71], cannot be successfully identified by the proposed approach. The relaxation time spectra of the fractional Maxwell model and the elementary fractional Scott–Blair model also lose the zero boundary condition at zero relaxation time, see [61] (Proposition 2, Equation (19)).

3.12. Direct Identification of the Relaxation Spectra of Viscoelastic Solid Materials

For isotropic viscoelastic solids [31]

$$\underset{t\to \mathrm{\infty }}{\mathit{lim}}G\left(t\right)=G_{\infty }>0,$$

where $G_{\infty }$ is the material equilibrium modulus. Then, Equation (1) takes the form presented below [31]:

$$G\left(t\right)={\int }_0^{\mathrm{\infty }}{\displaystyle \frac{\mathcal{H}\left(\tau \right)}{\tau }}{e}^{-t/\tau }d\tau +G_{\infty }.$$

(79)

Analogously, Equation (5)—basic for the direct approach and related to the modified frequency spectrum $H^M\left(v\right)$—can be rewritten as follows:

$$G\left(t\right)={\int }_0^{\mathrm{\infty }}{H}^M\left(v\right)e^{-tv}dv+G_{\infty }.$$

(80)

The relaxation spectra models $H_K^M\left(v\right)$ (9), $H_K\left(v\right)$ (10), and ${\mathcal{H}}_K\left(\tau \right)$ (14) do not require modification, while the related relaxation modulus model $G_K\left(t\right)$ (11) should be replaced by the following:

$$G_K\left(t\right)={\int }_0^{\mathrm{\infty }}{H}_K^M\left(v\right)e^{-tv}dv+G_{\infty }={\sum }_{k=1}^K{g}_k{\varphi }_k\left(t\right)+G_{\infty },$$

which, however, does not affect the identification procedure itself.

The square integral index $\mathcal{J}\left({\mathit{g}}_K\right)$, given by Equation (16) for the model $H_K^M\left(v\right)$ and by (20) for ${\mathcal{H}}_K\left(\tau \right)$, is defined as above. However, by (80) and (7) we have:

$${\int }_0^{\mathrm{\infty }}{H}^M\left(v\right)h_k\left(v\right)dv={\int }_0^{\mathrm{\infty }}{H}^M{e}^{-\alpha kv}dv=G\left(\alpha k\right)-G_{\infty }=∆G\left(\alpha k\right).$$

(81)

Therefore, by (17) and (81), the index $\mathcal{J}\left({\mathit{g}}_K\right)$ is given by the following expression (compare to (19)):

$$\mathcal{J}\left({\mathit{g}}_K\right)={\int }_0^{\mathrm{\infty }}{H^M\left(v\right)}^{2}dv-2{\sum }_{k=1}^K{g}_k∆G\left(\alpha k\right)+{\sum }_{k=1}^K{{\sum }_{m=1}^K{g}_{k}g}_m{\phi }_{km}.$$

As above, from (79) and (20), the analogous formula results for the relaxation time spectrum model. As a consequence, the integral-empirical index ${\overline{\mathcal{J}}}_K\left({\mathit{g}}_K\right)$ (24) is now as follows

$${\overline{\mathcal{J}}}_K\left({\mathit{g}}_K\right)={\int }_0^{\mathrm{\infty }}{H^M\left(v\right)}^{2}dv-2{\sum }_{k=1}^K{g}_k∆\overline{G}\left(\alpha k\right)+{\sum }_{k=1}^K{{\sum }_{m=1}^K{g}_{k}g}_m{\phi }_{km},$$

where, compare $∆G\left(\alpha k\right)$ (81), the relaxation modulus increment is defined as follows:

$$∆\overline{G}\left(\alpha k\right)=\overline{G}\left(\alpha k\right)-G_{\infty }.$$

(82)

Since real materials may relax over a very long time, two cases can occur.

• Case 1. If the duration of the relaxation test can be extended so as to experimentally record a time-constant relaxation modulus (in practice, constant stress), then $G_{\infty }$ is experimentally evaluated and the proposed identification algorithm can be simply applied by replacing the measurements $\overline{G}\left(\alpha k\right)$ with their increments $∆\overline{G}\left(\alpha k\right)$ (82) in relation to known $G_{\infty }$.
• Case 2. For identification purposes, only time-varying relaxation modulus measurements are available, i.e., the steady-state stress was not recorded during the experiment. In such a situation, non-negative $G_{\infty }$ is an additional model parameter that should be extrapolated beyond the experiment time horizon limited by the upper bound $t_K=\alpha K$. The linear-quadratic problem (32) of optimal identification needs to be reformulated, re-regularized and solved, which creates a new research problem.

However, for many materials, the equilibrium modulus is accessible by experiment; then, the algorithm of direct relaxation spectra identification can be applied with the simple modification as described above.

4. Conclusions

Summarizing, this paper addresses the relaxation spectrum identification problem in a new original way. The novelty of the paper is that it directly takes into account the unknown spectrum in the model quality index being minimized. The main result is based only on the definition of the relaxation spectrum, which relates the spectrum to the measurable relaxation modulus, and on the fact that the set of exponential functions, i.e., a kernel of the Lagrange transform constitute a basis of the space of square-integrable functions. The analytical and numerical studies demonstrated that by applying the proposed relaxation spectra models and identification algorithm, it is possible to determine the spectra models for a wide range of relaxation times and frequencies of real materials.

The concept of direct relaxation spectrum identification can be applied both for viscoelastic fluids and viscoelastic solids; however, for solid materials, a respective modification of the algorithm may be required whenever the equilibrium relaxation modulus is not available by measurement, the development of which will be the subject of further research.

It is generally accepted that the choice of respective regularization parameters is important to identify the best model. The well-studied techniques for computing a good regularization parameter such as the discrepancy principle, generalized cross-validation, and the L-curve technique have been developed for classical least-squares task and hence they cannot be directly applied here. Therefore, the regularized minimization problem (32) should be reformulated to the classic form of the linear least-squares problem. Then, the applicability of the known techniques can be verified. An alternative approach is to develop a new method of selecting the regularization parameter, specifically addressing the problem of direct spectrum identification. Although the numerical studies have shown that the simple rule based on the condition number of the basic matrix for the linear-quadratic identification problem is sufficient in many cases, the example of a two-mode Gaussian-like spectrum motivates the search for a better rule for the regularization parameter selection, dedicated for this specific identification task. This will be the subject of further research.

The impact of the molecular weight distributions (MWD) on the viscoelastic properties is intensively studied in polymer rheology. Generic analytical formulas describing the relationship between MWD and the relaxation time spectrum are known. Future research directions may include the determination of the MWD, which can be obtained from the relaxation time spectrum model and recovered from experimental results by the proposed method.