Numerical solution of Fredholm integral equation of the first kind with collocation method and estimation of error bound

Numerical solution of Fredholm integral equation of the first kind with collocation method and estimation of error bound

Applied Mathematics and Computation 179 (2006) 352–359 www.elsevier.com/locate/amc Numerical solution of Fredholm integral equation of the first kind ...

161KB Sizes 0 Downloads 38 Views

Applied Mathematics and Computation 179 (2006) 352–359 www.elsevier.com/locate/amc

Numerical solution of Fredholm integral equation of the first kind with collocation method and estimation of error bound K. Maleknejad *, N. Aghazadeh, R. Mollapourasl School of Mathematics, Iran University of Science and Technology, Narmak, Tehran 16846, Iran

Abstract Fredholm integral equation of the first kind is one of the inverse problems that arise in many engineering fields such as image processing and electromagnetic. Most inverse problems are ill-posed, so that solving discretized system of such problems with large condition number has a lot of difficulties. In this paper with using wavelet basis and collocation method the integral equation reduce to a linear system of equation. Then for solving the system we use the CG method. Furthermore, we get an estimation of error bound for this method. Ó 2005 Elsevier Inc. All rights reserved. Keywords: Fredholm integral equation; Integral equation of the first kind; Wavelet basis; Collocation method; Ho¨lder space; CG method

1. Introduction For problems in mathematical physics, in particular for initial and boundary value problems, [4] expressed three properties: 1. Existence of the solution. 2. Uniqueness of the solution. 3. Continuous dependence of the solution on the data. Hence, one wants to make sure that small errors in the data will cause only small errors in the solution. So that every problem which satisfies in the above three conditions is called well-posed problem. If we define K as a operator that is denoted by

*

Corresponding author. E-mail addresses: [email protected] (K. Maleknejad), [email protected] (N. Aghazadeh), [email protected] (R. Mollapourasl). 0096-3003/$ - see front matter Ó 2005 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2005.11.159

K. Maleknejad et al. / Applied Mathematics and Computation 179 (2006) 352–359

K: X ! Y ;

Kðf ðtÞÞ ¼

Z

353

b

kðs; tÞf ðtÞ dt a

so that we have following definition: Definition 1. Let K : X ! Y be an operator from a normed space X into a normed space Y, the equation Kðf Þ ¼ g

ð1Þ

is called well posed if K is onto, one-to-one and the inverse operator K1 : Y ! X is continuous. Otherwise the equation is called ill-posed [5]. According to this definition we may distinguish three types of ill-posedness [1]. 1. If K is not onto then Eq. (1) is not solvable for all g 2 Y (non-existence). 2. If K is not one-to-one then Eq. (1) may have more than one solution (non-uniqueness). 3. If K1 exist but is not continuous then the solution f of Eq. (1) does not depend continuously on the data g (instability).

2. Discretization by collocation method Consider the first kind Fredholm integral equation of the form: Z b kðs; tÞf ðtÞ dt ¼ gðsÞ; 1 < a 6 s 6 b < 1.

ð2Þ

a

For numerical solving of Eq. (2) we should choose a finite dimensional family of functions which the exact solution can be estimated by them. Methods that use this strategic are called projection methods, because the exact solution of equation is projected to the space with finite dimension. One of the most famous of this methods, is collocation method. For introducing this method we write the equation in the operator form: Kf ¼ g. We choose a sequence of finite dimensional subspaces X n  L2 ðRÞ for n P 1, with Xn having dimension dn. Assume that Xn has a basis of form {/1, /2, . . . ,P /d} with d  dn for notational simplicity and fn is a function d belongs to Xn, so that we can write it as fn ðtÞ ¼ j¼1 cj /j ðtÞ. By substitution into Eq. (2) we have rn ðsÞ ¼

Z

b

kðs; tÞfn ðtÞ dt  gðsÞ ¼

a

Z

b

kðs; tÞ a

d X

cj /j ðtÞ dt  gðsÞ a 6 s 6 b;

j¼1

where rn is called the residual in the approximation of the equation when using f  fn. In the operator form we have rn ¼ Kfn  g. Now to determine unknown coefficients we impose the following requirements: rn ðsi Þ ¼ 0;

i ¼ 1; 2; . . . ; d;

where si are the collocation node points. These coefficients are determined uniquely if and only if /j(t)s are being independent P [1]. In this paper, because we use Haar wavelet family which are orthogonal wavelet basis, so that fn ðtÞ ¼ dj¼1 cj /j ðtÞ is uniquely determined.

354

K. Maleknejad et al. / Applied Mathematics and Computation 179 (2006) 352–359

3. Wavelet families We can introduce wavelet basis of L2 ðRÞ space so that wjk ðtÞ ¼ 2j=2 wð2j t  kÞ;

t 2 R; j; k 2 Z;

where j, k are scaling and shifting parameters. For every f 2 L2 ðRÞ there is a unique expansion X f ðtÞ ¼ hf ; wjk iwjk ðtÞ; j;k2Z

which converges in the L2 -norm [3]; and inner product is defined by Z hf ; wjk i ¼ f ðtÞwjk ðtÞ dt. R

Moreover, if Z wjk ðtÞwj0 k0 ðtÞ dt ¼ djj0 dkk0 ;

8j; j0 ; k; k 0 2 Z

R

and wavelet basis (wjk)jk satisfy the stability condition X 2 2 2 8f 2 L2 ðRÞ; c1 kf k2 6 jhf ; wjk ij 6 c2 kf k2 j;k2Z

for some constants c1, c2 > 0, then (wjk)jk is said to be an orthonormal wavelet basis of L2 ðRÞ [2]. Also wavelets with compact supports and scaling function /(t) satisfy in refinement equation: /ðtÞ ¼

n1 X

an /ð2t  nÞ;

ð3Þ

n¼n0

where [n0, n1] is support of scaling function and (an) is finite sequence of real numbers. If /, w be scaling function and mother wavelet function respectively, then X n wðtÞ ¼ ð1Þ a1n /ð2t  nÞ. n

Definition 2 (Haar wavelet). The oldest and simplest scaling function is Haar. This scaling function is defined by  1; 0 6 x < 1; /ðxÞ ¼ 0; otherwise. We can write refinement equation as X /ðxÞ ¼ ak /ð2x  kÞ ¼ /ð2xÞ þ /ð2x  1Þ. k2Z

Haar scaling function has compact support so that there is only finite numbers of ak that are non-zero, and its support of it is [0, 1]. Definition 3 (Ho¨lder space). Ho¨lder space of order 0 < s < 1 is defined as   jf ðx þ hÞ  f ðxÞj C s ðRÞ ¼ f 2 L1 ðRÞ; sup < 1 s jhj and if s ¼ n þ s : 0 < s < 1 then   dn s s 1 n C ðRÞ ¼ f 2 L ðRÞ \ C ðRÞ; n f 2 C ðRÞ . dx

K. Maleknejad et al. / Applied Mathematics and Computation 179 (2006) 352–359

355

Theorem 1. Consider the integral equation of the first kind Z b kðs; tÞf ðtÞ dt ¼ gðsÞ; 1 < a 6 s 6 b < 1 a

assume that k(s, t) is continuous on the square [a, b]2,and the solution of the equation belong to ðC a \ L2 Þð½a; bÞ for some a > 12. And assume that the scaling function is Ho¨lder continuous of order r > a and scaling function / is compact support. For J that is positive integer SJ is a set of all non-negative integers less than 2J. Now by using P num num collocation method and projection operator P J ðf ÞðtÞ ¼ k2S J aJk /Jk ðtÞ and collocation nods si ¼ 2iJ ; i ¼ 1; 2; . . . ; 2J we have system of linear equation AJX = bJ where Z b  AJ ¼ kðsi ; tÞ/Jk ðtÞ dt ; a

k2S J

bJ ¼ ½gðsi Þ; i ¼ 1; 2; . . . ; 2J ; X T ¼ ½anum Jk k2S J if A is invertible then       X 3J   J a 2 kA1 k 6 c2 sup f ðtÞ  anum / ðtÞ 1 þ 2 .  Jk 1 J Jk  t2½a;b  k2S J

P

ðf ÞðtÞ ¼ Proof. Let f ðtÞ ’ P num J X aJk /Jk ðtÞ. P J ðf ÞðtÞ ¼

num k2S J aJk /Jk ðtÞ

and

k2S J

Now, if we substitute the approximation of f(t) with wavelet basis in the integral equation, then the right-hand side of integral equation is exchanged by a new function that we denote it by g^ðsÞ. So that, Z b gðsÞ ¼ kðs; tÞP num ðf ÞðtÞ dt; ð4Þ J a Z b kðs; tÞP J ðf ÞðtÞ dt. ð5Þ g^ðsÞ ¼ a

If we solve Eq. (5); determine the {aJk, k 2 SJ} by ðaJk Þk2S J ¼ A1 gðsi ÞÞ and if solve Eq. (4), determine the J ð½^ 1 num fanum ; k 2 S g by ða Þ ¼ A ð ½gðs Þ Þ. Consequently we have J i J Jk Jk k2S J 1 sup jaJk  anum gðsi Þ  gðsi ÞjÞ Jk j 6 kAJ k1 maxðj^

i ¼ 1; 2; . . . ; 2J .

ð6Þ

k2S J

For finding a bound for maxsi 2½a;b j^ gðsi Þ  gðsi Þj, we need to estimate the g^ðsÞ. For this result we have Rb gðsÞ ¼ a kðs; tÞf ðtÞ dt. So that Z b Z b kðs; tÞ½f ðtÞ  P J ðf ÞðtÞ dt þ kðs; tÞP J ðf ÞðtÞ dt ¼ gðsÞ; a a Z b Z b kðs; tÞP J ðf ÞðtÞ dt ¼ gðsÞ  kðs; tÞ½f ðtÞ  P J ðf ÞðtÞ dt; a a Z b kðs; tÞðf ðtÞ  P J ðf ÞðtÞÞ dt. g^ðsÞ ¼ gðsÞ  a

Then

Z  sup j^ gðsÞ  gðsÞj ¼ sup  s2½a;b

a

b

  kðs; tÞðf ðtÞ  P J ðf ÞðtÞÞ dt 6 ðb  aÞ sup ðjkðs; tÞjjf ðtÞ  P J ðf ÞðtÞjÞ.

We let M ¼ supsi ;t2½a;b jkðsi ; tÞj and have jf(t)  PJ(f)(t)j 6 cf2 max j^ gðsi Þ  gðsi Þj 6 Mcf ðb  aÞ2J a ¼ M 2 ðb  aÞ2J a

si 2½a;b

s;t2½a;b

Ja

[6], then

356

K. Maleknejad et al. / Applied Mathematics and Computation 179 (2006) 352–359

and with substitute this bound in the inequality (6), we have J a 1 . sup jaJk  anum Jk j 6 kAJ k1 M 2 ðb  aÞ2 k2S J

Also we need to determine a bound for jP J ðf ÞðtÞ  P num ðf ÞðtÞj, hence we have J    X X   num jP J ðf ÞðtÞ  P num ðf ÞðtÞj ¼  ðaJk  anum j/Jk ðtÞj J Jk Þ/Jk ðtÞ 6 kaJk  aJk k1 sup  k2S t2½a;b k2S J J X J =2 3J =2 sup j/ð2J t  kÞj ¼ kaJk  anum M 1. 6 kaJk  anum Jk k1 2 Jk k1 2 t2½a;b k2S J

Now, with using this inequalities the result of this theorem is determined. So kf ðtÞ  P num ðf ÞðtÞk1 6 kf ðtÞ  P J ðf ÞðtÞk1 þ kP J ðf ÞðtÞ  P num ðf ÞðtÞk1 J J 3

J ða2Þ 6 cf 2J a þ kA1 . J k1 M 1 M 2 ðb  aÞ2

Then   3J J a 2 kA1 k kf ðtÞ  P num ðf ÞðtÞk 6 2 c þ M M ðb  aÞ2 f 1 2 1 ; J J C ¼ maxfcf ; M 1 M 2 ðb  aÞg. 3

ðf ÞðtÞk1 6 C2J a f1 þ 22J kA1 Hence kf ðtÞ  P num J k1 g and the proof is completed. J

h

The drawback of this theorem is the error bound of the scheme contains the quantity kA1 J k1 . Hence, in the following lemma by using extra condition we found a bound for kA1 k and condition number of AJ. 1 J 0

Lemma 1. Consider the previous theorem assume that there exists J > J such that 0

k22J AJ  I S J k ¼ M J 0 < 1; 2J 0

2 where k.k is the maximum norm of rows and I S J is a identity matrix of order jSJj. Then kA1 J k 6 1M 0 also J

2ðJ J 0 Þ

condðAJ Þ 6

2 . 1  M J0

Proof. First determine a bound for kA1 J k, let 0

kBJ 0 k ¼ k22J AJ  I S J k ¼ M J 0 < 1. Considering of geometric series theorem we have 1

kðI þ BJ 0 Þ k 6

1 ; 1  kBJ 0 k 0

0

0

kðI S J þ BJ 0 Þ1 k ¼ kðI S J þ 22J AJ  I S J Þ1 k ¼ kð22J AJ Þ1 k ¼ 22J kA1 J k 6

1 1  kBJ 0 k

so that 0

kA1 J k 6

22J . 1  M J0

Now, we need a bound for kAJk.   Z X Z b    kðs ; tÞ/ ðtÞ dt ¼ max kAJ k ¼ max i Jk   i

k2S J

a

i

b 2

jkðsi ; tÞj dt a

12 X Z k2S J

a

b 2

j/Jk ðtÞj dt

12 .

ð7Þ

K. Maleknejad et al. / Applied Mathematics and Computation 179 (2006) 352–359

Because k(s, t) is continuous on the square [a, b]2, then Z b

1=2 2 M 1 ¼ max jkðsi ; tÞj dt i

k2S J

ð8Þ

a

is finite. X Z

357

b

j/Jk ðtÞj2 dt

1=2 6

a

X

0 @ðb  aÞ1=2

!1=2 1 A sup j/Jk ðtÞj2 t2½a;b

k2S J

¼ ðb  aÞ

1=2

X

"

k2S J

#1=2 sup ðj2

J =2

2

J

/ð2 t  kÞj Þ

t2½a;b

!1=2 ¼ ðb  aÞ

1=2 3J =2

sup j/ðtÞj

2

2

¼ 23J =2 M 3

ð9Þ

t2½a;b

* numerical solution - exact solution 1 0.95 0.9 0.85 0.8 0.75 0.7 0.65 0.6 0.55 0.5

0

2

4

6

8

10

12

14

16

18

Fig. 1. Results for Example 1 with scaling parameter J = 4.

* numerical solution - exact solution 1 0.95 0.9 0.85 0.8 0.75 0.7 0.65 0.6 0.55 0.5

0

10

20

30

40

50

60

70

Fig. 2. Results for Example 1 with scaling parameter J = 6.

358

K. Maleknejad et al. / Applied Mathematics and Computation 179 (2006) 352–359

with substituting (8) and (9) in Eq. (7), we have kAJ k 6 22J M 1 M 2 M 3 ¼ 22J M. Hence 0

condðAJ Þ ¼ kAJ kkA1 J k 6 M and the proof is completed.

22ðJ J Þ 1  M J0

h

4. Numerical examples In this section, for showing the accuracy and efficiency of the described method we present some examples. In the following examples for solving introduced system, we use CG method. sin s

Example 1. In this example we solve (2) with a = 0, b = 1, k(s, t) = et sin s and gðsÞ ¼  e with the exact solution f(t) = cos t by using Haar wavelet. Results showed in Figs. 1 and 2. * numerical solution

- exact solution

1

0.9

0.8

0.7

0.6

0.5

0.4

0

2

4

6

8

10

12

14

16

18

Fig. 3. Results for Example 2 with scaling parameter J = 4.

* numerical solution - exact solution 1

0.9

0.8

0.7

0.6

0.5

0.4

0

10

20

30

40

50

60

70

Fig. 4. Results for Example 2 with scaling parameter J = 6.

ð0:54 sin sþ0:84Þsin s cos s2 2

K. Maleknejad et al. / Applied Mathematics and Computation 179 (2006) 352–359

359

Example 2. In this example we solve (2) with a = 0, b = 1, k(s, t) = et(sin(s  t + 1) + 1) and g(s) = 1 + cos(s)  cos(s + 1) with the exact solution f(t) = et by using Haar wavelet with scaling parameter J = 4, 6. Results showed in Figs. 3, 4. 5. Conclusion We know that Fredholm integral equation of the first kind is a kind of importance in many engineering fields. In this paper, we reduce the integral equation to a linear system of equation using collocation method with wavelet basis. We give a bound for cond(AJ) and error of numerical solution. We show that the numerical solution of integral equation converge to the exact solution. Figs. 1–4 show the results for some examples. References [1] K.E. Atkinson, The Numerical Solution of Integral Equations of the Second Kind, Cambridge University Press, 1997. [2] A. Cohen, I. Daubechies, J.C. Feauveau, Biorthogonal basis of compactly supported wavelets, Comm. Pure Appl. Math. 45 (1992) 485–560. [3] I. Daubechies, Ten Lectures on Wavelets, vol. 61, CBMS-NSF Regional Conference Series in Applied Mathematics, SIAM, Philadelphia, 1992. [4] J. Hadamard, Lectures on Cauchy’s Problem in Linear Partial Differential Equations, Yale University Press, New Haven, 1923. [5] R. Kress, Linear Integral Equations, Springer-Verlag, Berlin, 1989. [6] A. Karoui, Wavelets: properties and approximate solution of a second kind integral equation, Comput. Math. Applicat. 46 (2003) 263– 277.