Skip to content
Publicly Available Published by De Gruyter November 8, 2019

Regularized Collocation in Distribution of Diffusion Times Applied to Electrochemical Impedance Spectroscopy

  • Sergei V. Pereverzev ORCID logo , Sergiy G. Solodky , Vitalii B. Vasylyk ORCID logo EMAIL logo and Mark Žic

Abstract

This paper is inspired by recently proposed approach for interpreting data of Electrochemical Impedance Spectroscopy (EIS) in terms of Distribution of Diffusion Times (DDT). Such an interpretation requires to solve a Fredholm integral equation of the first kind, which may have a non-square-integrable kernel. We consider a class of equations with above-mentioned peculiarity and propose to regularize them in weighted functional spaces. One more issue associated with DDT-problem is that EIS data are available only for a finite number of frequencies. Therefore, a regularization should unavoidably be combined with a collocation. In this paper we analyze a regularized collocation in weighted spaces and propose a scheme for its numerical implementation. The performance of the proposed scheme is illustrated by numerical experiments with synthetic data mimicking EIS measurements.

MSC 2010: 65R32; 92E99

1 Introduction

Electrochemical Impedance Spectroscopy (EIS) has become a powerful technique to determine resistive and dielectric properties of materials [1]. In EIS, the angular excitation frequency ω of the applied potential is varied over a range of frequencies [ ω 1 , ω m ] , while the impedance (both resistance and inductance) is measured.

The traditional interpretation of EIS data is given in terms of a distribution of relaxation times (DRT) that maps impedance spectra on a function containing the timescale characteristics of the system under consideration [3]. As it is pointed out in [17], DRT analysis is intended only to represent high-frequency region of the impedance spectra (i.e., interfacial charging and Faradaic reactions).

At the same time, low-frequency impedance may contain geometrical information that can be used to infer microstructural statistics of heterogeneous and nanostructural materials. A theory of such impedance for random heterogeneous materials has been recently proposed in [17]. It is based on the concept of a function of distribution of diffusion times (DDT) and is given by Fredholm integral equations

(1.1) y ( ω ) = 0 K ( ω , τ ) p ( τ ) 𝑑 τ ,

where p ( τ ) is DDT, y ( ω ) is the collective diffusion admittance, and K ( ω , τ ) is the so-called finite-length diffusion model that represents individual diffusion paths depending on boundary conditions and symmetries.

When inverting the impedance spectrum y ( ω ) to obtain DDT p ( τ ) , one should be aware that this class of inversion problems is known to be mathematically ill-posed and, thus, should be treated with a regularization.

Moreover, since only a finite number of spectrum measurements y ( ω i ) , i = 1 , 2 , , m , may actually be performed, the above-mentioned equations (1.1) appear as after application of a collocation method that needs to be combined with the regularization.

In EIS literature such a combination has been discussed in the context of inversion of DRT equations [15, 16, 5, 22, 2, 4, 24].

One of the essential differences between DRT and DDT approaches is that, in contrast to DRT equation, the kernels K ( ω , τ ) in (1.1) may not be square-integrable on [ ω 1 , ω m ] × [ 0 , ) . For example, the real and imaginary parts of the kernel K ( ω , τ ) = i ω τ tanh ( i ω τ ) suggested in [17] for blocking boundary conditions and planar symmetry, asymptotically behave like O ( ω τ ) and clearly are not square-integrable on [ ω 1 , ω m ] × [ 0 , ) . Then for equations (1.1) with such kernels K ( ω , τ ) a straightforward application of Tikhonov regularization combined with the discretization by the trapezoidal rule, as it is adopted in [17], is not justified by the Regularization theory.

Therefore, the goal of the present study is to theoretically analyze the regularized collocation of DDT-problems. In the next section we consider somewhat more general class of equations (1.1), containing DDT equations as special cases, and study their regularization in weighted L 2 -spaces. Then we interpret our results in the context of DDT-problem and illustrate them by numerical examples with synthetic data.

2 A Class of Integral Equations with Non-integrable Kernels and Their Regularization

Let L 2 , = L 2 ( 0 , ) be the Hilbert space of functions f that are square-integrable over half line [ 0 , ) ,

f L 2 , 2 = 0 | f ( τ ) | 2 𝑑 τ < .

By L 2 β = L 2 β ( 0 , ) , β > 0 , we denote a weighted subspace of L 2 , with the norm

f L 2 β := ( 0 ( 1 + τ ) 2 β | f ( τ ) | 2 𝑑 τ ) 1 2 < .

In the sequel the integral operators

K p ( ω ) = 0 K ( ω , τ ) p ( τ ) 𝑑 τ , ω [ ω 1 , ω m ] ,

in (1.1) will be considered as acting from L 2 β into the Hilbert space L 2 , Ω = L 2 ( ω 1 , ω m ) of functions y that are square-integrable over the interval ( ω 1 , ω m ) ,

y L 2 , Ω 2 = ω 1 ω m | y ( ω ) | 2 𝑑 ω < .

We will need also an embedding operator J of L 2 β into L 2 , such that for any f L 2 β we have

J f ( τ ) = ( 1 + τ ) β f ( τ ) L 2 , .

The adjoint of J is the operator J * : L 2 , L 2 β assigning to a function g L 2 , the function

J * g ( τ ) = ( 1 + τ ) - β g ( τ ) L 2 β .

Then equation (1.1) can be rewritten as

(2.1) H x = y ,

where

x ( τ ) = J p ( τ ) = ( 1 + τ ) β p ( τ ) L 2 ,

and H is an integral operator

(2.2) H x ( ω ) = 0 H ( ω , τ ) x ( τ ) 𝑑 τ , ω [ ω 1 , ω m ] ,

with a kernel

(2.3) H ( ω , τ ) = K ( ω , τ ) ( 1 + τ ) - β .

In what follows we assume that β can be chosen in such a way that the kernel H ( ω , τ ) is square-integrable on [ ω 1 , ω m ] × [ 0 , ) . For example, as it will be shown in the next section, the kernel K ( ω , τ ) = i ω τ tanh ( i ω τ ) considered in the paper [17] providing the motivation for the present study, allows the choice β > 3 2 to make the corresponding kernels (2.3) square-integrable. Note that the operator (2.2)–(2.3) with the above chosen β is a compact operator from L 2 , into L 2 , Ω . Therefore, potentially one could employ standard regularization techniques to treat the equation (2.1), but since in the application we have in mind, only a finite number of measurements y ( ω j ) , j = 1 , 2 , , m , is available, equation (2.1) is given only in a collocated form.

Formally, it can be written in terms of the sampling operator S Ω g = ( g ( ω 1 ) , g ( ω 2 ) , , g ( ω m ) ) m as

(2.4) S Ω H x = S Ω y ,

where the image space m of the operator S Ω is the Euclidean space of vectors u = ( u 1 , , u m ) and v = ( v 1 , , v m ) , equipped with the inner product

u , v m = j = 1 m γ j u j v j , γ j > 0 , j = 1 , , m .

Regularizing (2.4) by Tikhonov method [19], we arrive at the equation

(2.5) α x + H Ω * H Ω x = H Ω * S Ω y ,

where

H Ω = S Ω H = S Ω K J * : L 2 , m , H Ω x = ( 0 K ( ω i , τ ) ( 1 + τ ) - β x ( τ ) 𝑑 τ ) i = 1 m .

The solution x α , m of the regularized equation (2.5) admits the representation

x α , m = ( α I + H Ω * H Ω ) - 1 H Ω * S Ω y .

Moreover, it is easy to check that

H Ω * u ( τ ) = j = 1 m γ j u j K ( ω j , τ ) ( 1 + τ ) - β .

Then the regularized solution has the form

(2.6) x α , m ( τ ) = j = 1 m x α j K ( ω j , τ ) ( 1 + τ ) - β ,

where the coefficients x α j solve the following system of linear algebraic equations:

(2.7) α x α j + γ j l = 1 m x α l 0 K ( ω j , v ) K ( ω l , v ) ( 1 + v ) - 2 β 𝑑 v = γ j y ( ω j ) , j = 1 , , m .

In practice, however, instead of

S Ω y = ( y ( ω 1 ) , y ( ω 2 ) , , y ( ω m ) )

only noisy measurements ( y δ ( ω 1 ) , y δ ( ω 2 ) , , y δ ( ω m ) ) of y ( ω ) are available. Formally, a vector of such measurements can be seen as the image S Ω y δ of some perturbed version y δ ( ω ) of data function y ( ω ) . Then the level of perturbations can be estimated as follows:

S Ω ( y - y δ ) m 2 := j = 1 m γ j | y ( ω j ) - y δ ( ω j ) | 2 δ 2 .

The regularized approximation x α , m δ based on noisy data is given as

(2.8) x α , m δ = ( α I + H Ω * H Ω ) - 1 H Ω * S Ω y δ .

It also has the form (2.6) with the coefficients that solve system (2.7), where y ( ω j ) is substituted for y δ ( ω j ) .

Recall now that in the regularization theory [10, 11, 12, 13] the concept of source conditions is used to quantify the accuracy of an approximate solution. According to this concept [11] the Moore–Penrose generalized solution x = x of equation (2.1), assuming it exists, can always be represented in terms of an operator-valued non-decreasing index function φ such that φ ( 0 ) = 0 and

(2.9) x = φ ( H * H ) v , v L 2 , ρ .

Moreover, as it is explained in [14], when considered collocated Tikhonov regularization, such as (2.5), it is natural to assume that in (2.9) the index function φ is operator monotone on an interval containing the spectra of operators H * H , H Ω * H Ω . Then from [10, Proposition 2.22] it follows that

(2.10) φ ( H * H ) - φ ( H Ω * H Ω ) L 2 , L 2 , d φ ( H * H - H Ω * H Ω L 2 , L 2 , ) ,

where the coefficient d depends only on φ. Moreover, there exists a φ-dependent coefficient χ φ such that

(2.11) sup λ > 0 α α + λ φ ( λ ) χ φ φ ( α ) .

Note that examples of operator monotone index functions include the functions φ ( λ ) = λ ν , 0 < ν 1 , and φ ( λ ) = log - μ ( 1 λ ) , λ ( 0 , 1 ) , μ > 0 , which traditionally used in the regularization theory to describe the solution smoothness in terms of source conditions.

Now, using the argument similar to that in [14], we can prove the following statement

Theorem 1.

Assume that source condition (2.9) with an operator monotone index function is satisfied. Then

x - x α , m δ L 2 , δ 2 α + C ( φ ( α ) + φ ( H * H - H Ω * H Ω L 2 , L 2 , ) ) ,

where the coefficient C depends only on φ and ρ involved in condition (2.9).

Proof.

To estimate the error x - x α , m δ L 2 , , we write the decomposition

x - x α , m δ = ( x - x α , m ) + ( x α , m - x α , m δ ) .

For the first difference we have

x - x α , m = x - ( α I + H Ω * H Ω ) - 1 H Ω * S Ω y = ( I - ( α I + H Ω * H Ω ) - 1 H Ω * H Ω ) x = T 1 + T 2 ,

where

T 1 = ( I - ( α I + H Ω * H Ω ) - 1 H Ω * H Ω ) φ ( H Ω * H Ω ) v ,
T 2 = ( I - ( α I + H Ω * H Ω ) - 1 H Ω * H Ω ) ( φ ( H * H ) - φ ( H Ω * H Ω ) ) v .

By virtue of (2.10), (2.11) we have

T 1 L 2 , ( I - ( α I + H Ω * H Ω ) - 1 H Ω * H Ω ) φ ( H Ω * H Ω ) L 2 , L 2 , v L 2 ,
ρ sup λ α α + λ φ ( λ )
χ φ ρ φ ( α ) ,
T 2 L 2 , I - ( α I + H Ω * H Ω ) - 1 H Ω * H Ω L 2 , L 2 , φ ( H * H ) - φ ( H Ω * H Ω ) L 2 , L 2 , v L 2 ,
d ρ sup λ α α + λ φ ( H * H - H Ω * H Ω L 2 , L 2 , )
d ρ φ ( H * H - H Ω * H Ω L 2 , L 2 , ) .

It remains to estimate the difference x α , m - x α , m δ . To this end, we use the polar decomposition (see, for example, [21, p. 35])

(2.12) H Ω * = | H Ω | U * ,

where | H Ω | = ( H Ω * H Ω ) 1 2 : L 2 , L 2 , is the operator module of H Ω and U * : m L 2 , is the partial isometry operator, U * m L 2 , = 1 . Then (2.12) allows us to continue as follows:

x α , m - x α , m δ L 2 , = ( α I + H Ω * H Ω ) - 1 H Ω * S Ω ( y - y δ ) L 2 ,
( α I + H Ω * H Ω ) - 1 H Ω * m L 2 , S Ω ( y - y δ ) m
( α I + H Ω * H Ω ) - 1 | H Ω | L 2 , L 2 , U * m L 2 , S Ω ( y - y δ ) m
δ 2 α .

Thus, we obtain the desired estimate. ∎

Note that the term φ ( H * H - H Ω * H Ω L 2 , L 2 , ) of the error bound given by Theorem 1 reflects a contribution of the collocation of operator H to the total error. In the sequel it is assumed that

(2.13) H * H - H Ω * H Ω L 2 , L 2 , α .

In view of Theorem 1 such an assumption ensures that the error caused by the collocation of H does not dominate the regularization error.

To balance all components of the error bound of Theorem 1 we consider the function

Θ ( λ ) = λ φ ( λ ) , λ [ 0 , H L 2 , L 2 , Ω 2 ] .

Theorem 2.

Let α = Θ - 1 ( δ ) and relation (2.13) hold. Then under the conditions of Theorem 1 we have

(2.14) x - x α , m δ L 2 , = O ( φ ( Θ - 1 ( δ ) ) ) .

Proof.

It is clear that

δ Θ - 1 ( δ ) = φ ( Θ - 1 ( δ ) ) .

Then, by virtue of Theorem 1 and (2.13), we obtain

x - x α , m δ L 2 , φ ( Θ - 1 ( δ ) ) 2 + 2 C φ ( Θ - 1 ( δ ) ) = O ( φ ( Θ - 1 ( δ ) ) ) .

Note that, as it follows, for example, from [12], for given noise level δ and the source condition (2.9) the order (2.14) cannot be improved in general.

From the above discussion it follows that condition (2.13) is important for the performance of the regularized collocation. The fulfillment of (2.13) depends on the kernel (2.3) and the sampling operator S Ω . The properties of the later one are encoded in the distribution of sampling points ω j , j = 1 , 2 , , m , and the weights γ j , j = 1 , 2 , , m , inducing the inner product , m . We will interpret them as knots and coefficients of a numerical integration formula and adopt the following assumption.

Assumption 1.

Assume that the quadrature formula

(2.15) ω 1 ω m f ( ω ) 𝑑 ω j = 1 m γ j f ( ω j )

is exact for polynomials of the 1st degree and for any function f having continuous derivatives f ( r ) in [ ω 1 , ω m ] it holds

(2.16) | ω 1 ω m f ( ω ) 𝑑 ω - j = 1 m γ j f ( ω j ) | c r , Ω h r max ω 1 ω ω m | f ( r ) ( ω ) | ,

where the coefficient c r , Ω > 0 depends on r, Ω and

h = sup ω [ ω 1 , ω m ] min ω j | ω - ω j |

is the so-called grid norm.

Remark 1.

The simplest example of formulae (2.15) satisfying Assumption 1 is the trapezoidal formula (for r = 2 ). More sophisticated constructions can be found in [8].

In the sequel we restrict ourselves to equations (1.1) with the kernels K ( ω , τ ) admitting the representation

(2.17) K ( ω , τ ) = ω 1 2 ε 1 ( τ ) + ε 2 ( ω , τ ) ,

where the functions ε 1 and ε 2 are such that

  1. ( 1 + τ ) - β ε 1 ( τ ) L 2 , ,

  2. i ε 2 ( ω , τ ) ω i , i ( ω 1 2 ε 2 ( ω , τ ) ) ω i , i = 0 , 1 , , r , are continuous on [ ω 1 , ω m ] × [ 0 , ) and

    sup ω , τ | i ε 2 ( ω , τ ) ω i | < , sup ω , τ | i ( ω 1 2 ε 2 ( ω , τ ) ) ω i | < .

In the next section we will show that representation (2.17) is relevant in the context of DDT-problem motivating the present study.

Consider the following auxiliary functions:

g 0 ( τ ) = ( 1 + τ ) - β max ω | r ( ω 1 2 ε 2 ( ω , τ ) ) ω r | ,
g 1 ( ω ) = i = 0 r C i r | r - i ω 1 2 ω r - i | ( 0 | i ε 2 ( ω , u ) ω i ( 1 + u ) - β | 2 𝑑 u ) 1 2 ,
g 2 ( ω , τ ) = i = 0 r C i r | r - i ε 2 ( ω , τ ) ω r - i | ( 0 | i ε 2 ( ω , u ) ω i ( 1 + u ) - β | 2 𝑑 u ) 1 2 ,

where C i r are binomial coefficients. It is clear that g 0 , g 1 , g 2 are continuous and bounded on [ ω 1 , ω m ] × [ 0 , ) . Moreover, it is easy to see that

( 1 + τ ) - β max ω g 2 ( ω , τ ) L 2 , .

Theorem 3.

Let Assumption 1 and conditions (a)–(b) be fulfilled. Then for operators defined by (2.2), (2.3), (2.17) the following bound holds:

H * H - H Ω * H Ω L 2 , L 2 , c ¯ c r , Ω h r ,

where

c ¯ 2 = ( 1 + ) - β ε 1 L 2 , 2 ( g 0 L 2 , 2 + max ω | g 1 ( ω ) | 2 ) + ( 1 + ) - β max ω g 2 ( ω , ) L 2 , 2 .

Proof.

First of all we find a presentation for H * H - H Ω * H Ω . Recall that

H x ( ω ) = K J * x ( ω ) = 0 K ( ω , u ) ( 1 + u ) - β x ( u ) 𝑑 u

and

H * g ( τ ) = J K * g ( τ ) = ( 1 + τ ) β K * g ( τ ) = ( 1 + τ ) - β ω 1 ω m K ( ω , τ ) g ( ω ) 𝑑 ω .

Then

H * H x ( τ ) = H * ( 0 K ( τ , u ) ( 1 + u ) - β x ( u ) 𝑑 u )
= ( 1 + τ ) - β ω 1 ω m K ( ω , τ ) 0 K ( ω , u ) ( 1 + u ) - β x ( u ) 𝑑 u 𝑑 ω .

On the other hand

H Ω x = ( 0 K ( ω i , u ) ( 1 + u ) - β x ( u ) 𝑑 u ) i = 1 m ,
H Ω * g ( τ ) = j = 1 m γ j g j K ( ω j , τ ) ( 1 + τ ) - β ,
H Ω * H Ω x ( τ ) = j = 1 m γ j K ( ω j , τ ) ( 1 + τ ) - β 0 K ( ω j , u ) ( 1 + u ) - β x ( u ) 𝑑 u .

Therefore,

( H * H - H Ω * H Ω ) x ( τ ) = ( 1 + τ ) - β ( ω 1 ω m K ( ω , τ ) 0 K ( ω , u ) ( 1 + u ) - β x ( u ) d u d ω
- j = 1 m γ j K ( ω j , τ ) 0 K ( ω j , u ) ( 1 + u ) - β x ( u ) d u ) .

To simplify the calculations, we introduce the notation

(2.18) F ( ω , τ ) = K ( ω , τ ) 0 K ( ω , u ) ( 1 + u ) - β x ( u ) 𝑑 u .

Then

(2.19) H * H - H Ω * H Ω L 2 , L 2 , 2 = 0 ( 1 + τ ) - 2 β ( ω 1 ω m F ( ω , τ ) 𝑑 ω - j = 1 m γ j F ( ω j , τ ) ) 2 𝑑 τ .

By virtue of (2.17), (2.18), we have

F ( ω , τ ) = ( ω 1 2 ε 1 ( τ ) + ε 2 ( ω , τ ) ) 0 ( ω 1 2 ε 1 ( u ) + ε 2 ( ω , u ) ) ( 1 + u ) - β x ( u ) 𝑑 u
= F 1 ( ω , τ ) + F 2 ( ω , τ ) + F 3 ( ω , τ ) + F 4 ( ω , τ ) ,

where

F 1 ( ω , τ ) = ω ε 1 ( τ ) 0 ε 1 ( u ) ( 1 + u ) - β x ( u ) 𝑑 u ,
F 2 ( ω , τ ) = ω 1 2 ε 2 ( ω , τ ) 0 ε 1 ( u ) ( 1 + u ) - β x ( u ) 𝑑 u ,
F 3 ( ω , τ ) = ω 1 2 ε 1 ( τ ) 0 ε 2 ( ω , u ) ( 1 + u ) - β x ( u ) 𝑑 u ,
F 4 ( ω , τ ) = ε 2 ( ω , τ ) 0 ε 2 ( ω , u ) ( 1 + u ) - β x ( u ) 𝑑 u .

Thus, we have

( H * H - H Ω * H Ω ) x ( τ ) = ( 1 + τ ) - β l = 1 4 ( ω 1 ω m F l ( ω , τ ) 𝑑 ω - j = 1 m γ j F l ( ω j , τ ) )

with

ω 1 ω m F 1 ( ω , τ ) 𝑑 ω - j = 1 m γ j F 1 ( ω j , τ ) = ε 1 ( τ ) 0 ε 1 ( u ) ( 1 + u ) - β x ( u ) 𝑑 u ( ω 1 ω m ω 𝑑 ω - j = 1 m γ j ω j ) ,
ω 1 ω m F 2 ( ω , τ ) 𝑑 ω - j = 1 m γ j F 2 ( ω j , τ ) = 0 ε 1 ( u ) ( 1 + u ) - β x ( u ) 𝑑 u ( ω 1 ω m ω 1 2 ε 2 ( ω , τ ) 𝑑 ω - j = 1 m γ j ω j 1 2 ε 2 ( ω j , τ ) ) ,
ω 1 ω m F 3 ( ω , τ ) d ω - j = 1 m γ j F 3 ( ω j , τ ) = ε 1 ( τ ) ( ω 1 ω m ω 1 2 0 ε 2 ( ω , u ) ( 1 + u ) - β x ( u ) d u d ω
- j = 1 m γ j ω j 1 2 0 ε 2 ( ω j , u ) ( 1 + u ) - β x ( u ) d u ) ,
ω 1 ω m F 4 ( ω , τ ) 𝑑 ω - j = 1 m γ j F 4 ( ω j , τ ) = ω 1 ω m ε 2 ( ω , τ ) 0 ε 2 ( ω , u ) ( 1 + u ) - β x ( u ) 𝑑 u 𝑑 ω
- j = 1 m γ j ε 2 ( ω j , τ ) 0 ε 2 ( ω j , u ) ( 1 + u ) - β x ( u ) 𝑑 u .

Let us estimate all the above components. Since the quadrature formula (2.15) is exact for polynomials of the first degree, we have

ω 1 ω m F 1 ( ω , τ ) 𝑑 ω - j = 1 m γ j F 1 ( ω j , τ ) = 0 .

With the help of relation (2.16) and conditions (a)–(b), we find

(2.20) ω 1 ω m F 2 ( ω , τ ) d ω - j = 1 m γ j F 2 ( ω j , τ ) c r , Ω h r max ω | r ( ω 1 2 ε 2 ( ω , τ ) ) ω r | ( 1 + ) - β ε 1 L 2 , 2 x L 2 , 2 .

Further, taking into account relation (2.16) and conditions (a)–(b), we obtain

ω 1 ω m F 3 ( ω , τ ) 𝑑 ω - j = 1 m γ j F 3 ( ω j , τ ) c r , Ω h r | ε 1 ( τ ) | max ω | r ( ω 1 2 0 ε 2 ( ω , u ) ( 1 + u ) - β x ( u ) 𝑑 u ) ω r |
(2.21) c r , Ω h r | ε 1 ( τ ) | max ω g 1 ( ω ) x L 2 , 2 .

Finally, using relation (2.16) and condition (b), we get

ω 1 ω m F 4 ( ω , τ ) 𝑑 ω - j = 1 m γ j F 4 ( ω j , τ ) c r , Ω h r max ω | r ( ε 2 ( ω , τ ) 0 ε 2 ( ω , u ) ( 1 + u ) - β x ( u ) 𝑑 u ) ω r |
(2.22) c r , Ω h r max ω g 2 ( ω , τ ) x L 2 , 2 .

Substituting (2.20)–(2.22) into formula (2.19), we arrive at the desired estimate. ∎

3 Application to DDT-Problem

In this section we interpret the outcome of our analysis for the context of DDT-problem. To be specific, we consider equation (1.1) with the kernel

(3.1) K ( ω , τ ) = i ω τ tanh i ω τ ,

corresponding to DDT-problem for blocking boundary conditions and planar symmetry [17]. By separating real and imaginary parts, we obtain that

K ( ω , τ ) = K 1 ( ω , τ ) 2 ω τ 2 + i K 2 ( ω , τ ) 2 ω τ 2 ,

with

K 1 ( ω , τ ) = sinh 2 ω τ - sin 2 ω τ cosh 2 ω τ + cos 2 ω τ ,
K 2 ( ω , τ ) = sinh 2 ω τ + sin 2 ω τ cosh 2 ω τ + cos 2 ω τ .

Since we are interested in a real-valued solution p ( τ ) of (1.1), this equation can be reformulated into a system of two integral equations of the form (1.1) with the kernels K = K 1 , K = K 2 and y ( ω ) = Y 1 ( ω ) , y ( ω ) = Y 2 ( ω ) , where Y 1 , Y 2 are real and imaginary parts of Y ( ω ) = Y 1 ( ω ) + i Y 2 ( ω ) . Note that a similar situation takes place in DRT-problem, where one often uses only a part of the impedance data Y 1 or Y 2 to reconstruct the quantity of interest (see, e.g., [3, 2, 23]). Here, in the same spirit, we consider an equation (1.1) with K = K 1 and y = Y 1 corresponding to the real part of DDT-problem. Moreover, it is convenient for us to change the variables as τ τ , such that we arrive at an equation with operator H of the form (2.2)–(2.3) with

(3.2) K ( ω , τ ) = 2 ω K 1 ( ω , τ 2 ) τ 2 , x ( τ ) = ( 1 + τ ) β p ( τ 2 ) .

Taking into account the form of K 1 ( ω , τ 2 ) , we have the representation

K ( ω , τ ) = sinh τ 2 ω - sin τ 2 ω cosh τ 2 ω + cos τ 2 ω 2 ω τ 2
= 2 ω τ 2 ( 1 - [ 1 - sinh τ 2 ω - sin τ 2 ω cosh τ 2 ω + cos τ 2 ω ] )
= 2 ω τ 2 ( 1 - cosh τ 2 ω - sinh τ 2 ω + cos τ 2 ω + sin τ 2 ω cosh τ 2 ω + cos τ 2 ω )
= 2 ω τ 2 ( 1 - e - τ 2 ω + 2 sin ( π 4 + τ 2 ω ) 1 + e - 2 τ 2 ω + 2 e - τ 2 ω cos τ 2 ω 2 e - τ 2 ω ) ,

telling us that K ( ω , τ ) has the form (2.17), where

ε 1 ( τ ) = 2 τ 2 ,
ε 2 ( ω , τ ) = 2 ω τ 2 e - τ 2 ω + 2 sin ( π 4 + τ 2 ω ) 1 + e - 2 τ 2 ω + 2 e - τ 2 ω cos τ 2 ω 2 e - τ 2 ω .

So, ε 1 ( τ ) increases with τ as O ( τ 2 ) . Therefore, the multiplier p ( τ 2 ) of the unknown x ( τ ) should belong to L 2 β , β > 5 2 to compensate such an increase. Therefore condition (a) is satisfied with β > 5 2 . Conditions (b) can easily be verified directly.

To find a regularized solution of the integral equation (2.2), (2.3) with the above determined kernel we have the following system of linear algebraic equations, due to (2.7):

(3.3) α x α j + γ j l = 1 m a j l x α l = γ j y ( ω j ) , j = 1 , , m ,

with

a j l = ω j ω l 0 [ 1 - Q ( ω j , v ) ] [ 1 - Q ( ω l , v ) ] v 4 ( 1 + v ) - 2 β 𝑑 v ,

where

Q ( ω j , v ) = e - v 2 ω + 2 sin ( π 4 + v 2 ω ) 1 + e - 2 v 2 ω + 2 e - v 2 ω cos v 2 ω 2 e - v 2 ω .

For the above integrand we have

[ 1 - Q ( ω j , v ) ] [ 1 - Q ( ω l , v ) ] v 4 ( 1 + v ) - 2 β = q c ( v ) + q j ( v ) + q l ( v ) + q j l ( v ) ,

where

q c ( v ) = v 4 ( 1 + v ) 2 β , q j ( v ) = v 4 ( 1 + v ) 2 β Q ( ω j , v ) , q j l ( v ) = v 4 ( 1 + v ) 2 β Q ( ω j , v ) Q ( ω l , v ) , j , l = 1 , , m .

It is clear that

(3.4) I 1 = 0 q c ( v ) 𝑑 v = 24 ( 2 β - 5 ) ( 2 β - 4 ) ( 2 β - 3 ) ( 2 β - 2 ) ( 2 β - 1 ) .

To estimate the integrals of other components q j , q j l , we recall [18] that if f ( x ) is such that

(3.5) | f ( z ) | C ( | z | 1 + | z | ) α e - β Re ( z ) , 0 < α 1 ,  0 < β 1 ,

for

z { z : | arg ( sinh z ) | < d } , 0 < d < π ,

then (see e.g. [18, Example 4.2.11, p. 194 and Corollary 4.2.7, pp. 188–189])

(3.6) I ( f ) = 0 f ( x ) 𝑑 x I N ( f ) = h k = - N N f ( x k ) 1 + e - 2 k h , x k = ln ( e k h + 1 + e 2 k h ) , h = 1 N ,

and there exist positive constants C, γ such that

| I ( f ) - I N ( f ) | C exp ( - γ N ) .

One can check that the functions q j ( x ) , q j l ( x ) satisfy the above condition (3.5). For example, for real values z > 0 we have

| q j ( z ) | 2 | e - z 2 ω j + 2 sin ( π 4 + z 2 ω j ) | | 1 + e - 2 z 2 ω j + 2 e - z 2 ω j cos z 2 ω j | e - z 2 ω j z 4 ( 1 + z ) 2 β
2 ( 1 + 2 ) 0.87 e - z 2 ω j z 4 ( 1 + z ) 2 β 5.55 e - z 2 ω j z 4 ( 1 + z ) 2 β ,

that corresponds to (3.5) with α = β = 1 . Similar bounds are valid for q j l ( v ) .

Thus, using (3.4) and (3.6), we estimate the coefficients of system (3.3) as follows:

a j l I 1 + I N ( q j ) + I N ( q l ) + I N ( q j l ) .

After solving system (3.3), the regularized solution of the integral equation (2.2), (2.3) is defined by (2.6). Remark that system (3.3) possesses a symmetry matrix. Therefore it can be solved using one of the efficient method.

The next step is to select a value of the regularization parameter α to provide an accurate approximation for the solution of the problem (2.1)–(2.3), (3.2). The optimal value of α is defined in Theorem 2, but it requires a knowledge of the smoothness index function φ and the perturbation level δ. In the considered application such information is usually not available. To resolve this issue we propose to use the well-known quasi-optimality criterion [7, 20, 6, 9]. Recall that in this criterion, we consider a geometric sequence of values of the regularization parameters α n = α 0 q n , n = 0 , 1 , , N , with a fixed q < 1 and α 0 > 0 . Then for each α = α n we compute the corresponding regularized approximation x α n , m δ given by (2.6), (3.3), and select the value α = α n * with

n * = arg min { x α n , m δ - x α n - 1 , m δ L 2 , : n = 1 , 2 , , N } .

In the numerical tests below we use the quasi-optimality criterion with α 0 = 1 , q = 0.5 , N = 50 , and we stop to compute x α n , m δ at n = k if

x α k , m δ - x α k - 1 , m δ L 2 , > x α k - 1 , m δ - x α k - 2 , m δ L 2 , .

4 Numerical Experiments

We consider a DDT-problem formulated as the integral equation (1.1) with the kernel (3.1). First we simulate EIS data y ( ω ) in such a way that the exact solution of (1.1) is given as p ( τ ) = exp ( - ln 2 ( τ ) ) τ (see Figure 1). Note that such p ( τ ) mimics some of DDT-functions considered in [17].

Figure 1 
               Synthetic DDT-function 
                     
                        
                           
                              
                                 p
                                 ⁢
                                 
                                    (
                                    τ
                                    )
                                 
                              
                              =
                              
                                 
                                    exp
                                    ⁡
                                    
                                       (
                                       
                                          -
                                          
                                             
                                                ln
                                                2
                                             
                                             ⁡
                                             
                                                (
                                                τ
                                                )
                                             
                                          
                                       
                                       )
                                    
                                 
                                 τ
                              
                           
                        
                        
                        {p(\tau)=\frac{\exp(-\ln^{2}(\tau))}{\tau}}
                     
                  .
Figure 1

Synthetic DDT-function p ( τ ) = exp ( - ln 2 ( τ ) ) τ .

Figure 2 
               The error of the reconstruction from non-perturbed data:

                     
                        
                           
                              
                                 
                                    ∥
                                    
                                       
                                          p
                                          ⁢
                                          
                                             (
                                             τ
                                             )
                                          
                                       
                                       -
                                       
                                          
                                             p
                                             α
                                          
                                          ⁢
                                          
                                             (
                                             τ
                                             )
                                          
                                       
                                    
                                    ∥
                                 
                                 
                                    L
                                    
                                       2
                                       ,
                                       ∞
                                    
                                 
                              
                              =
                              
                                 0.581
                                 *
                                 
                                    10
                                    
                                       -
                                       4
                                    
                                 
                              
                           
                        
                        
                        {\|p(\tau)-p_{\alpha}(\tau)\|_{L_{2,\infty}}=0.581*10^{-4}}
                     
                  .
Figure 2

The error of the reconstruction from non-perturbed data: p ( τ ) - p α ( τ ) L 2 , = 0.581 * 10 - 4 .

EIS data sampling is simulated in the range Ω = [ ω 1 , ω m ] with ω 1 = 10 - 2 , ω m = 10 2 . Moreover, to be close to real EIS data set we suppose that sampling points ω i are uniformly distributed in four decades Ω j = [ 10 j - 2 , 10 j - 1 ] , j = 0 , 1 , 2 , 3 , with five points in each decade. To illustrate the algorithm described in the previous section, we apply it to the non-perturbed simulated data ( δ = 0 , Example 1), and to noisy data y δ ( ω i ) = y ( ω i ) + δ i , i = 1 , 2 , , m , where δ i are the realizations of zero-mean normally distributed random variables corresponding to noise-to-signal ratio of 1 % (Example 2) and 10 % (Example 3) respectively. The error observed in Example 1 is displayed in Figure 2. The reconstructions obtained for Examples 2 and 3 are displayed in Figures 3 and 4. The relative errors in L 2 , -norm observed in these examples are 4.3 % and 23 % , respectively. Keeping in mind the corresponding noise-to-signal ratios, one can conclude that in the considered examples the proposed algorithm demonstrates a stable and reliable behavior.

Figure 3 
               Example 2. Reconstruction from noisy data.
Figure 3

Example 2. Reconstruction from noisy data.

Figure 4 
               Example 3. Reconstruction from noisy data.
Figure 4

Example 3. Reconstruction from noisy data.

We perform one more test, where EIS data y ( ω ) correspond to the exact solution of (1.1) given as p ( τ ) = exp ( - ln ( τ ) 2 ) + 1.3 exp ( - 2 ( 2 - ln ( τ ) ) 2 ) (see Figure 5). This solution mimics a bimodal behavior of DDT-functions. Data simulations are performed in the same way as above. The algorithm is applied to non-perturbed data (Example 4) and to noisy data corresponding to noise-to-signal ratio of 2 % (Example 5). The corresponding reconstructions are displayed in Figures 6 and 7.

Figure 5 
               Bimodal synthetic DDT-function.
Figure 5

Bimodal synthetic DDT-function.

Figure 6 
               Example 4. Reconstruction from non-perturbed data.
Figure 6

Example 4. Reconstruction from non-perturbed data.

Figure 7 
               Example 5. Reconstruction from noisy data.
Figure 7

Example 5. Reconstruction from noisy data.

Note that even in the case of non-perturbed data some reconstruction error unavoidably appears due to data discretization. In the considered case the discretization causes a relative error of 0.16 % . The reconstruction from noisy data is given with a relative error of 5.2 % . The results of this test demonstrate that the proposed algorithm is able to reconstruct important geometric features of DDT-functions such as the number of modes (picks) and their locations.

Remark 2.

In the framework of this article, we restrict our consideration to solutions (2.9) of low smoothness, i.e. φ ( λ ) = λ ν for ν 1 , because of the known Tikhonov method’s saturation effect in the case ν > 1 . This may lead to a loss of the optimal order of accuracy. It is necessary to use a regularization method of higher qualification than that of the Tikhonov method to expand the set of problems to be solved in the case of smoother solutions without loss of accuracy. This problem will be investigated in the future works.

Remark 3.

In Theorem 2 the order-optimal accuracy estimate was established in the case of a priori choice of the regularization parameter λ. In Section 4 the calculations were performed for a posteriori rule of choosing λ, namely, the quasi-optimality criterion, which does not require knowledge of the quantity δ. In turn, this leads to the fact that with such a choice of λ it is impossible to estimate the approximation error. To find the (order-optimal) error estimate in the a posteriori case, it is necessary to apply another stopping rule, for example, the discrepancy principle or the balancing principle. The authors plan to establish similar results in their future investigations.

5 Conclusions

In this article we have analyzed a regularization scheme for solving Fredholm integral equations of the first kind with possibly non-square-integrable kernels. Such equations arise in EIS data processing performed in terms of DDT-functions. In our analysis we use a weighted space setting that allows an application of methods from the Regularization theory. The analyzed scheme is illustrated in application to DDT-problem for blocking boundary conditions and planar symmetry. In numerical implementation of the considered scheme we use fast (exponentially convergent) quadrature rules to find coefficients of a system of linear algebraic equations determining regularized solutions. Moreover, we use quasi-optimality criterion for selecting a value of the regularization parameter. Numerical experiments with synthetic data demonstrate reliability of the proposed scheme.

Note, however, that the exact solution of DDT-problem (1.1) is expected to be non-negative, while the regularized approximation (2.8) may not possess such a property (see, e.g., Figure 3). Of course, a negative part of the regularized solution can be attributed to the influence of noise and therefore ignored. Another way to deal with this issue is to incorporate non-negativity constraint into Tikhonov regularization that can be a subject of future research.

Funding statement: This research was done when the second and the third co-authors performed their secondments at the Johann Radon Institute for Computational and Applied Mathematics (RICAM) of the Austrian Academy of Sciences within the framework of EU-Horizon2020-MSCA-RISE-project AMMODIT, while the forth co-author was hosted at RICAM within the framework of Joint Excellence in Science& Hummanities (JESH) programme of the Austrian Academy of Sciences. Inspiring working conditions of RICAM, as well as, the support of AMMODIT-consortium and JESH-programme are gratefully acknowledged.

References

[1] E. Barsoukov and J. R. Macdonald, Impedance Spectroscopy: Theory, Experiment, and Applications, John Wiley & Sons, New York, 2005. 10.1002/0471716243Search in Google Scholar

[2] B. A. Boukamp and A. Rolle, Analysis and application of distribution of relaxation times in solid state ionics, Solid State Ionics 302 (2017), 12–18. 10.1016/j.ssi.2016.10.009Search in Google Scholar

[3] F. Dion and A. Lasia, The use of regularization methods in the deconvolution of underlying distributions in electrochemical processes, J. Electr. Chem. 475 (1999), 28–37. 10.1016/S0022-0728(99)00334-4Search in Google Scholar

[4] A. L. Gavrilyuk, D. A. Osinkin and D. I. Bronin, The use of Tikhonov regularization method for calculating the distribution function of relaxation times in impedance spectroscopy, Russian J. Electr. 53 (2017), 575–588. 10.1134/S1023193517060040Search in Google Scholar

[5] J. K. Hansen, J. D. Hogue, G. K. Sander, R. A. Renaut and S. C. Popat, Non-negatively constrained least squares and parameter choice by the residual periodogram for the inversion of electrochemical impedance spectroscopy data, J. Comput. Appl. Math. 278 (2015), 52–74. 10.1016/j.cam.2014.09.017Search in Google Scholar

[6] S. Kindermann and A. Neubauer, On the convergence of the quasioptimality criterion for (iterated) Tikhonov regularization, Inverse Probl. Imaging 2 (2008), no. 2, 291–299. 10.3934/ipi.2008.2.291Search in Google Scholar

[7] S. Kindermann, S. Pereverzyev, Jr. and A. Pilipenko, The quasi-optimality criterion in the linear functional strategy, Inverse Problems 34 (2018), no. 7, Article ID 075001. 10.1088/1361-6420/aabe4fSearch in Google Scholar

[8] D. Levin, Stable integration rules with scattered integration points, J. Comput. Appl. Math. 112 (1999), no. 1–2, 181–187, 10.1016/S0377-0427(99)00218-6Search in Google Scholar

[9] S. Lu and P. Mathé, Heuristic parameter selection based on functional minimization: optimality and model function approach, Math. Comp. 82 (2013), no. 283, 1609–1630. 10.1090/S0025-5718-2013-02674-9Search in Google Scholar

[10] S. Lu and S. V. Pereverzev, Regularization Theory for Ill-posed Problems. Selected Topics, Inverse Ill-posed Probl. Ser. 58, De Gruyter, Berlin, 2013. 10.1515/9783110286496Search in Google Scholar

[11] P. Mathé and B. Hofmann, How general are general source conditions?, Inverse Problems 24 (2008), no. 1, Article ID 015009. 10.1088/0266-5611/24/1/015009Search in Google Scholar

[12] P. Mathé and S. V. Pereverzev, Moduli of continuity for operator valued functions, Numer. Funct. Anal. Optim. 23 (2002), no. 5–6, 623–631. 10.1081/NFA-120014755Search in Google Scholar

[13] P. Mathé and S. V. Pereverzev, Geometry of linear ill-posed problems in variable Hilbert scales, Inverse Problems 19 (2003), no. 3, 789–803. 10.1088/0266-5611/19/3/319Search in Google Scholar

[14] M. T. Nair and S. V. Pereverzev, Regularized collocation method for Fredholm integral equations of the first kind, J. Complexity 23 (2007), no. 4–6, 454–467. 10.1016/j.jco.2006.09.002Search in Google Scholar

[15] R. A. Renaut, R. Baker, M. Horst, C. Johnson and D. Nasir, Stability and error analysis of the polarization estimation inverse problem for microbial fuel cells, Inverse Problems 29 (2013), no. 4, Article ID 045006. 10.1088/0266-5611/29/4/045006Search in Google Scholar

[16] M. Saccoccio, T. H. Wan, C. Chen and F. Ciucci, Optimal regularization in distribution of relaxation times applied to electrochemical impedance spectroscopy: Ridge and Lasso regression methods - A theoretical and experimental study, Electr. Acta 147 (2014), 470–482. 10.1016/j.electacta.2014.09.058Search in Google Scholar

[17] J. Song and M. Z. Bazant, Electrochemical impedance imaging via the distribution of diffusion times, Phys. Rev. Lett. 120 (2018), no. 11, Article ID 116001. 10.1103/PhysRevLett.120.116001Search in Google Scholar

[18] F. Stenger, Numerical Methods Based on Sinc and Analytic Functions, Springer Ser. Comput. Math. 20, Springer, New York, 1993. 10.1007/978-1-4612-2706-9Search in Google Scholar

[19] A. N. Tikhonov and V. Y. Arsenin, Solutions of Ill-posed Problems, John Wiley & Sons, New York, 1977. Search in Google Scholar

[20] A. Tikhonov and V. Glasko, Use of the regularization method in non-linear problems, USSR Comput. Math. Math. Phys. 5 (1965), no. 3, 93–107. 10.1016/0041-5553(65)90150-3Search in Google Scholar

[21] G. M. Vaĭnikko and A. Y. Veretennikov, Iteration Procedures in Ill-Posed Problems (in Russian), “Nauka”, Moscow, 1986. Search in Google Scholar

[22] T. H. Wan, M. Saccoccio, C. Chen and F. Ciucci, Influence of the discretization methods on the distribution of relaxation times derconvolution: Implementing radial basis functions with DRT tools, Electr. Acta 184 (2015), 483–499. 10.1016/j.electacta.2015.09.097Search in Google Scholar

[23] Y. Zhang, Y. Chen, M. Li, M. Yan, M. Ni and C. Xia, A high-precision approach to reconstruct distribution of relaxation times from electrochemical impedance spectroscopy, J. Power Sources 308 (2016), 1–6. 10.1016/j.jpowsour.2016.01.067Search in Google Scholar

[24] M. Zic and S. Pereverzyev, Jr., Adaptive multi-parameter regularization in electrochemical impedance spectroscopy, RICAM Report 2018-16, Austrian Academy of Sciences, 2018. Search in Google Scholar

Received: 2019-07-18
Revised: 2019-10-11
Accepted: 2019-10-16
Published Online: 2019-11-08
Published in Print: 2020-07-01

© 2020 Walter de Gruyter GmbH, Berlin/Boston

Downloaded on 24.4.2024 from https://www.degruyter.com/document/doi/10.1515/cmam-2019-0111/html
Scroll to top button