5
Numerical method for solving an elliptic
Equation
OLGA MARTIN
Department of Mathematics
University “Politehnica” of Bucharest
Splaiul Independentei 313,Bucharest 16
ROMANIA
Abstract: - In this paper we provide a periodical solution of a Dirichlet problem using the Fast Fourier Transform (FFT). The presented algorithm allows a considerable reduction of the operations number necessary to solve a partial differential equation, which appear in the potential theory. A detailed analysis of the approximation of periodical solutions obtained for the numerical examples is achieved.
Key-Words: periodical solution, finite difference methods, Fast Fourier Transform (FFT), interpolation polynomial, interpolation error.
5
1. Introduction
We consider the Dirichlet problem on a square domain D for the equation:
(1)
, (2)
where ¶D is the boundary of D. To obtain a periodical solution with respect to variable x, we use an algorithm which relies on the finite difference method and Fast Fourier Transform (FFT).
We study the dependence of solution on the form of f (x,y), on the step of the network and compare the approximate and true solutions of the numerical examples.
This numerical method is an extension of the theoretical results which was obtained by Cooley, Tuckey [1], Marciuk [6] and Rombaldi [7], for stationary problems. Unlike, the other papers which used the FFT, the numerical examples prove that the algorithm is more flexible and computationally more efficient for the solution of difference equations and to define an interpolation polynomial.
2. Problem formulation
The network on D has a step h = 1/N, nodes xk = k/ N, yl = l/ N and
DN = {(xk, yl) | 0£ k £ N, 0£ l£ N }. Let ¶D be the set of nodes of the form (x0, yl), (xl, y0), (xN, yl), (xl, yN) and let function j(x) continuous on D be the solution of (1) - (2) which is approximated by a reticular function j h. Defining , the Laplacian operator in node (xk, yl) is
and (3)
Once the step h was established we denote the reticular function jh with j and problem (3) becomes
(4)
. (5)
Introducing the notations
(6)
and defining the matrix of (N - 1) order
we obtain the following set of equations
(8)
with B = A + (2 + mh2)IN - 1, IN - 1 being the identity matrix of (N - 1) order.
Since the eigenvalues of the matrix A (respectively B) are distinct, the corresponding eigenvectors are linearly independent. Using the relationship between matrices A and B a common basis is identified consisting of the eigenvectors
{v(1), v(2),..., v(N -1)}, where , being the k component of the vector v (m). The factor has been chosen since the basis is orthonormal in the space RN-1. In this space the vectors j l and fl , l = 1,2,…, N-1 are represented as:
(9)
(10)
Using these expressions in (8) and multiplying with the vector v (m), one obtains for each fixed m a system of equations with a band coefficients matrix
(11)
where
.
To obtain the solution of (8), the Fourier coefficients of the vectors fl, , Fm,l, with
l = 1, 2, ..., N-1 were computed by FFT and (11) was solved (N - 1) times (for each ). Since the vectors v (m) are orthonormal, the Fourier coefficients Fm,l have the following form:
(12)
For the periodic function f having the period Tx=1, the relations (12) coincide with those which define the bj coefficients of Fourier interpolation polynomial for f. Its values in the nodes xj are
(13)
where
(14) .
Using (12) and (14) it follows that
, m = 1,2,…., N (15)
for each fixed l.
Now, to obtain the coefficients bm, a network on [0,1] with 2N nodes was considered for each fixed l and f0,l = fN,l =.....= f2N-1,l = 0. To halve the computations number one separates the components of fj,l with even index from those with odd index for each l. Let
uj = f2j + i.f2j+1, j = 0, 1,…, N - 1 (16)
and define Direct Discret Fourier Transform (DDFT)
(17)
with , j Î{0,1,.., N-1}. DDFT can be expressed matriceally as:
U = T × C (18)
where
and t – the operation of transposition.
where .
Since one obtains
, where is the conjugate
matrix of T. It follows from the form of T that it is a
nonsingular matrix so that . Thus, using (17) one obtains
(19)
which define the Inverse Discrete Fourier Transform and they can be written in a matrix form as
(20)
Considering that and (20) we get
(21)
where P is a permutation matrix.
With the next theorem,[1], we find the coefficients from (21) and (15) for every l. Theorem . Let D2N be a network on the interval [0, 1], {fk}, kÎ{0, 1, ..., 2N-1} be the reticular function, where N = 2p, pÎN* and uk defined by (17). Then the coefficients of trigonometrically polynomial (13) are computed from the formulas
(22) In general, the complex products number from FFT algorithm has value N log2N.
In case when u is not defined by (16), the number of the computations decreases if the matrix T of NxN dimensions is replaced by a product of sparse matrices,[8],[9].
After we determine Fm,l , m, l =1, 2,..., N-1 we will solve (11) for each lm and from (9)
(23)
Now, we define with values (23) an interpolation polynomial p(x), xÎ[0,1], for each fixed l and
j = 1,2,..., N-1. This polynomial will be of the form (13), where
,
(24)
To obtain Ak and Bk we use FFT. For each l = 1,2,..,,N-1, let
(25)
where N* = N/2. Using the same algorithm which was presented by relations (17) – (21), where the matrix U is now , we obtain
(26)
To obtain the coefficients Ak, Bk of Pl(x) we may use (22).
4. Numerical examples
I. Let’s define a quadratic network with
N = 4, h = 1/4, m = 2 and a periodic function with period
Tx =1. In this case, the system of equations (8) can be written as
(33)
. For we obtain
.
The eigenvalues of the matrix A are
(34)
where: lk(B) = lk(A) + 2 + m.h2 = lk(A) +2,125,
k = 1, 2, 3.
Let us now find the eigenvector components corresponding to eigenvalues lk with k = 1, 2, 3. These are
(35)
and (11) will be
(36)
To obtain the solutions of (36), the Fourier coefficients of the vectors fl, lÎ{1, 2, 3} have to be computed. For each fixed l we use 2N = 8 nodes on [0,1], xj = j/2N, j = 0,1,..., 2N-1 and in this aim we define the matrix
,
.
In our case, the distribution of fk,l, k, l = 1, 2, 3 is
Tabel 1
xk / 0 / 1/4 / 1/2 / 3/4l = 1
(y = 1/4) / 0 / 81. / 0 / -81.
l = 2
(y = 1/2) / 0 / 0 / 0 / 0
l = 3
(y = 3/4) / 0 / -81. / 0 / 81.
Since the permutation matrix is obtained as
(37)
and using (21) one obtains for l = 1
.
From (22) and (15) we find for l = 1, that
b0 = b1 = b3 = 0 and b2 = 40.5, a i =0, i = 0,1,2,3 as well as
(38) Then, for l = 2 and
Now, we will identify the solutions of the system (36) for mÎ{1, 2,3}, where: l1 = 2.71;
l2 = 4.125; l3 = 5.54 and these will be the Fourier coefficients of the vectors jl.
For l1 = 2.71 the system (36) has the following form
with the solutions:.
For l2 =4.125:
,
and for l3 = 5.54
, .
From (9) and (35), the approximate solution,j and the exact solution, je of the Dirichlet problem are
Table 2
xk / 1/4 / 1/2 / 3/4j / je / j / je / j / je
l =1 / 1.19 / 1 / 0 / 0 / -1.19 / -1
l =2 / 0 / 0 / 0 / 0 / 0 / 0
l =3 / -1.19 / -1 / 0 / 0 / 1.19 / 1
II. Now will be considered m = 40 in (1) and
f (x,y) = sin2px sin2py , (x, y)ÎD (39)
In this case the Dirichlet problem is
(40)
On a given D use the same network with h =1/8
(N =8). The solution of the problem (40) is presented in Table 3
5
Table 3
k / 1 / 2 / 3 / 4 / 5 / 6 / 7l=1 / 0.0043 / 0.006 / 0.0043 / 0 / -0.004 / -0.006 / -0.004
l=2 / 0.006 / 0.0084 / 0.006 / 0 / -0.006 / -0.008 / -0.006
l=3 / 0.004 / 0.006 / 0.004 / 0 / -0.004 / -0.006 / -0.004
jl (xk) / l=4 / 0 / 0 / 0 / 0 / 0 / 0 / 0
l=5 / -0.004 / -0.006 / -0.004 / 0 / 0.0043 / 0.006 / 0.0043
l=6 / -0.006 / -0.008 / -0.006 / 0 / 0.006 / 0.0084 / 0.006
l=7 / -0.004 / -0.006 / -0.004 / 0 / 0.004 / 0.006 / 0.004
5
5
5
Now we find the interpolation polynomials for each function jl, l = 1, 2, 3 represented by the table 3. For example, when l = 1, let
and (26) becomes for N* = 8/2 = 4
From (22) we obtain the coefficients Ak, Bk
The interpolation polynomial for l = 1 is of the form
III. If in the first case on replaced f with
f(x,y) = ((4p2+2)y(y-1) – 2)sin2px (41)
we have distribution of f as follows:
Table 4
xk / 0 / 1/4 / 1/2 / 3/4l = 1 / 0 / -9.8 / 0 / 9.8
l = 2 / 0 / -12.4 / 0 / 12.4
l = 3 / 0 / -9.8 / 0 / 9.8
5
The solution of the problem (41) is presented in Table 4
Table 5
xk / 1/4 / 1/2 / 3/4j / je / j / je / j / je
l =1 / -0.21 / -0.19 / 0 / 0 / 0.21 / 0.19
l =2 / -0.29 / -0.25 / 0 / 0 / 0.29 / 0.25
l =3 / -0.21 / -0.19 / 0 / 0 / 0.21 / 0.19
6. CONCLUSIONS
If we denote g(x) = j (x,yl) for fixed l and xÎ[0, 1], then it may be expanded in a trigonometric Fourier series at any point x of this interval:
(42)
where
.
Let us now consider PN(x), the Fourier interpolation polynomial for g and we get
where
.
The increase of nodes number N of the network leads to a good approximation of the solution, which
may be verified with discrete Parseval’s equality
When N increases, the sum from the right part tends to . If the integral is bounded, the sum is bounded. Therefore, the sum of the left converge and An,Bn will be small. Thus, the trigonometric interpolation will converge for continuous function g or with a finite number of points of discontinuity.
Then, comparing the shape of the surface , where with the shape of we observe that they coincide, but theirs amplitudes are differently
References
[1]. Brigham. E. O., The Fast Fourier Transform,
Prentice-Hall, Englewood Cliffs, 1974.
[2] Dang Q., A approximate method for solving an
elliptic problem with discontinuous
coefficients, J. Comput. Appl.Math., vol.51,
No.2, 1994, pp.193-203.
[3] Godunov S.K., Reabenki V.S., The theory of
difference schemes, Editura Tehnica, Bucharest,
1977.
[4] Marciuk G. I., Shayduov V., Refining of the
solutions of the difference schemes, Mir
Publishers, Moscow, 1983.
[5] Martin O., Une méthode de résolution de
l’équation du transfert des neutrons, Rev. Roum. Sci. Tech. Méc. Apll., 37, 1992, 623 - 649.
[6] Marciuk G. I., Numerical Methods, Edit.
Academiei, Bucuresti, 1983.
[7] Rombaldi J. E., Problèmes corrigés d’analyse
numérique, Masson, Paris, 1996.