Comment on “a spline-based method for experimental data deconvolution” by I. Beniaminy and M. Deutsch

5

Computer Physics Communications 24 (1981) 5—7 North-Holland Publishing Company

COMMENT ON “A SPLINE-BASED METHOD FOR EXPERIMENTAL DATA DECONVOLUTION” BY I. BENIAMINY AND M. DEUTSCH Henk FREDRIKZE and Peter VERKERK Interuniversitafr Reactor Instituut, 2629 JB Deift, The Netherlands Received 7 July 1981

Beniaminy and Deutsch [1] present an algorithm and a computer program to solve the Fredholm integral equation of the first kind with translational invariant kernel, which may be written as J(t)

=

f

K(t

—

u)Jo(u) du.

(1)

J(t) could be for instance the result of a measurement of J0 (t) with a resolution function K(t). The algo-

rithm is based on the following argument: Let J0 (1) be represented by a piece-wise cubic spline n—i J 0 (t)

=

E

and coefficients Ak’ Bk, Ck and Dk, that are related to the coefficients ak, bk, Ck and dk of J0 (t) by Ak Bk

=

Ck

akMo bkMo

—

CkMO

—

3akMl 2bkMi + 3akM2

—

ckMi + bkM2 — akM3,

Dk = dkMo

(6)

where Mm are the moments of the resolution function mK(t) dt. (7) Mm = t

f

(2)

Sk(t)

However, because substitution of (2) and (3) into (1) yields the equation

k= 1

with 10,

t<~k,

2 +Ckt+dk,

sk(t)Jakt3 +bkt

10,

E j’tk+1 (aku3 +bku2 +CkU +dk)K(t_u)du, (8)

~

~k~t~<~k+i,

J(t)

t > ~+i.

which is not a cubic spline m general, this argument is

k1

The coefficients in (3) are such that J 0(t) is continuous and has continuous first and second derivatives at the knots ~k• Beniaminy and Deutsch fmd that by substitution of (2)spine and (3) into J(t) knots: is also represented by a cubic with the(1) same

n1 1(t) =

~I)Sk(t)

(4)

k1

with 10, Sk(t)Akt~+Bkt+Ckt+Dk,

10,

t
~k~t~k+1,

t>

false. Moreover, it is clear from (8) that for t E [Ek, ~k+1] J(t) depends not only on the coefficients ak, bk, Cktoo, andifdkt lies but close at least on ak_i, Ck_<“width’1 1 and 4_i enough to ~kbk_i, (t ~k of K(t)). This observation contradicts eq. (6). The function 1(t) defined by (4) through (7), from now on denoted by Rt) to differentiate from 1(t) defmed by (8), and its derivatives are in general discontinuous at the knots ~k of J0 (t). Define d’1 d’~ ~lim —~Y(t±e)] (9) —

,

(5)

~k+1’

0010-4655/81/0000—0000/$02.75 © 1981 North-Holland

d~J~t~

Lsiod

then it is easily shown, using the fact that J~(t) and its

H. Fredrikze, P. Verkerk / Comment on a spline-based method

6

f 10(1

first and second derivativesare continuous at ~k’ that Y+(tk)—L(t~k)—M

3(ak

dt

7÷ (t = ~k) =

—6M1 (ak

3M2(ak—ak_1), (t~k)

—

— -~--~

7

(t =

=

(10)

~k)

ak_i).

—

From (10) we conclude that in the special case of a symmetrical resolution function 7(t) and its second derivative are continuous at t = ~k’ whereas its first derivative is discontinuous. Although the theoretical basis of the algorithm of Beniaminy and Deutsch is incorrect, the application of their deconvolution program to the example shown in ref. [1] gives surprisingly good results. The exampie used is characterized by 2, (11) 10(t) = (sin 2irt/2irt) I

lOiti),

ti >0.1.

tamed from numerical integration of (1) and those For this example compared the values for 1(t) obfor 7(t), using for we 1 0(t) the cubic spline interpolation according to the algorithm of Reinsch [2] for equidistant knots (~= —2.5, ~ = 0.05). The differences were found to be of the order iO~,which explains the results of the numerical test in ref. [1]. However, severe differences were observed using a five times broader resolution function (2(1 K(t) = ~o,

—

21t1),

ti ~ 0.5 t> 0.5,

(13)

keeping the other parameters as above. In fig. 1 we show the spline representation of 10(t), the result for 1(t) from numerical integration, and the result 7(t) as calculated from (4) through (7). In this case 7(t) is continuous at t = ~k as expected, the discontinuities in Y’(t) at t = ~k can be observed clearly. I

A

I’

0.75

I I I:

I

1.00

I

A

I’

I ....

I

-

II

-

i I

0.50

0.00

(12)

I

1.00~

0.25

ti ~ 0.1,

K(t)=

ak_i),

—

dd2 d2 7+ (t = ~k )——J

—

—

0.75

0.50

ii

-

\.~

•1

~

-ao

I

Q00 0.25 I

-I.0

I

I

0.0

—

—

:1

I

.0

20

I

-2.0 Fig. 1. The spline representations of Jo(t) (dashed line), J(r) (dotted line), and 1(t) (solid line). J(t) is obtained by numerical integration of (1) for J 0(t) given by (11) and for the resolution function given by (13). The latter is indicated by the triangle. .1(t) is given by (4)—(7).

I

I

-1.0

0.0

I

1.0

t

20

Fig. 2. The spline representations of JQ(t) (dashed line), .1(t) (dotted line), both as in fig. 1, and of J~(t)(solid line) obtallied by (2), (3), (6) and (7). The resolution function (13) is also indicated.

H. Fredrikze, P. Verkerk / Comment on a spline-based method

From the above we expect that the deconvolution algorithm of Beniaminy and Deutsch applied to 1(t) obtained by numerical integration is bound to give results that deviate from J0 (t). The inversion of (6), which is essentially the deconvolution method of Beniaminy and Deutsch, applied to the Reinsch spine interpolation [2] (~ —2.0, = 0.05) of the 7o= (t). Wez~ show this func. latter 1(t), gives a fur ction tion 7 0(t) together with 10(t) and 1(t) in fig. 2. indeed a large deviation of 70(t) from J0(t) is observed. From these examples we conclude that great care

7

is to be taken if the program DECONV of Beniaminy and Deutsch, based on the algorithm given by (2) through (7), is used. The deconvolution program developed recently by one of us [3], applied to J(t) of the latter example gave satisfactory results. References [1] I.21(1980) Beniaminy and M. Deutsch, Comput. Phys. Commun. 271. [21C.H. Reinsch, Numer. Math. 10 (1967) 177. [3] P. Verkerk, Comput. Phys. Commun., to be submitted.