# Comparative Analysis of Methods for Regularizing an Initial Boundary Value Problem for the Helmholtz Equation

Authors: Sergey Igorevich Kabanikhin; M. A. Shishlenin; D. B. Nurseitov; A. T. Nurseitova; S. E. Kasenov
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

We consider an ill-posed initial boundary value problem for the Helmholtz equation. This problem is reduced to the inverse continuation problem for the Helmholtz equation. We prove the well-posedness of the direct problem and obtain a stability estimate of its solution. We solve numerically the inverse problem using the Tikhonov regularization, Godunov approach, and the Landweber iteration. Comparative analysis of these methods is presented.


---

## Body

## 1. Introduction

Let us consider the initial boundary value problem (continuation problem) for the Helmholtz equation in the domainΩ=(0,l)×(0,π):
(1)uxx+uyy+k2u=0,(x,y)∈Ω,ux(0,y)=0,y∈[0,π],u(0,y)=f(y),y∈[0,π],uy(x,0)=uy(x,π)=0,x∈[0,l],
where k is a given constant. To find the function u(x,y) in Ω from f(y) is required.The continuation problem is ill-posed problem; its solution is unique, but it does not depend continuously on the Cauchy data [1–10]. Note that the problem was studied by many authors. For example, Tuan and Quan [11] considered the case 0<k<1 and proposed a regularization technique which allows one to obtain a stable solution in a two-dimensional domain. Regińska and Regiński [12] showed that if k satisfies a certain condition, then the Cauchy problem for the Helmholtz equation has a stable solution in a three-dimensional domain. Isakov and Kindermann [13] used the singular value decomposition to prove that in a simple domain the considered problem becomes more stable with increasing k. The same result was obtained numerically for the general case. The uniqueness of the solution of the investigated problem was proved, for example, by Arendt and Regińska [14], where the concept of weak normal derivative was introduced in formulating the problem. In [15, 16] singular values of the continuation problem were obtained for the two-dimensional Helmholtz equation with complex wave number for simple geometry.We consider two approaches to the numerical solution of the problem (1). The first consists of formulating problem (1) in an operator form Aq=f and minimizing the coast functional J(q)=〈Aq-f,Aq-f〉 by the Landweber iteration [7]. In the second approach, problem (1) is reduced to the system of linear algebraic equations which is solved using the Tikhonov regularization and Godunov approach. In this work we present a comparative analysis of the proposed methods for the numerical solution of the problem (1).


## 2. The Direct and Inverse Problems

Let us consider the direct (well-posed) problem of finding the functionu(x,y) from the relations
(2)uxx+uyy+k2u=0,(x,y)∈Ω,(3)ux(0,y)=0,y∈[0,π],(4)u(l,y)=q(y),y∈[0,π],(5)uy(x,0)=uy(x,π)=0,x∈[0,l]
and note that the continuation problem (1) can be reduced to the inverse problem of finding function q(y) from (2)–(5) using the additional information
(6)u(0,y)=f(y),y∈[0,π].
Let us consider some theoretical results [7, 17].Definition 1.
A functionu∈L2(Ω) is called a generalized solution of the direct problem (2)–(5), if for any ω∈H2(Ω) such that
(7)ωx(0,y)=0,y∈[0,π],ω(l,y)=0,y∈[0,π],ωy(x,0)=ωy(x,π)=0,x∈[0,l],
the following equality holds:
(8)∫0l‍∫0π‍u·(ωxx+ωyy+k2ω)dydx-∫0π‍q(y)ωx(l,y)dy=0.Theorem 2 (existence of a generalized solution of the direct problem).
Ifq∈L2(0,π) and k2l2<1, then the direct problem (2)–(5) has a unique generalized solution u∈L2(Ω) and the following estimate is true:
(9)∥u∥L2(Ω)≤∥q∥L2(0,π)l1-k2l2.Proof.
Let us introduce the auxiliary problem(10)ωxx+ωyy=u-,(x,y)∈Ω,(11)ωx(0,y)=0,y∈[0,π],(12)ω(l,y)=0,y∈[0,π],(13)ωy(x,0)=ωy(x,π)=0,x∈[0,l].
Integrating the identity
(14)ωu-=ω(ωxx+ωyy)=(ωxω)x-ωx2+(ωyω)y-ωy2
over the domain Ω and considering (10)–(13), we obtain
(15)∫0l‍∫0π‍ωu-dydx=∫0π‍ωxω|0,yl,ydy+∫0l‍ωyω|x,0x,πdx-∫0l‍∫0π‍(ωx2+ωy2)dydx,
whence
(16)∫0l‍∫0π‍(ωx2+ωy2)dydx≤∥ω∥L2(Ω)·∥u∥L2(Ω).
Taking into account (12) and the equality
(17)ω(l,y)=ω(x,y)+∫xl‍ωξ(ξ,y)dξ|ω(x,y)|2=∫xl‍ωξ2(ξ,y)dξ≤∫0l‍1·ωξ2(ξ,y)dξ≤∫0l‍1dξ·∫0l‍ωξ2(ξ,y)dξ=l·∫0l‍ωξ2(ξ,y)dξ,
we have
(18)|ω(x,y)|2≤l·∫0l‍ωξ2(ξ,y)dξ,∫0l‍∫0π‍|ω(x,y)|2dydx≤∫0l‍∫0π‍(l·∫0l‍ωξ2(ξ,y)dξ)dydx∥ω∥L2(Ω)≤l2·∥ωξ∥L2(Ω).
Combining (16) and (18) yields
(19)∥ωx∥L2(Ω)≤l·∥u-∥L2(Ω),∥ω∥L2(Ω)≤l2·∥u-∥L2(Ω).
From the identities
(20)ωxu-=ωx(ωxx+ωyy)=12(ωx2)x+ωxωyy=12(ωx2)x+(ωxωy)y-12(ωy2)x,
it follows that
(21)12(ωx2)x+(ωxωy)y=ωxu-+12(ωy2)x.
Integrating (21) on Ω, we get
(22)12∫0π‍ωx2(l,y)dy=∫0l‍∫0π‍ωxu-dydx+12∫0π‍ωy2(l,y)dy-12∫0π‍ωy2(0,y)dy≤∥ωx∥L2(Ω)·∥u-∥L2(Ω)≤l∥u-∥L2(Ω)2.
Due to (8), we obtain
(23)∫0l‍∫0π‍(u-2+k2ωu-)dydx≤∥q(y)∥L2(0,π)·∥ωx(l,y)∥L2(0,π).
Hence,
(24)∥u-∥L2(Ω)2≤∥q(y)∥·∥ωx(l,y)∥+k2∥ω∥·∥u-∥≤∥q(y)∥·l∥u-∥+k2∥u-∥·l2∥u-∥,∥u-∥L2(Ω)≤∥q(y)∥·l1-k2l2.
Thus, we have proved the well-posedness of the direct problem, which allows us to apply well-elaborated computational methods. Also, a stability estimate has been obtained inL2.


## 3. Landweber Iteration

### 3.1. Formulation in the Operator Form and Description of the Algorithm

Let us reduce the inverse problem (2)–(6) to the operator equation. Let us consider the operator A such that
(25)A:q(y)⟼u(0,y),
where u(x,y) is the solution of the direct problem (2)–(5). Then the inverse problem (2)–(6) takes the form
(26)Aq=f.
We will find the solution of the problem (26) by minimizing the following functional [7, 18–20]:
(27)J(q)=∥Aq-f∥2=∫0π‍[u(0,y;q)-f(y)]2dy
using the Landweber iteration
(28)qn+1=qn-αJ′qn,
where α∈(0,1/∥A∥2) is the descent parameter [7].Let us describe the iterative algorithm. First we choose the initial approximationq0 and we suppose that we calculated successively qn by formula (28). Assuming that we have found qn, we show below how to calculate qn+1.(1)
Solve the direct problem (2)–(5) with the known qn.(2)
CalculateJ(qn) by formula (27).(3)
Check the stopping criterionJ(qn)≤ηε2 [21]. Finish if qn meets this inequality.(4)
Solve the adjoint problem(29)ψxx+ψyy+k2ψ=0,(x,y)∈Ω,ψ(l,y)=0,y∈[0,π],ψx(0,y)=2(u(0,y)-f(y)),y∈[0,π],ψy(x,π)=ψy(x,0)=0,x∈[0,l].(5)
Calculate the gradientJ′(qn) by the formula
(30)J′q=ψx(l,y).(6)
Calculate the next approximationqn+1=qn-αJ′qn and proceed to step 1.


### 3.2. The Numerical Solution of the Direct and Adjoint Problems

The direct and adjoint problems are solved using the direct finite-difference method. For discretizing the direct problem, we construct a grid inΩ with steps hx=l/Nx, hy=π/Ny, where Nx,Ny are positive integers. Let us denote the grid by ωh={x=ihx,y=jhy;i=0,Nx¯,j=0,Ny¯}. After exchanging derivatives by finite-difference analoges with the second order, we obtain the following discrete direct problem (2)–(5):
(31)ui+1,j-2ui,j+ui-1,jhx2+ui,j+1-2ui,j+ui,j-1hy2+k2ui,j=0,i=1,Nx-1¯,j=1,Ny-1¯,u1,j-u0,jhx=0,j=1,Ny-1¯,uNx,j=qj,j=1,Ny-1¯,ui,1-ui,0hy=ui,Nx-ui,Nx-1hy=0,i=0,Nx¯.
By introducing the parameters a=1/hx2,b=1/hy2, c=k2-2a-2b, we get
(32)aui-1,j+bui,j-1+cui,j+bui,j+1+aui+1,j=0,i=1,Nx-1¯,j=1,Ny-1¯,u1,j-u0,j=0,j=1,Ny-1¯,uNx,j=qj,j=1,Ny-1¯,ui,1-ui,0=ui,Nx-ui,Nx-1=0,i=0,Nx¯.Thus, we obtain the system of algebraic equations(33)ATXT=BT,
where AT is a matrix of size (Nx+1)(Ny+1), XT is an unknown vector of the form
(34)XT=(u0,0,u0,1,u0,2,…,u0,Ny,u1,0,u1,1,u1,2,…,u1,Ny,…,uNx,0,uNx,1,uNx,2,…,uNx,Ny),
and BT is the data vector (boundary and additional conditions).Similarly, the discrete adjoint problem (29) has the form
(35)aψi-1,j+bψi,j-1+cψi,j+bψi,j+1+aψi+1,j=0,i=1,Nx-1¯,j=1,Ny-1¯,ψ1,j-ψ0,j=2(u0,j-fj),j=1,Ny-1¯,ψNx,j=0,j=1,Ny-1¯,ψi,1-ψi,0=ψi,Nx-ψi,Nx-1=0,i=0,Nx¯.
As above, this problem can be reformulated in a matrix form
(36)ATYT=B-T,
where
(37)YT=(ψ0,0,ψ0,1,ψ0,2,…,ψ0,Ny,ψ1,0,ψ1,1,ψ1,2,…,ψ1,Ny,…,ψNx,0,ψNx,1,ψNx,2,…,ψNx,Ny)
is an unknown vector and B-T is the data vector (boundary and additional conditions of the adjoint problem).


### 3.3. Results of the Numerical Experiment

Letl=1, Nx=Ny=20. We choose the parameter k=0.9. In order to test the algorithm, we assume that the exact solution has the form q(y)=1-cos⁡(2y) and calculate the corresponding additional information f. Then, let q0(y)=0.1 be the initial approximation; we try to restore the original exact solution using the Landweber iteration with α=0.01. If the data are given with an error ε, we choose the following one as a stopping criterion: J(qn)≤ηε2 [21].The computational experiment was carried out for different noise levels. Tables1, 2, and 3 show the calculation results obtained using PC Intel(R)Core(TM) i7 processor with a frequency 3.9 GHz.Table 1
Calculation results in the case of no noiseε=0.
Number
The number of iterationsn
∥
q
T
-
q
~
∥
J
(
q
)
Runtime
1
10
0.8243
0.094
4 sec
2
100
0.0632
9.42
·
1
0
-
5
40 sec
3
1000
0.0414
9.55
·
1
0
-
7
6 min 40 sec
4
5000
0.0311
1.74
·
1
0
-
8
33 min 30 sec
5
7318
0.0304
9.29
·
1
0
-
9
49 minTable 2
Calculation results in the case of noise within 1%ε=0.01.
Number
The number of iterationsn
∥
q
T
-
q
~
∥
J
(
q
)
Runtime
1
10
0.83
0.09
4 sec
2
100
0.077
4.91
·
1
0
-
4
40 sec
3
1000
0.047
3.38
·
1
0
-
4
6 min 33 sec
4
1508
0.051
3.32
·
1
0
-
4
10 minTable 3
Calculation results in the case of noise within 5%ε=0.05.
Number
The number of iterationsn
∥
q
T
-
q
~
∥
J
(
q
)
Runtime
1
10
0.86
9.07
·
1
0
-
2
4 sec
2
100
0.165
9.20
·
1
0
-
3
40 sec
3
1000
0.210
8.33
·
1
0
-
3
6 min 32 sec
4
2000
0.297
8.15
·
1
0
-
3
13 min 4 sec
5
2418
0.326
8.10
·
1
0
-
3
15 min 48 secThe approximate solution in the case ofε=0.05 is shown in Figure 1.Figure 1
▾ indicates exact solution;  ▴ indicates solution obtained by the Landweber iteration.We observe that in the case of no noise the functionalJ(q) decreases monotonically, while in the other cases the decrease stops after 100 iterations. This phenomenon can be explained by the error that arises in solving the direct and adjoint problems. Note that in the case of noise, the stopping criterion does not guarantee the minimal error in the solution of the inverse problem. However, the criterion ensures that the error is of the same order as that of the minimal one, since further calculation leads to an increase in the error.


## 3.1. Formulation in the Operator Form and Description of the Algorithm

Let us reduce the inverse problem (2)–(6) to the operator equation. Let us consider the operator A such that
(25)A:q(y)⟼u(0,y),
where u(x,y) is the solution of the direct problem (2)–(5). Then the inverse problem (2)–(6) takes the form
(26)Aq=f.
We will find the solution of the problem (26) by minimizing the following functional [7, 18–20]:
(27)J(q)=∥Aq-f∥2=∫0π‍[u(0,y;q)-f(y)]2dy
using the Landweber iteration
(28)qn+1=qn-αJ′qn,
where α∈(0,1/∥A∥2) is the descent parameter [7].Let us describe the iterative algorithm. First we choose the initial approximationq0 and we suppose that we calculated successively qn by formula (28). Assuming that we have found qn, we show below how to calculate qn+1.(1)
Solve the direct problem (2)–(5) with the known qn.(2)
CalculateJ(qn) by formula (27).(3)
Check the stopping criterionJ(qn)≤ηε2 [21]. Finish if qn meets this inequality.(4)
Solve the adjoint problem(29)ψxx+ψyy+k2ψ=0,(x,y)∈Ω,ψ(l,y)=0,y∈[0,π],ψx(0,y)=2(u(0,y)-f(y)),y∈[0,π],ψy(x,π)=ψy(x,0)=0,x∈[0,l].(5)
Calculate the gradientJ′(qn) by the formula
(30)J′q=ψx(l,y).(6)
Calculate the next approximationqn+1=qn-αJ′qn and proceed to step 1.


## 3.2. The Numerical Solution of the Direct and Adjoint Problems

The direct and adjoint problems are solved using the direct finite-difference method. For discretizing the direct problem, we construct a grid inΩ with steps hx=l/Nx, hy=π/Ny, where Nx,Ny are positive integers. Let us denote the grid by ωh={x=ihx,y=jhy;i=0,Nx¯,j=0,Ny¯}. After exchanging derivatives by finite-difference analoges with the second order, we obtain the following discrete direct problem (2)–(5):
(31)ui+1,j-2ui,j+ui-1,jhx2+ui,j+1-2ui,j+ui,j-1hy2+k2ui,j=0,i=1,Nx-1¯,j=1,Ny-1¯,u1,j-u0,jhx=0,j=1,Ny-1¯,uNx,j=qj,j=1,Ny-1¯,ui,1-ui,0hy=ui,Nx-ui,Nx-1hy=0,i=0,Nx¯.
By introducing the parameters a=1/hx2,b=1/hy2, c=k2-2a-2b, we get
(32)aui-1,j+bui,j-1+cui,j+bui,j+1+aui+1,j=0,i=1,Nx-1¯,j=1,Ny-1¯,u1,j-u0,j=0,j=1,Ny-1¯,uNx,j=qj,j=1,Ny-1¯,ui,1-ui,0=ui,Nx-ui,Nx-1=0,i=0,Nx¯.Thus, we obtain the system of algebraic equations(33)ATXT=BT,
where AT is a matrix of size (Nx+1)(Ny+1), XT is an unknown vector of the form
(34)XT=(u0,0,u0,1,u0,2,…,u0,Ny,u1,0,u1,1,u1,2,…,u1,Ny,…,uNx,0,uNx,1,uNx,2,…,uNx,Ny),
and BT is the data vector (boundary and additional conditions).Similarly, the discrete adjoint problem (29) has the form
(35)aψi-1,j+bψi,j-1+cψi,j+bψi,j+1+aψi+1,j=0,i=1,Nx-1¯,j=1,Ny-1¯,ψ1,j-ψ0,j=2(u0,j-fj),j=1,Ny-1¯,ψNx,j=0,j=1,Ny-1¯,ψi,1-ψi,0=ψi,Nx-ψi,Nx-1=0,i=0,Nx¯.
As above, this problem can be reformulated in a matrix form
(36)ATYT=B-T,
where
(37)YT=(ψ0,0,ψ0,1,ψ0,2,…,ψ0,Ny,ψ1,0,ψ1,1,ψ1,2,…,ψ1,Ny,…,ψNx,0,ψNx,1,ψNx,2,…,ψNx,Ny)
is an unknown vector and B-T is the data vector (boundary and additional conditions of the adjoint problem).


## 3.3. Results of the Numerical Experiment

Letl=1, Nx=Ny=20. We choose the parameter k=0.9. In order to test the algorithm, we assume that the exact solution has the form q(y)=1-cos⁡(2y) and calculate the corresponding additional information f. Then, let q0(y)=0.1 be the initial approximation; we try to restore the original exact solution using the Landweber iteration with α=0.01. If the data are given with an error ε, we choose the following one as a stopping criterion: J(qn)≤ηε2 [21].The computational experiment was carried out for different noise levels. Tables1, 2, and 3 show the calculation results obtained using PC Intel(R)Core(TM) i7 processor with a frequency 3.9 GHz.Table 1
Calculation results in the case of no noiseε=0.
Number
The number of iterationsn
∥
q
T
-
q
~
∥
J
(
q
)
Runtime
1
10
0.8243
0.094
4 sec
2
100
0.0632
9.42
·
1
0
-
5
40 sec
3
1000
0.0414
9.55
·
1
0
-
7
6 min 40 sec
4
5000
0.0311
1.74
·
1
0
-
8
33 min 30 sec
5
7318
0.0304
9.29
·
1
0
-
9
49 minTable 2
Calculation results in the case of noise within 1%ε=0.01.
Number
The number of iterationsn
∥
q
T
-
q
~
∥
J
(
q
)
Runtime
1
10
0.83
0.09
4 sec
2
100
0.077
4.91
·
1
0
-
4
40 sec
3
1000
0.047
3.38
·
1
0
-
4
6 min 33 sec
4
1508
0.051
3.32
·
1
0
-
4
10 minTable 3
Calculation results in the case of noise within 5%ε=0.05.
Number
The number of iterationsn
∥
q
T
-
q
~
∥
J
(
q
)
Runtime
1
10
0.86
9.07
·
1
0
-
2
4 sec
2
100
0.165
9.20
·
1
0
-
3
40 sec
3
1000
0.210
8.33
·
1
0
-
3
6 min 32 sec
4
2000
0.297
8.15
·
1
0
-
3
13 min 4 sec
5
2418
0.326
8.10
·
1
0
-
3
15 min 48 secThe approximate solution in the case ofε=0.05 is shown in Figure 1.Figure 1
▾ indicates exact solution;  ▴ indicates solution obtained by the Landweber iteration.We observe that in the case of no noise the functionalJ(q) decreases monotonically, while in the other cases the decrease stops after 100 iterations. This phenomenon can be explained by the error that arises in solving the direct and adjoint problems. Note that in the case of noise, the stopping criterion does not guarantee the minimal error in the solution of the inverse problem. However, the criterion ensures that the error is of the same order as that of the minimal one, since further calculation leads to an increase in the error.


## 4. Regularization Methods

In this section we consider a discrete analog of problem (1) and study the stability of its solution. Tikhonov regularization method and Godunov approach method are applied.


### 4.1. The Discretization of Problem (1)

We reduce the continuation problem (1) to the system of linear algebraic equations as follows [22]:
(38)AX=B,
where A is a matrix of size (Nx+1)(Ny+1), B is the data vector, and X is the desired vector of the form
(39)X=(u0,0,u0,1,u0,2,…,u0,Ny,u1,0,u1,1,u1,2,…,u1,Ny,…,uNx,0,uNx,1,uNx,2,…,uNx,Ny).Assumingl=1, Nx=Ny=50, k=0.9, and q(y)=1-cos⁡(2y), we calculate the norm and the condition number μ(A) of the matrices A and AT corresponding to the original problem (1) and the direct problem (2)–(5), respectively.The matrixA is ill-conditioned [23] (see decreasing of its singular values in Figure 4). The condition number and matrix norma of the discrete direct problem are presented in Table 4 and Figure 5. We see that the direct problem is well-posed.Table 4
The characteristics of the matricesA and AT.
Matrix
μ
(
A
)
∥
A
∥
A
6.02
·
1
0
17
321098.0
A
T
3.42
·
1
0
6
321098.0In view of the ill-conditioning of the matrixA, that is, the ill-posedness of the original problem, we will use regularization methods.


### 4.2. Tikhonov Regularization

The Tikhonov regularization consists of replacing the systemAX=B by the system (αI+A*A)X~=A*B [7]. We choose the regularization parameter α, minimizing the discrepancy ∥Xαi∥-∥Xαi-1∥ according to [24].As above, we putNx=Ny=20,k=0.9, and q(y)=1-cos⁡(2y) and calculate with ε=0, ε=0.01, and ε=0.05. The approximate solution in the case of ε=0.05 is shown in Figure 2.Figure 2
▾ indicates exact solution q(y) ;  ▴ indicates solution obtained by the Tikhonov regularization.


### 4.3. Godunov Approach

S. K. Godunov proposed considering the extended system(40)[(1-α)AαD]X=[(1-α)f0],
whereas D contains some a priori information concerning the inverse problem solution. We take as a priori information the existence of the second derivative of the solution [25].We chooseα to minimize ∥Xαi∥-∥Xαi-1∥.ForNx=Ny=20, k=0.9, and q(y)=1-cos⁡(2y), we calculate with different values of ε. The result is shown in Figure 3.Figure 3
▾ indicates exact solution q(y);  ▴ indicates solution obtained by the Godunov regularization.Figure 4
Singular values of matrixA for N=50.Figure 5
Singular values of matrixAT for N=50.


## 4.1. The Discretization of Problem (1)

We reduce the continuation problem (1) to the system of linear algebraic equations as follows [22]:
(38)AX=B,
where A is a matrix of size (Nx+1)(Ny+1), B is the data vector, and X is the desired vector of the form
(39)X=(u0,0,u0,1,u0,2,…,u0,Ny,u1,0,u1,1,u1,2,…,u1,Ny,…,uNx,0,uNx,1,uNx,2,…,uNx,Ny).Assumingl=1, Nx=Ny=50, k=0.9, and q(y)=1-cos⁡(2y), we calculate the norm and the condition number μ(A) of the matrices A and AT corresponding to the original problem (1) and the direct problem (2)–(5), respectively.The matrixA is ill-conditioned [23] (see decreasing of its singular values in Figure 4). The condition number and matrix norma of the discrete direct problem are presented in Table 4 and Figure 5. We see that the direct problem is well-posed.Table 4
The characteristics of the matricesA and AT.
Matrix
μ
(
A
)
∥
A
∥
A
6.02
·
1
0
17
321098.0
A
T
3.42
·
1
0
6
321098.0In view of the ill-conditioning of the matrixA, that is, the ill-posedness of the original problem, we will use regularization methods.


## 4.2. Tikhonov Regularization

The Tikhonov regularization consists of replacing the systemAX=B by the system (αI+A*A)X~=A*B [7]. We choose the regularization parameter α, minimizing the discrepancy ∥Xαi∥-∥Xαi-1∥ according to [24].As above, we putNx=Ny=20,k=0.9, and q(y)=1-cos⁡(2y) and calculate with ε=0, ε=0.01, and ε=0.05. The approximate solution in the case of ε=0.05 is shown in Figure 2.Figure 2
▾ indicates exact solution q(y) ;  ▴ indicates solution obtained by the Tikhonov regularization.


## 4.3. Godunov Approach

S. K. Godunov proposed considering the extended system(40)[(1-α)AαD]X=[(1-α)f0],
whereas D contains some a priori information concerning the inverse problem solution. We take as a priori information the existence of the second derivative of the solution [25].We chooseα to minimize ∥Xαi∥-∥Xαi-1∥.ForNx=Ny=20, k=0.9, and q(y)=1-cos⁡(2y), we calculate with different values of ε. The result is shown in Figure 3.Figure 3
▾ indicates exact solution q(y);  ▴ indicates solution obtained by the Godunov regularization.Figure 4
Singular values of matrixA for N=50.Figure 5
Singular values of matrixAT for N=50.


## 5. Comparative Analysis of Methods

Tables5, 6, and 7 present the results of the numerical solution of the problem (1) with different levels of noise. We see that the Godunov method is more accurate compared to the other methods. The same is demonstrated in Figure 6.Table 5
Comparative analysis of methods in the case of no noise.
Number
Name methods
∥
q
T
-
q
~
∥
Runtime
1
Landweber iteration
0.030
49 min
2
Tikhonov regularization
0.021
23 sec
3
Godunov regularization
0.019
46 secTable 6
Comparative analysis of methods in the case of noise within 1%.
Number
Name methods
∥
q
T
-
q
~
∥
Runtime
1
Landweber iteration
0.051
10 min
2
Tikhonov regularization
0.190
15 sec
3
Godunov regularization
0.055
17 secTable 7
Comparative analysis of methods in the case of noise within 5%.
Number
Name methods
∥
q
T
-
q
~
∥
Runtime
1
Landweber iteration
0.5349
18 min 24 sec
2
Tikhonov regularization
0.4088
7 sec
3
Godunov regularization
0.2989
12 secFigure 6
▾ indicates exact solution q(y); ● indicates Landweber iteration; ▪ indicates Tikhonov regularization; ▴ indicates Godunov regularization.



---

*Source: 786326-2014-09-08.xml*