Introduction

Integral equations occur naturally in many areas of mathematical physics. Many engineering and applied science problems arising in water waves, potential theory and electrostatics are reduced to solving integral equations. The problem of finding out the crack energy and distribution of stress in the vicinity of a cruciform crack leads to the integral equation

$$\begin{aligned} \chi (x)+\int _{0}^{1} {\mathcal {L}}(x, t)\chi (t)\mathrm{d}t=f(x) ~~ ,~~~~~~0<x\le 1, \end{aligned}$$
(1.1)

where

$$\begin{aligned} {\mathcal {L}}(x,t) =\frac{4xt^{2}}{\pi (x^{2}+t^{2})^{2}}. \end{aligned}$$
(1.2)

This is an integral equation of some special type since the kernel \({\mathcal {L}}(x,t)\) has singularity at (0, 0) only. f(x) is a prescribed function relating to the internal pressure given by

$$\begin{aligned} f(x) =\frac{\Gamma (\frac{\sigma }{2})}{\sqrt{\pi }\Gamma (\frac{\sigma +1}{2})}x^{\sigma -1} . \end{aligned}$$
(1.3)

Since the cracks are in the shape of a cross, the problem is known as the cruciform crack problem. Of interest here, is the stress intensity factor \(\chi (1)\) which is directly proportional to stress intensity at the crack trip. How the integral equation (1.1) occurs in the problem of cruciform crack is explained by Stallybrass [8] who solved the integral equation in a closed form using Wiener–Hopf technique and provided the numerical results for the stress intensity factor. Rooke and Sneddon [6] solved this integral equation approximately by using an expansion in terms of Legendre functions and obtained numerical results which are very close to those of Stallybrass [8] although the convergence is slow. The two methods appear to be somewhat elaborate. The integral equation (1.1) has also been solved numerically by various other methods from time to time. For example, Elliot [2] employed the method of Sigmoidal transformation to obtain approximate solution for the case \(f(x)=1\). It is not obvious if this method is useful for other forms of f(x). Tang and Li [9] solved the integral equation approximately by employing Taylor series expansion for the unknown function and obtained very accurate numerical estimates for the stress intensity factor. They made use of Cramer’s rule in the mathematical analysis so that if one increases the number of terms in the approximation the calculation becomes unwieldy so as to make the method unattractive. Bhattacharya and Mandal [1] solved the integral equation approximately by two different methods, one is based on expansion of the unknown function in terms of Bernstein polynomials and the other is based on expansion in terms of rationalized Haar functions. Singh and Mandal [7] also solved it by using Legendre multi-wavelets. All these methods provide numerical results for \(\chi (1)\) which are very close to exact results given by Stallybrass [8]. Expansion in terms of Bernstein polynomials or Haar functions or Legendre multi-wavelets suggest expansion in terms of other functions such as Daubechies scale functions since these provide a somewhat new tool in the numerical solution of integral equations.

In this paper, Daubechies scale functions are employed to expand the unknown function \(\chi (x)\). K-Daubechies scale function is employed to find approximate solution of integral equation taking \(K=3\). It may be noted that \(K=1\) corresponds to Haar wavelets. As the result can be improved taking larger value of K, so the results obtained by using K-Daubechies scale function are better than the results using the rationalized Haar functions. Though Legendre multi-wavelets give satisfactory results, K-Daubechies scale function has some interesting features like compact support, fractal nature and no explicit form at all resolutions. Only the knowledge of the low-pass filter coefficients in two-scale relation is required throughout the calculation. For these reasons, Daubechies scale function is used as an efficient and new mathematical tool to solve integral equations. At \(x=1\), the expansion of \(\chi (x)\) reduces to a finite expansion because most of Daubechies scale functions vanish. Actually, the integral equation (1.1) produces a system of linear equations in the unknown coefficients. After solving this linear system, the unknown function \(\chi (x)\) is evaluated at \(x=1\) so as to obtain numerically the value of the stress intensity factor. For different values of \(\sigma\) in the expression of internal pressure f(x) given by (1.3), \(\chi (1)\) is obtained and compared to known results available in the literature. It is found that the method is quite accurate as the approximate values of \(\chi (1)\) obtained by the present method are seen to differ negligibly from exact values.

Basic properties of Daubechies scale function and wavelets

Daubechies discovered a whole new class of compactly supported orthogonal wavelets, which is generated from a single function \(\phi (x)\), known as Daubechies scale or refinable function. K-Daubechies scale function \((K\ge 1)\) has 2K scaling coefficients and has compact support \([0, 2K-1]\). It may be noted that \(K=1\) corresponds to the Haar wavelets. The two-scale relations of scale functions and wavelet functions are given by

$$\begin{aligned} \phi \left( x\right) =\sqrt{2} \sum _{l=0}^{2K-1} h_{l}\phi \left( 2x-l\right) \end{aligned}$$
(2.1)

and

$$\begin{aligned} \psi \left( x\right) =\sqrt{2}\sum _{l=0}^{2K-1} g_{l}\phi \left( 2x-l\right) \end{aligned}$$
(2.2)

where the low-pass filter coefficients \(h_{l}\) and the high-pass filter coefficients \(g_{l}\) are related by

$$\begin{aligned} g_{l}=\left( -1\right) ^{l}h_{2K-1-l}. \end{aligned}$$
(2.3)

A set of orthonormal basis in \({\mathbb {R}}\) is generated from the single scaling function \(\phi (x)\) using the repeated application of two operators, translation operator T and the scale transformation D, defined as

$$\begin{aligned} \phi _{n}\left( x\right) =T^{n}\phi \left( x\right) =\phi \left( x-n\right) \end{aligned}$$
(2.4)

and

$$\begin{aligned} \phi _{sn}\left( x\right) =T^{n}D^{s}\phi \left( x\right) =2^{\frac{s}{2}}\phi \left( 2^{s}x-n\right) . \end{aligned}$$
(2.5)

The scaling coefficients \(h_{l}~~(l=0, 1, 2,\ldots ,2K-1)\) for K-Daubechies scale function are determined using conditions (2.6), (2.7) and (2.8) given below,

Vanishing moment:

$$\begin{aligned} \int _{{\mathbb {R}}} x^{n}\psi \left( x\right) \mathrm{d}x=0~~~~~~n=0, 1, 2, \ldots , 2K-1. \end{aligned}$$
(2.6)

Orthonormal condition:

$$\begin{aligned} \left( \phi _{n}, \phi _{m}\right) =\int _{\mathbb {R}} \phi _{m}\left( x\right) \phi _{n}\left( x\right) \mathrm{d}x =\delta _{mn}. \end{aligned}$$
(2.7)

Normalization condition:

$$\begin{aligned} \int _{\mathbb {R}} \phi \left( x\right) \mathrm{d}x=1. \end{aligned}$$
(2.8)

It is evident that all the properties of scaling function described above are applicable on \({\mathbb {R}}\), some of which viz. translation property (2.4) is not applicable on \([a, b]\subset {\mathbb {R}}\), where a and b are integers. In order to apply the machinery of K-Daubechies scale function on a finite interval [a, b] , we have to modify the properties of scale function \(\phi _{sn}^{L ~\text {or}~ R}\) for \(2^{s}a-(2K-2)\le n \le 2^{s}a-1\) and \(2^{s}b-(2K-2)\le n\le 2^{s}b-1\), whose supports overlap partially with the finite interval [a, b]. The superscripts L and R stand for left and right partial overlaps with [a, b], respectively. By \(\phi _{sn}^{I}\), we mean the interior scale function whose supports are contained in [a, b].

Method of solution

Here, the interval [a, b] is [0, 1]. To find the approximate solution of (1.1), we approximate the unknown function \(\chi \left( x\right)\) defined in [0, 1] in terms of Daubechies scale functions in the form

$$\begin{aligned} \chi \left( x\right) =\sum _{n}\alpha _{n}\phi _{sn}\left( x\right) . \end{aligned}$$
(3.1)

Here, \(\phi _{sn}\left( x\right)\) represents the translate of scaling function on a sufficiently fine scale s.

Using (3.1) in the integral equation (1.1), we obtain

$$\begin{aligned} \sum _{n}\alpha _{n}\phi _{sn}\left( x\right) + \sum _{n}\alpha _{n}\int _{0}^{1}\phi _{sn}\left( t\right) {\mathcal {L}}(x,t) \mathrm{d}t=f\left( x\right) . \end{aligned}$$
(3.2)

Using the orthonormal condition of \(\phi _{sm}\) for different values of m, Eq. (3.2) is reduced to the form

$$\begin{aligned} \sum _{n}\alpha _{n}\left[ N_{mn}+J_{mn}\right] =B_{m}, \end{aligned}$$
(3.3)

with

$$\begin{aligned} B_{m}= & {} \int _{0}^{1}f\left( x\right) \phi _{sm}\left( x\right) \mathrm{d}x, \end{aligned}$$
(3.4)
$$\begin{aligned} N_{mn}= & {} \int _{0}^{1} \phi _{sn}\left( x\right) \phi _{sm}\left( x\right) \mathrm{d}x, \end{aligned}$$
(3.5)
$$\begin{aligned} J_{mn}= \int _{0}^{1}\int _{0}^{1} {\mathcal {L}}(x,t) \phi _{sn}\left( t\right) \phi _{sm}\left( x\right) \mathrm{d}x\mathrm{d}t. \end{aligned}$$
(3.6)

Here, the range of m and n depends on scales s and K . As the support of \(\phi \left( x\right)\) is \(\left[ 0,2K-1\right]\) and \(K(\ge 1)\) is a positive integer, so \(m, n \in \left[ -(2K-2),2^{s}-1\right]\).

After solving the system (3.3), the unknown constants \(\alpha _{n}\) are obtained. Now, from (3.1), we get

$$\begin{aligned} \chi \left( 1\right) =2^{\frac{s}{2}}\sum _{n}\alpha _{n}\phi \left( 2^{s}-n\right) . \end{aligned}$$
(3.7)

As the support of \(\phi \left( x\right)\) is \(\left[ 0,2K-1\right]\), the values of n in (3.7) satisfy the range \(1-2K+2^{s}\le n\le 2^{s}\). If we take \(( K=3)\)-Daubechies scale function, we need values of \(\phi \left( p\right)\) for \(p = 5, 4, 3, 2, 1, 0\). The detailed tricks for calculation of (3.4), (3.5) and (3.6) are described below.

Using the explicit form of f(x) in (1.3) and the two-scale relation (2.1), the expression in (3.4) reduces to the form

$$\begin{aligned} B_{m}=2^\frac{\left( 1-2\sigma \right) s }{2}\frac{\Gamma \left( \frac{\sigma }{2}\right) }{\sqrt{\pi }\Gamma \left( \frac{\sigma +1}{2}\right) }\int _{-m}^{2^{s}-m}\left( x+m\right) ^{\sigma -1} \phi \left( x\right) \mathrm{d}x. \end{aligned}$$
(3.8)

Now, using the Gauss-type quadrature rule with complex nodes and weights for integrals involving Daubechies scale function (cf. Panja and Mandal [5]), we obtain

$$\begin{aligned} B_{m}=2^\frac{\left( 1-2\sigma \right) s }{2}\frac{\Gamma \left( \frac{\sigma }{2}\right) }{\sqrt{\pi }\Gamma \left( \frac{\sigma +1}{2}\right) }\sum _{i=1}^{N}w_{i}b_{m}\left( x_{i}\right) \end{aligned}$$
(3.9)

where

$$\begin{aligned} b_{m}\left( x\right) =\left( x+m\right) ^{\sigma -1}. \end{aligned}$$
(3.10)

The determination of the nodes \(x_{i}\) and weights \(w_{i}\) is described by Panja and Mandal [5].

The basic trick for the calculation of the integral (3.5) is described by Kessler et al. [3] and Panja and Mandal [4]. If \(\phi (x)\) is the scale function with compact support \([0, 2K-1] ~ (K\ge 1)\), then it produces a system of orthonormal basis \(\phi _{sn}\) given by (2.5). From (3.5), we get

$$\begin{aligned} N_{mn}=\int _{0}^{2^{s}}\phi _{m}\left( x\right) \phi _{n}\left( x\right) \mathrm{d}x. \end{aligned}$$
(3.11)

If \(m, n = -(2K-2),-(2K-3),\ldots ,-1\), we denote \(N_{mn}\) by \(N_{mn}^{L}\) and \(\phi _{m}\) by \(\phi _{m}^{L}\). Again if \(m,n =2^{s}-(2K-2),2^{s}-(2K-3),\ldots ,2^{s}-1\), we denote \(N_{mn}\) by \(N_{mn}^{R}\) and \(\phi _{m}\) by \(\phi _{m}^{R}\). The two-scale relation (2.1) produces

$$\begin{aligned} \phi _{n}\left( x\right) =\sqrt{2} \sum _{l_{1}=0}^{2K-1} h_{l_{1}}\phi \left( 2x-2n-l_{1}\right) . \end{aligned}$$
(3.12)

From (3.11), one gets the recursion relation

$$\begin{aligned} N_{mn}^{L ~\text {or}~ R}=\sum _{l_{1}=0}^{2K-1}\sum _{l_{2}=0}^{2K-1}h_{l_{1}} h_{l_{2}} N_{2m+l_{1}2n+l_{2}}^{L~\text {or}~ R}. \end{aligned}$$
(3.13)

Here, \(h_{l}~(l=0, 1, 2,\ldots ,2K-1)\) are the low-pass filters. The recursion relation (3.13) together with \(N_{mn}^{L}=0\) when \(m~\text {or}~ n \le -(2K-1)\) or \(\mid m-n \mid \ge 2K-1\) and \(N_{mn}^{R}=0\) when m or \(n \ge 2^{s}\) or \(\mid m-n \mid \ge 2K-1\) gives all the values of \(N_{mn}^{L~\text {or}~R}\) for \(m,n= -(2K-2), -(2K-3), \ldots ,-1\) and \(2^{s}-(2K-2), 2^{s}-(2K-3), \ldots ,2^{s}-1\). Here, the superscripts L and R stand for left edge and right edge of [0, 1], respectively. Again if \(0\le m,n \le 2^{s}-(2K-1)\), we mean \(N_{mn}\) by \(N_{mn}= N_{mn}^{I}=\delta _{mn}\). For \(( K=3)\)-Daubechies scale function, the numerical values of \(N_{mn}^{L ~\text {or}~ R}\) and the corresponding values of the inverse of the matrix formed by the elements \(N_{mn}^{L~\text {or}~ R}\) are tabulated in different tables by Panja and Mandal [4].

Now, the calculation of (3.6) is described. Using the two-scale relation (2.1), from (3.6), we obtain

$$\begin{aligned} J_{mn}=\frac{4}{\pi }2^{s}\int _{0}^{1}\int _{0}^{1} \frac{xt^{2}}{(x^{2}+t^{2})^{2}}\phi \left( 2^{s}t-n\right) \phi \left( 2^{s}x-m\right) \mathrm{d}x\mathrm{d}t =\dfrac{4}{\pi } I_{m,n}^{s} \end{aligned}$$
(3.14)

where \(I_{m,n}^{s}\) is given by the relation

$$\begin{aligned} I_{m,n}^{s}=\int _{0}^{2^{s}}\int _{0}^{2^{s}} \frac{xt^{2}}{(x^{2}+t^{2})^{2}}\phi _{m}\left( x\right) \phi _{n}\left( t\right) \mathrm{d}x\mathrm{d}t. \end{aligned}$$
(3.15)

With the help of two-scale relation (2.1) and (3.12), the expression for \(I_{m,n}^{s}\) becomes

$$\begin{aligned} I_{m,n}^{s}= \sum _{l_{1}=0}^{2K-1} \sum _{l_{2}=0}^{2K-1}h_{l_{1}}h_{l_{2}}\int _{-2m-l_{1}}^{2^{s+1}-2m-l_{1}} \int _{-2n-l_{2}}^{2^{s+1}-2n-l_{2}}\Theta _{m,n,l_{1},l_{2}}\left( x,t\right) \phi \left( x\right) \phi \left( t\right) \mathrm{d}x\mathrm{d}t \end{aligned}$$
(3.16)

where

$$\begin{aligned} \Theta _{m,n,l_{1},l_{2}}\left( x,t\right) =\frac{\left( x+2m+l_{1}\right) \left( t+2n+l_{2}\right) ^{2}}{\left[ \left( x+2m+l_{1}\right) ^{2}+\left( t+2n+l_{2}\right) ^{2}\right] ^{2}}. \end{aligned}$$
(3.17)

Using the Gauss quadrature rule involving the Daubechies scale function, (3.16) is reduced to the form

$$\begin{aligned} I_{m,n}^{s}= \sum _{l_{1}=0}^{2K-1} \sum _{l_{2}=0}^{2K-1} \sum _{i_{1}=1}^{N}\sum _{i_{2}=1}^{N}h_{l_{1}}h_{l_{2}}w_{i_{1}}w_{i_{2}} \Theta _{m,n,l_{1},l_{2}}\left( x_{i_{1}},t_{i_{2}}\right) . ~~~~ \end{aligned}$$
(3.18)

If \(-\left( 2K-2\right) \le m,n\le 0\), then \(\Theta _{m,n,l_{1},l_{2}}\left( x,t\right)\) in (3.17 ) has singularity at \(\left( 0,0\right)\), and for these values of m, n, the values of \(I_{m,n}^{s}\) cannot be determined using the relation (3.18).

If \(-\left( 2K-2\right) \le m,n\le 0\), the recursion relation for \(I_{m,n}^{s}\) is obtained as

$$\begin{aligned} I_{m,n}^{s}= & {} \sum _{l_{1}=0}^{2K-1} \sum _{l_{2}=0}^{2K-1}h_{l_{1}}h_{l_{2}} I_{2m+l_{1},2n+l_{2}}^{s}\nonumber \\&+ \sum _{l_{1}=0}^{2K-1} \sum _{l_{2}=0}^{2K-1}h_{l_{1}}h_{l_{2}}\int _{2^{s}-2m-l_{1}}^{2^{s+1}-2m-l_{1}}\int _{2^{s}-2n-l_{2}}^{2^{s+1}-2n-l_{2}}\Theta _{m,n,l_{1},l_{2}}\left( x,t\right) \phi \left( x\right) \phi \left( t\right) \mathrm{d}x\mathrm{d}t\nonumber \\&+\sum _{l_{1}=0}^{2K-1} \sum _{l_{2}=0}^{2K-1}h_{l_{1}}h_{l_{2}}\int _{-2m-l_{1}}^{2^{s}-2m-l_{1}}\int _{2^{s}-2n-l_{2}}^{2^{s+1}-2n-l_{2}}\Theta _{m,n,l_{1},l_{2}}\left( x,t\right) \phi \left( x\right) \phi \left( t\right) \mathrm{d}x\mathrm{d}t\nonumber \\&+\sum _{l_{1}=0}^{2K-1} \sum _{l_{2}=0}^{2K-1}h_{l_{1}}h_{l_{2}}\int _{2^{s}-2m-l_{1}}^{2^{s+1}-2m-l_{1}}\int _{-2n-l_{2}}^{2^{s}-2n-l_{2}}\Theta _{m,n,l_{1},l_{2}}\left( x,t\right) \phi \left( x\right) \phi \left( t\right) \mathrm{d}x\mathrm{d}t. \end{aligned}$$
(3.19)

Again using the Gauss quadrature rule involving the Daubechies scale function, (3.19) is reduced to the form

$$\begin{aligned} I_{m,n}^{s}= & {} \sum _{l_{1}=0}^{2K-1} \sum _{l_{2}=0}^{2K-1} h_{l_{1}} h_{l_{2}} I_{2m+l_{1},2n+l_{2}}^{s}\nonumber \\&+ \sum _{l_{1}=0}^{2K-1} \sum _{l_{2}=0}^{2K-1}\sum _{i_{1}=1}^{N}\sum _{i_{2}=1}^{N}h_{l_{1}} h_{l_{2}}w^{'}_{i_{1}}w^{'}_{i_{2}}\Theta _{m,n,l_{1},l_{2}}\left( x^{'}_{i_{1}},t^{'}_{i_{2}}\right) \nonumber \\&+\sum _{l_{1}=0}^{2K-1}\sum _{l_{2}=0}^{2K-1}\sum _{i_{1}=1}^{N} \sum _{i_{2}=1}^{N}h_{l_{1}}h_{l_{2}}w^{''}_{i_{1}}w^{'}_{i_{2}} \Theta _{m,n,l_{1},l_{2}}\left( x^{''}_{i_{1}},t^{'}_{i_{2}}\right) \nonumber \\&+\sum _{l_{1}=0}^{2K-1}\sum _{l_{2}=0}^{2K-1} \sum _{i_{1}=1}^{N}\sum _{i_{2}=1}^{N}h_{l_{1}}h_{l_{2}}w^{'}_{i_{1}}w^{''}_{i_{2}} \Theta _{m,n,l_{1},l_{2}}\left( x^{'}_{i_{1}},t^{''}_{i_{2}}\right) . \end{aligned}$$
(3.20)

Here, \(h_{l_{j}}\) for \(j=1, 2\)\((l_{j}=0, 1, 2\ldots ,2K-1)\) are the low-pass filters. Basic trick for calculating weights \(w^{'}_{i_{j}}\), \(w^{''}_{i_{j}}\) for \(j=1, 2\) and nodes \(x^{'}_{i_{1}}\), \(x^{''}_{i_{1}}\) and \(t^{'}_{i_{2}}\), \(t^{''}_{i_{2}}\) with a program in MATHEMATICA has been discussed by Panja and Mandal [5]. Also, \(I_{m,n}^{s}=0\) for m or \(n \le -(2K-1)\) or m or \(n \ge 2^{s}\). We present here the numerical values of \(I_{m,n}^{s}\) for( \(K=3\))-Daubechies scale functions taking \(s =3\) for those values of m and n for which \(\Theta _{m,n,l_{1},l_{2}}\left( x,t\right)\) has singularity at (0, 0). Table 1 shows the values of \(I_{m,n}^{s}\) for \(N=5,\) whereas Table 2 shows the values of \(I_{m,n}^{s}\) for \(N=7\).

Table 1 Numerical values of \(I_{m,n}^{s}\) for \(K=3, N=5\)
Full size table
Table 2 Numerical values of \(I_{m,n}^{s}\) for \(K=3, N=7\)
Full size table

Numerical results

A comparison between the numerical values of \(\chi (1)\) obtained here by using Daubechies scale functions and exact results of Stallybrass [8] is given in Table 3 for different values of \(\sigma =1, 2, 3, \ldots ,10\).

Table 3 Approximate values of \(\chi \left( 1\right)\) for different \(\sigma\)
Full size table

The table shows the exact values of \(\chi (1)\)\((\sigma =1, 2, \ldots ,10)\) according to Stallybrass [8] and the results obtained by the present method with their relative errors.

The numerical results are displayed in Fig. 1(a–e). For the sake of clarity, five figures are drawn, wherein the stress intensity factor \(\chi (1)\) is depicted against the parameter \(\sigma\) for different integral values. In each figure, \(\chi (1)\) obtained from Stallybrass’s [8] exact result is denoted by the symbol “\(\square\),” and \(\chi (1)\) obtained from other approximate methods is shown. The figures are self-explanatory. However, as the result obtained by Sigmoidal transformation method is only for \(\sigma =1\), this is not shown here. From these figures, it is obvious that all the methods including the present method provide very accurate results.

Fig. 1
figure 1

\(\chi (1)\) against \(\sigma\)

Full size image

Conclusion

Here, a numerical scheme based on expansion in terms of K-Daubechies scale function is employed for obtaining approximate numerical estimates of an integral equation of some special type arising in the classical problem of cruciform crack in elasticity. Comparison between the numerical results obtained by the present method with the exact results obtained by Stallybrass [8] shows that the method is quite accurate. The method works nicely for moderate values of K (e.g., \(K=3\)). The results can be further improved taking larger values of \(K ~(K>3\)).