{'fraction_non_alphanumeric': 0.05449917652343165, 'fraction_numerical': 0.012227379348205818, 'mean_word_length': 4.179116050659085, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Identification of a relaxation kernel using two boundary measures

We consider a distributed system with persistent memory of a type which is often encountered in viscoelasticity or in the study of diffusion processes with memory. The relaxation kernel, i.e. the kernel of the memory term, is scarcely known from first principles, and it has to be inferred from experiments taken on samples of the material. We prove that two boundary measures give a linear Volterra integral equation of the first kind for the unknown kernel. Hence, with two measures, the identification of the kernel, which in principle is a nonlinear problem, is reduced to the solution of a mildly ill posed but linear problem.


Introduction
Let w = w(x, t), t > 0 and x ∈ Ω (a bounded region with smooth boundary). The equation with persistent memory (∆ is the laplacian) is encountered in thermodynamics of systems with memory (w is the temperature) in nonfickian diffusion (w is the concentration) or in viscoelasticity. In this case w(t) represents the displacement and the equation is usually written as (1. 2) The vector is the (density of the) flux at time t and position x in thermodynamics or nonfickian diffusion. In viscoelasticity, q ′ (x, t) is the traction. Let n(x) be the exterior normal to ∂Ω at the point x. The products q(x, t) · n(x) or q ′ (x, t) · n(x), x ∈ ∂Ω, represent respectively the flux throughout the boundary or the traction on the boundary. In both the cases, they are observable quantities. Especially in viscoelasticity, N(t) is called the relaxation kernel. In this paper we assume that it is a smooth function (see below for details) and in the applications we mentioned above, it must be N(0) > 0.
It is a fact that the kernel N(t) is a material property. Only its overall properties can be deduced from general principles (see for example [9]) and the specific function N(t) has to be identified from experiments taken on the system. The "experiments" usually considered consist in the observation of the traction/flux on the boundary of a suitable sample of the material, subject to known excitations. The memory kernel N(t) is identified from such observations.

Remark 1
The integrals in (1.1) should extend from the time the material has been "created", conceivably −∞. We assume that the measures are taken after a certain time t 0 sufficiently far from the "creation" of the material and that the material was kept undisturbed for a long time, so that the effect of the memory of the previous history is negligible. This time t 0 we rename 0.
The shape of sample of the material usually does not influence the memory kernel otherwise, N would depend also on the position x (see [4]) so that we can try to determine N(t) from samples which have simple geometries. As in [4] and most of the references (commented in Sect. 1.2), we shall use samples which are in the shape of a slab, so that we can assume that Ω is a segment and we can normalize Ω = (0, π) .

The identification of N(t) is a nonlinear problem since if N(t)
is not yet known, the function w(t) cannot be computed and the problem consists in the identification of the pair (N(t), w(x, t)). This pair enters in a nonlinear way in the equation. In spite of this, we shall see that using two measures we can construct a linear Volterra integral equation solely for N(t). So, the identification of N(t) is reduced to a linear problem. This is a first difference with existing literature. A second difference is as follows: in many papers, especially in engineering applications, the identification of N(t) is done as follows: it is assumed that the kernel belongs to a certain class of functions which depend on a finite number of parameters. The special kernel N(t) assumed to be the "true" kernel is the one which minimizes the discrepancy between a finite number of experimental data and the corresponding theoretically computed values. In contrast with this, we don't assume any known a priori information on the functional form of N(t), a part smoothness and the condition N(0) > 0.
Due to different applications and terminology, for consistency and clarity we state that w will be called "temperature" and q will be called the "flux" (note that the usual minus sign in the definition of the flux or the traction has no influence in the following arguments and it is ignored).

1.1
The assumption on N (t) and informal description of the procedure The first of our assumptions is that the kernel N(t) is of class C 3 (so that the memory kernel M(t) = N ′ (t) of the viscoelastic system is of class C 2 as often assumed in this kind of problems, see [1]). It is seen from the form (1.2) of Eq. (1.1) that this equation is a perturbation of the wave equation since we explicitly assume c 2 = N(0) > 0. A consequence of N(0) > 0 is that signals propagate in the body with finite velocity c = N(0). So, whether this assumption holds is easily checked and the value of c is easily computed: just measure the velocity of propagation of the signals in the sample of the material.
Our identification procedure applies in the case when c > 0. Solely to present simpler formulas, we assume that a first, simple, step is the identification of the value c = N(0). We assume that this step has already been done. Then, the change of the time scale t → t/c transform N(t) to the new kernel N 1 (t) = (1/c)N(t/c). Note that N 1 (0) = 1. So, we assume: • the velocity c > 0 of signal propagation has been identified.
Now we describe informally the identification procedure. We recall the goal: we fix any time T and we want to identify N(t) for t ∈ [0, T ] (we already know N(0) = 1).
1. We perform the following two measures, both on intervals of duration T (it is not restrictive to rename such intervals as [0, T ]): (a) the temperature w(x, t) at the ends x = 0 and x = π is kept equal zero, w(0, t) = 0, w(π, t) = 0. We impose a special initial condition w(x, 0) = ξ(x) and we read the flux y ξ (t) = q(π, t) (note the lower case y). We use a suitable initial condition ξ 0 (a ramp) in order to identify an intermediate function K(t) which is then used in the identification of the memory kernel.
(b) We impose null initial condition and w(π, t) = 0 but now we excite the system from the left end: we impose w(0, t) = f (t). We measure q(π, t). For clarity, the flux due to the boundary displacement f (with zero initial condition) is denoted Y (t) = Y f (t) (note the upper case).
In this second measure we don't need to use any special function f (t), in particular f (t) is not related to the initial condition used in the first measure. But in practice the function f (t) will be smooth. For continuity we have to take f (0) = 0 and f (t) "saturates" in a short time, i.e. it approaches lim t→+∞ f (t) (which in practice is finite) quite fast.
2. We derive an algorithm for the identification of N(t) using these two observations.
We note that the first step might be done by imposing w(x, t) = 0 and w ′ (x, 0) = η(x), with a special initial velocity η(x) (integrating (1.2) we see that in this case Eq. (1.1) must have an additive term η(x) on the right hand side). For definiteness, we study η = 0 and ξ = 0 but in practical applications to viscoelasticity ξ = 0 and η = 0 might be easier to impose. Now we describe the organization of the paper. In Sect. 1.2 we discuss few references on the identification problem and in Sect. 2 we present the preliminary properties of system (1.1) we shall needed. The identification algorithm is in Sect. 3.
It is possible to identify N(t) using a similar algorithm, which uses two measures with different initial conditions (and zero temperature at the boundary). We believe that this algorithm is not practical since to impose initial condition is more difficult then to impose a boundary temperature (or deformation). However, for completeness, this second algorithm is described in Sect. 3.3.

Comments on previous references
Due to the importance of this problem, reference on kernel identification is enormous and both Engineers and Mathematicians addressed the problem, along different and scarcely related lines. Most of the available results can be described as follows. An approach to kernel identification, seemingly first proposed in [13] and then followed by a wealth of related papers (see also the books [1,14] for references) is as follows: it is noted that the identification of N(t) on the basis of additional measures is in fact the identification of the pair (w(t), N(t)) and w(t) solves an integrodifferential equation. Now the trick is to compute derivatives of the measured output (of course when the data are sufficiently smooth) and to derive a second integrodifferential equation for N(t). In this way we get a nonlinear integrodifferential equation for the pair (w(t), N(t)) and, in general using fixed point theorems, it is possible to prove existence of solutions of this systems.
This leads to the identification in particular of N(t) (according to the assumptions, either locally, for a short time, or globally, as in [15]).
The identification requires only one measure, taken assuming known (and smooth) initial condition (or affine term).
The approach that is most common in engineering literature is as follows. It is assumed that the kernel N(t) belongs to a certain class of functions, identified by a finite number of parameters. The classes that are most often considered are Dirichlet sums, i.e.
as for example in [4,8,24,25] (here a ≥ 0, b ≥ 0) or a combination of Abel kernels sometimes a combination of both (see for example [5,10]). The class of the kernels (1.3) or (1.4) have a "physical" sense since they are encountered when solving ordinary differential equations of integer or fractional order and these classes of kernels are reported to fit sufficiently well experimental data, but not in every application, see [23].
In order to identify the parameters, a sample of the material, usually in the form of a slab, is excited using a known initial or boundary deformation. The stress relaxation is measured on (a part of the) boundary and the parameters (a k , b k , sometime also the exponents β k ) are chosen so to fit better the experimental measures (usually minimizing a quadratic index). This can be done using a single measure. Multiple measures can be used to improve accuracy or to remove spurious frequencies.
Different algorithms works in the frequency domain, using Fourier or Laplace transformations, as for example in [12].
The use of a single measure is common to the algorithm presented in [13,15] and related papers. Here instead we assume that two different measure are performed and, from these measures, the kernel N(t) is reconstructed. We don't assume that the kernel belongs to any specific class of functions, a part regularity and the physical property N(0) > 0. In fact, no further properties of the kernel are in principle needed for the reconstruction algorithms (not even positivity for every time!) in spite of the fact that such conditions are imposed by thermodynamics and can be used to improve the results, see Sect. 3.4. Note however that we are assuming that N(t) is smooth for t ≥ 0 and this rules out Abel kernels. The extension of our algorithm to Abel kernels will be pursued in future researches.
We consider the sequence of the functions It is known form the theory of Fourier series that {φ n } is an orthonormal basis of L 2 (0, π) and that We expand the initial condition ξ and the solution in sine series: This is a scalar integrodifferential equation whose solution can be represented as w n (t) = ξ n z n (t) + n The transformation (ξ, f ) → w is linear and continuous from We need the following observation: and every n we have: Proof. This inequality follows for example from the equality (5.4) in [17] (see also [18,16]). We must be careful with the notation. With α = −N ′ (0)/2 letÑ (t) = e 2αt N(t) so thatÑ(0) = 1 andÑ ′ (0) = 0. The functions denoted z n (t) and N(t) in [17, formula (5.4)] are respectively e 2αt z n (t) andÑ (t) = e 2αt N(t) (see [17, formula (5.1)]). So, we have the following equality for the function z n (t) in (2.1) (here β n = √ n 2 − α 2 ): We compute the derivatives of both the sides, usingÑ ′ (0) = 0. We get the equality The result follows by dividing both the sides with n and nothing that if t ∈ [0, T ] then there exists M such that | sin β n t − sin t| ≤ M n .
We are going to use these properties of the system for the identification of the kernel N(t), assumed unknown. Note that these same properties have also been used to solve a problem of source identification, see [19,20].

The identification procedure
Now we give the details of the two steps of the identification procedure described in Sect. 1.1.

So we have
Using the inequality (2.2), we see that y ξ (t) ∈ L 2 (0, T ) for every T > 0 and every ξ ∈ L 2 (0, π). If in particular we choose then the series in (3.2) is The function we wanted to identify is This is the first measure we need in the identification process.

Measure of the flux due to a boundary excitation
Now we combine the identification of K(t) obtained in the previous section with an observation of the flux at x = π due to a nonzero temperature applied to the left hand x = 0 (and zero initial condition). Hence we consider the following system: (3.5) We choose any C 1 input f such that f (0) = 0, (this is consistent with the fact that the initial condition is zero). So we have and we explicitly assume We use Here δ(x) is the Dirac's delta (supported at x = 0) and so when observing the flux at x = π we get where, from (3.4), We have seen in the previous section that K(t) can be estimated from the measure of one boundary flux. So, at this point, K(t) is a known function. Formula (3.7) explains its significance: K(t) − 1/2 t 0 N(s) ds is the impulse response of the system with boundary input f at x = 0 and the flux observed at x = π.
Remark 3 In fact, the δ-function is an artifact, due to the fact that formula (3.6) gives an odd extension of the solution which is periodic of period 2π and so the formula introduces a jump at t = 0 (because f (t) = 0) and this produces a δ function at x = 0 and also at every multiple of 2π. This fact has been discussed in [2,3].
The function N(t) can be identified from formula (3.7) which is a Volterra integral equation of the first kind for Although this may not be the best way from the point of view of numerical computation, we may note that if g ∈ C 1 , hence f ∈ C 2 , then Y (t) is differentiable and upon computing the derivatives of both the sides of (3.7) we get: We have chosen an input signal f such that f ′ (0) = g(0) = 0. So, this is a Volterra integral equation of the second kind (but note that the right hand side requires the numerical computation of Y ′ (t)). Once this Volterra integral equation has been solved, a last step of numerical differentiation identifies N(t).

Remark 4 We note that
• If g(0) = 1 and g ′ (t) = 0 then we get an explicit formula for M(t) but in this case f (t) = t. A steadily increasing temperature cannot be applied. In practice, f (t) has to saturate, as for example f (t) = t/(t+1) or 1 − e −t .
• In practice it does not seems difficult to impose a variable temperature (or traction) on the boundary. It may be more difficult to impose a special initial condition but we note that the initial condition we need in the identification of K is a ramp, which may be achieved by letting the equilibrium be reached between two different temperatures of the ends of the sample.

A different identification procedure
We present now a second identification algorithm which uses two different measures, obtained using two different, but related, initial conditions. In the algorithm we proposed above, the most questionable step seems the use of a special initial condition. It is even more difficult to impose also a second, related, initial temperature. So, the procedure we present now may not have a practical interest, but it seems mathematically interesting.
We impose an initial temperature ξ(x) = +∞ n=1 ξ n sin nx and we use (3.2): Let us compute the convolution N * y ξ (we recall that the kernel N(t) is still unidentified): The functions y ξ and y η can be measured and so we find a Volterra integral equation for the unknown function N(t). The Volterra integral equation is of the second kind, i.e. its solution is a well posed problem, if ξ is such that +∞ n=1 (−1) n ξ n n = 0 .
In spite of this, the identification problem is ill posed since the procedure requires the numerical computation of the derivative of y η (t).

A toy simulation
Whether the identification procedure we presented is really useful is a matter to be settled by the Engineers on the basis of concrete experiments. Here we present a toy simulation on the basis of theoretical computations. We assume that N(t) is a linear combination of exponentials, as it is often assumed in applications, and for simplicity we confine ourselves to the case that N(t) is the combination of three exponentials, where the coefficients c k and the exponents τ k are positive. Then, the equation of z n (t) is It is easy to reduce the computation of z n (t) to the solution of an ordinary differential equation. In fact, z n (t) is the first component of the solution (z, w 1 , w 2 , w 3 ) of the following problem: So, the functions z n (t) can be easily computed and we can compute also Y f (t) and K(t). Then we choose f (t) = 1 − e −t and we reconstruct the kernels M(t) = t 0 N(s) ds and N(t). Before we present the simulation we note that: • the reconstruction of the relaxation kernel is in two steps: first we identify M(t) and then a next step of numerical differentiation gives N(t). The identification of M(t) is a deconvolution problem which can be solved, essentially with the same results, with Tikonov regularization (as described in [11]) or Lavrientev method (see [6,7]). Both these methods introduce fast oscillations but it is known in this problem that the oscillations are an artifact due to the amplification of the noise. In fact, M(t) is increasing (since for physical reasons N(t) ≥ 0) and, in the physical case that N(t) is decreasing, M(t) is concave. Hence these oscillations can be attenuated for example by computing the average of the reconstruction of M(t) on nearby steps. This has been done in the simulation below.
• in the simulation we present below we choose to reconstruct M(t) using Lavrentev method (followed by an average, to remove fast oscillations). Instead we use Tikonov regularization in the next step of the reconstruction of N(t) from M(t) since Lavrentev method is more efficient when it is known that the function to be reconstructed is 0 at the initial time (and since Lavrentev method is more sensitive to disturbances).
• Also the function N(t) can be averaged to reduce the effect of errors since in practice we know that it is a decreasing (and often a convex) function. Taking into account the known value N(0) this last step (not shown in the plots) will also improve the reconstruction near t = 0.
The special case we test is as follows: the time interval is [0, 5] and we assume 10 measures every time unit. We choose N(t) = (1/10)e −t/2 + (1/5)e −2t + (1/2)e −3t (3.12) (i.e. we normalize so to have N(0) = 1). The functions K(t) and Y f (t) are approximated (and corrupted by an error of 1% at every "measure"). Then, we "reconstruct" both M(t) and N(t). Figure 1 (left) presents the plots of M(t) and its reconstruction from Lavrientev method (the plot where the amplification of the disturbances is evident) and its average (at every step tree values on the left and on the right have been averaged). The averaged reconstruction of M(t) is fed to the Tikonov algorithm in order to compute the numerical derivative of M(t), i.e. N(t), which is plotted on the right. The penalization parameters of the Lavrientev and Tikonov algoritms are respectively 0.1 (equal to the interval of time between consecutive measures) and 0.01.