A numerical solution of one class of Volterra integral equations of the first kind in terms of the machine arithmetic features

The paper is devoted to the studies of a special class of Volterra integral equations of the first kind

t

У K n ( t - s ) Ф ( s ) d s = y ( t ) ,                              (1)

N

K n ( t - s ) = ^ ( - 1) ^{q +1} q ² e- ^{-n 2 q 2(} t ^s ) .                           (2)

q =1

The specific feature of kernel K N ( t — s ) in (2) lies in that K N = 0 in some neighborhood of zero. The qualitative theory and numerical methods of solving the Volterra equations of the first kind are dealt with in many studies (for example, [1-4] and the eqreferences given in them).

The goal of this paper is to consider the application of numerical methods for solving the equations of form (1), (2) taking into account the mechanisms of error occurrence in the computer calculations. The work is a continuation of the research started in [5,6].

The Volterra equation of the convolution type (1), (2) was first obtained in [7]. The authors of [7] suggest a method of searching for a solution u (1 ,t ) = ф ( t ), t >  0, of an inverse boundary-value problem

U t = U xx , X e (0 , 1) , t > 0 ,                               (3)

u (x, 0) = 0, u (0, t) = 0, Ux (0, t) = gо (t)

by reducing (3), (4) to the Volterra integral equation of convolution type:

t

У K(t - s)Ф(s)ds = y(t), 0 < s < t < T,(5)

о

∞

K(t - s) = £ (-1)q+1 q2e-2q2(t-s), y(t) = 5-5g0 (t).(в)

2 n ²

q =1

Instead of g ₀( t ) we norm;ally know g _s ( t ): || g _s ( t ) — g ₀( t ) || _c 6 5. 5 >  0 .

Problem (3), (4) plays an important part in the applied problems, including those related to the research of non-stationary thermal processes. Solving of the inverse problems is, as a rule, complicated by the instability of these problems with respect to the initial data errors. The application of the methods which employ the Fourier and Laplace transforms in combination with the theory of ill-posed problems found wide application in the construction of stable solutions to the inverse problems of heat conductance. In particular, to solve the problem similar to (3), (4), the authors of [8] used a stabilizing functional after applying the Fourier transform. In [9] to regularize and assess the convergence of the obtained solutions the authors used a method of conjugate gradients. In [10] and [11] the Laplace transform was considered for solving the Cauchy problem. In [12] the Laplace transform was used in a two-dimensional problem. The existing approaches, as a rule, after taking the Laplace transform, apply regularization methods to the obtained equations and then perform the inverse transform.

Taking into account the ideas from [13,14] the authors of [7] approximated u x (0 ,t ) by the sum of N first summands:

U x (0 ,t ) = 2 n ²

£ ( — 1) ’ ⁺¹ q ² / e ^{- 2 q 2( t-s} ’ Ф ( s ) ds, q -¹             0

where N is a positive integer. Then (5), (6) are reduced to the form (1), (2). The performance of this approach was discussed in [7]. According to the conclusion made by the author, such a way of solving of the inverse problem makes it possible to reduce the initial problem to the Volterra integral equation of the first kind and exclude the components of the operator calculus from the regularization process.

It is known that the Volterra integral equations of the first kind belong to the class of conditionally-correct problems, and the discretization procedure has a regularizing feature with a regularization parameter, a step of mesh, which is in a certain manner connected to the level of disturbances of the initial data 5.

This paper presents an algorithm for numerical solution of (1), (2) at an exactly specified right-hand part. Note that when solving (1), (2) we face three types of errors related first of all to the approximation of the initial problem (5), (6), secondly to the accuracy of a numerical method, and finally to the computation errors in the machine arithmetic operations with real floating-point numbers. For the research, of greatest interest is the first of the indicated cases. However, to pass to the problem of assessing the parameter N in (2) it is necessary to develop an algorithm, for computation of K N max. which takes into account the specific features of machine arithmetic and provides the desired (specified) number of valid digits in the significand.

1.    The Specific Features of the Numerical Calculations

The computational experiment in [5] shows that with an increase in the number of summands in (2) the roots X* of the equation K N ( X ) = 0 decrease (at the same time monotonisity is observed only separately in even and odd N ). The values X* will be used further to limit the magnitude of the mesh step h from above, for the value of the mesh function K N at the first node to be non-zero. As is known, the condition K N (0) = 0 is necessary for (1), (2) to be correct on the pair ( C,C [1_{0 T} ] ), where y (0) = 0, y‘ ( t ) G C [0 , _T ].

Note, that there can be computational errors in the calculation of X* , they can be related to the application of a fixed mesh in the machine number representation.

Let us consider the known cases of systematic error accumulation [15] which appear in the calculation of kernel values (2). Use the system of computer algebra MaplelO. Following [16] we will include parameter f in the generally accepted representation of the real number. The parameter is equal to the number of valid digits in the significand (starting from the left). Assume that the real number x = s • M • 10 - L + p is specified by the set ( s, M,p, f ). where s G { — 1 , 0 , +1 } is the sign of the number. M G { 10 L - ¹ , 10 L - ¹ + 1 ,..., 10 L — 1 } U { 0 } is significand cT the number. L is number of significand positions, p is exponent part of the number.

Let us illustrate the details of the calculations when summing up the numbers with different exponent parts in (2), using the example.

Example 1. Let N = 50. X о = 10 - 3. L >  8. Take

34                                    50

x i = E ( — 1) q ⁺¹ q ² e- n ² q ^{2 x} 0 , x 2 = E ( — 1) q ⁺¹ q ² e- ² q ^{2 x} 0 q =11                              q =35

and find x 3 = x ₁ + x2.

Assuming

10p s-fs = 10p 1 -f 1 + 10p 2 -f2, p 3 > p 1 > p 2, according to [16], it is easy to obtain that f3 > fs = [f 1 — lg(1 + 10-p 1+f 1+p2-f2)],                          (7)

where the symbol [ ... ] means the greatest integer. Table 1 presents the parameters (1 ,M ₁ , 2 ,f 1). (1 ,M 2 , — 2 ,f 2) and (1 ,M ₃ , 2 ,f ₃). which define x ₁. x 2 anid x ₃ respectively. The last column of the table shows the values of minorant fs, that are calculated using (7).

The next example illustrates the situation arising in the calculation of the difference between the numbers which have coinciding exponent parts and several high-order digits of the significand.

Example 2. Let N = 50. X 0 = 10 3. L >  8. Introduce 10                                     50

x з = E ( — 1) q ⁺¹ q ² e- n ² q ^{2 x} 0 , x 4 = E ( — 1) q ⁺¹ q ² e- ² q ^{2 x} 0 . q =1                               q =11

Define x д = |x 4 1 — |x 3 1.

Table 1

L	M 1	f 1	M 2	f 2	M Σ	f Σ	f s
8	18652239	6	44981421	6	18656737	6	5
9	186522441	8	449814458	7	186567422	8	7
10	1865224455	9	4498144699	8	1865674269	8	8
11	18652244592	11	44981446726	8	18656742737	11	10
12	186522445926	11	449814466957	10	186567427373	11	10
13	1865224459248	12	4498144669376	10	1865674273715	12	11
14	18652244592468	13	44981446694089	11	18656742737137	13	12

A allies M and f for x 1. x 2 aiid x 3

Suppose that

|M 4 — M 3 1 <  10 L ¹ — 1 , p д < p 3 = p 4

and, following [16] use the empiric estimate:

f Д > f r = (

0 , if f min — L + lg( A ) <  0 ,

[fmin — L + lg(A)] , if fmin — L + lg(A) > 0, where A = |M4 — M31 + 1, fmin = min{f3, f4}.

Below are the parameters ( — 1 ,M ₃ , 2 ,f ₃), (1 ,M ₄ , 2 ,f ₄) and ( — 1 ,M _д , 2 ,f _д), which specify the values x ₃. x ₄ arid x _д. The estimation1 from below of fr. obtained using (8). is given in the last column of Table 2.

A allies M and f for x ₃. x ₄ arid x _д

Table 2

L	M 3	f 3	M 4	f 4	M ∆	f ∆	f r
8	18656743	7	18656737	6	00000006	0	0
9	186567428	8	186567424	8	000000004	0	0
10	1865674274	9	1865674268	8	0000000006	0	0
11	18656742750	11	18656742736	10	00000000014	1	0
12	186567427505	12	186567427372	11	000000000133	3	1
13	1865674275054	12	1865674273715	12	0000000001339	4	2
14	18656742750534	14	18656742737138	13	00000000013396	4	3

It is obvious that when several high-order digits are set to zero, there appears the number with lower quantity of significant digits in the significand. Figure 1 illustrates an instantaneous loss of high-order digits, which takes place in this situation. The calculations are made using the software developed by the authors in PASCAL.

The next section is devoted to the problem of approximate solving of (1), (2) in terms of the specific features of machine arithmetic.

Fig. 1. Instantaneous loss of high-order digits for K N (0 . 001) , where N = 12

2.    Results of Solving of (1), (2)

Let us introduce the uniform meshes of nodes t i = ih, t i - 1 = ( i — 2 ) h, i = 1 ,n,nh = T in [0 ,T ] and, by approximating the integral in (1) by the middle rectangle quadrature and product integration method [17], write the corresponding mesh analogues to (1), (2) as

N

i

h £ ( — 1) q ⁺¹ q ²E

q =1

j = ¹

e-n2q2(j-2)h фЬ    x = yh i-J + 2      i

and

N

E ( — 1) q ⁺¹ q ² q =1

jh

^∑ ϕ j ^h - 1 2

j = ¹         ( J - 1) h

e- ^{n 2} q ² ( ,h - s ) ds = y h .

Designate their solutions by (j)^h and ф\ respectively. Conduct a numerical experiment with the guaranteed number of valid digits in the significand not less than f 5 = 8.

Example 3. Assume ф ( t ) from [18] as the reference function:

^{ф (} t ) =

-

- ^t e ^α

-

- ¹ e ^α

- t,

where a = 10 ¹ , 10 ²

Set in (9), (10) N = 2 , 5; 10; 15. Table 3, 4 give the results of numerical calculations in the integration interval [0 , 1]. Here, the notation is as follows:

Y 1 = log₂

P h ¹ Il C h

I P2 Ik

Y 2 = log 2

EJk

I P 2 Ik ’

where h 1 = 2 h2, \\e ^h \\ c h is maximum absolute difference between the precise solution and the approximate solution at the nodal points (the approximate solution is obtained using the middle rectangle quadrature); \ e h \ c _h is maximum absolute difference between the accurate solution and the approximate solution at the nodal points (the approximate solution is obtained using the product integration method); symbol * means that the error norm is higher than max 1ф ( t ) |.

Table 3

Errors in the mesh solution for function ф, where a = 10 ¹

h	||ε ˇ || C N h =2	γ 1	||ε|| C N h =2	γ 2	||ε ˇ || C N h =3	γ 1	||ε|| C N h =3	γ 2
1 / 256	0,005001	2,009	0,000499	1,996	0,003815	2,002	0,001171	1,994
1 / 512	0,001242	2,002	0,000125	1,998	0,000952	2,000	0,000294	1,998
1 / 1024	0,000310	2,000	0,000031	1,999	0,000238	2,000	0,000074	1,999
1 / 2048	0,000078	1,981	0,000008	1,999	0,000059	1,989	0,000018	2,001
h	||ε ˇ || C N h =4	γ 1	||ε ˆ || C N h =4	γ 2	||ε ˇ || C N h =5	γ 1	||ε ˆ || C N h =5	γ 2
1 / 256	0,065009	2,107	0,002056	1,987	0,025682	2,010	0,003130	1,973
1 / 512	0,015090	2,025	0,000519	1,996	0,006377	2,002	0,000797	1,993
1 / 1024	0,003707	2,006	0,000129	1,999	0,001591	2,000	0,000200	1,998
1 / 2048	0,000923	1,999	0,000032	2,000	0,000398	1,996	0,000050	2,000
h	||ε ˇ || C N h =10	γ 1	||ε ˆ || C N h =10	γ 2	||ε ˇ || C N h =15	γ 1	||ε ˆ || C N h =15	γ 2
1 / 256	∗	—	0,009531	1,744	∗	—	0,013378	1,394
1 / 512	∗		0,002845	1,922	0,485248	2,149	0,005092	1,719
1 / 1024	0,137360	2,212	0,000751	1,979	0,109395	2,033	0,001547	1,910
1 / 2048	0,029650	2,047	0,000190	1,995	0,026724	2,008	0,000411	1,957

Table 4

Errors in the mesh solution for function ф, where a = 10 ²

h	||ε ˇ || C N h =2	γ 1	||ε|| C N h =2	γ 2	||ε ˇ || C N h =3	γ 1	||ε|| C N h =3	γ 2
1 / 256	0,006113	2,009	0,000402	1,995	0,007855	2,004	0,006159	1,875
1 / 512	0,001518	2,002	0,000101	1,998	0,001958	2,001	0,001679	1,939
1 / 1024	0,000379	1,999	0,000025	1,992	0,000489	2,000	0,000438	1,970
1 / 2048	0,000095	1,985	0,000006	1,927	0,000122	1,998	0,000112	1,985
h	||ε ˇ || C N h =4	γ 1	||ε ˆ || C N h =4	γ 2	||ε ˇ || C N h =5	γ 1	||ε ˆ || C N h =5	γ 2
1 / 256	0,080125	2,117	0,014295	1,859	0,051629	2,022	0,024140	1,840
1 / 512	0,018474	2,027	0,003940	1,934	0,012716	2,006	0,006744	1,928
1 / 1024	0,004532	2,006	0,001031	1,968	0,003167	2,001	0,001772	1,966
1 / 2048	0,001128	2,000	0,000264	1,984	0,000791	2,000	0,000453	1,985
h	||ε ˇ || C N h =10	γ 1	||ε ˆ || C N h =10	γ 2	||ε ˇ || C N h =15	γ 1	||ε ˆ || C N h =15	γ 2
1 / 256	∗	—	0,081670	1,584	∗	—	0,114747	1,214
1 / 512	∗		0,027232	1,849	∗		0,049456	1,637
1 / 1024	0,170835	2,232	0,007556	1,945	0,2227532	2,073	0,015899	1,874
1 / 2048	0,036370	2,052	0,001962	1,978	0,0529318	2,018	0,004338	1,959

Tables show that both finite difference methods have the second order of convergence.

Figures 2, and 3 illustrate the behavior of functions Is h i = | ф_i _i n h 1 1, |ё ^h | = |^ _i — ф^к 1 1

2 i - ₂                    2      i - ₂

on a unified mesh with the step h = 27 at fixed values N = 2, a = 10 - 1, T =1. Plot "Line 1" corresponds to function |s h | , plot "Line 2" corresponds to function |S ^h | and plot "Line

3" corresponds to function S h min | = min {|S ^h |, |S ^h |}.

Fig. 2. Absolute values of errors in mesh solutions at exactly specified initial data

Figure 2 was obtained at functions y ^h accurately specified in (9), (10). Here, the maximum values of errors made up

||ε ˇ || _C ^N h ⁼² = 0 , 0795 , ||ε ˆ || _C ^N h ⁼² = 0 , 0424 .

Fig. 3. Absolute values of errors in mesh solutions at perturbed initial data

Figure 3 is obtained at a saw-tooth perturbance of mesh function yh: y ˜( t i ) = y ( t i ) + ( - 1) ⁱ 10 ^{- 2} , i = 1 , 27 , nh = 1 .

In this case, the maximum values of errors made up

\ H \N = ² = 0 , 0638 , \\S\\N=² = 0 , 0609 , \\e min \\N = ² = 0 , 0582 .

Remark. Step h for the fixed level of the initial data perturbances was chosen using the Fibonacci method with ten trials.

The comparison of p hHcy and p hHcy makes it obvious that the use of the product integration method is more preferable. The computational experiment conducted for the given example at a = 10 - ³ shows the convergence for h <  10 - 3. In this connection, further studies are supposed to use the numerical methods of higher order, in particular the third- and fourth-order Runge-Kutta methods. Further, it is planned to construct stability regions of the considered algorithms by the analogy with [19].

Conclusion

The paper presents the research of the approximate solution to the Volterra integral equation of the first kind of convolution type, which occurs in the inverse boundary-value heat conduction problem, by the second-order finite difference methods. The calculation results obtained using the system Maple 10 are presented. The computational experiment is conducted in terms of the error occurrence mechanism in computer calculations. Typical cases of systematic error accumulation are demonstrated by the test examples. Software for the calculation of kernels, tracking the valid digits in the significand, is developed in PASCAL.

Acknowledgements. The work was supported by the RFBR, project No. 15-01-0Ц25-а.