1
VARIATIONAL SOLUTION OF THE
TWO-DIMENSIONAL TRANSPORT EQUATION
Olga MARTIN
Department of Mathematics
University “Politehnica” of Bucharest
Splaiul Independentei 313,Bucharest 16
ROMANIA
Abstract: - The paper presents a numerical method for two-dimensional transport equations, which use the variational calculus. In order to find out the solution of the stationary regime we consider a function u, which satisfies the boundary conditions and minimizes a given functional J(u). The numerical solution is in a good concordance with the analytical solution.
Key-Words: transport equation, variational calculus,Euler-Lagrange equation.
1 Problem Formulation
1
Let us now consider the integro-differential equation of the transport theory
1
(1)
with
- the boundary of D1(2)
1
where u( x,y,,) is the density of the neutrons number and , are the components of the velocity in the x and y directions, respectively.
In order to solve the problem
(1) – (2) with the help of the variational methods we should replace this with the equation obtained by addition of the partial derivatives of (1) with respect to the variable x and y, respectively, multiplied by and
1
(3)
or
(4)
where
(5)
1
in accordance with (1).
The numerical method for solving the problem (4) with boundary condition (2) is based on the discretization D of the domain D. This means it is splitted into small(elementary) regions
1
.
For example
1
and we assumed that the and have in the points of Dij the constant values: i, j.
Next, we shall show that the function, u(x,y,i,j), minimizing a functional
1
(6)
1
is a solution of the relation (4) (or (1)).
Indeed, the function u which leads to an extremum of J(u) is a solution of the Euler-Lagrange equation, which may be written for (4) in the form
1
(7)
1
and thus we obtain the equation (4) for .
Let us now consider for D a square network with step h and then each square will be splitted into two isosceles triangles (fig.1).
The solution of the equation
1
(8)
1
will be presented in the form of the serie
1
(9)
1
where the basis is formed by pyramidal functions, , which correspond to the
N – hexagonal regions, with the center in the point Pkl of the network. It should be noted that the points Pkldon’t belong to the boundary contour. Each hexagonal subdomain is the union of six triangles, Dkl,j, j1,2,…,6 and their numbering is shown in fig.1.
1
OO
Fig.1 Fig.2
1
Then the function can be written in the following form
1
(10)
1
The geometric interpretation of the function is a pyramid with the height equals to unity and P is the projection of its apex (fig.2).
1
In order to minimize the functional J(u) into D
1
we find uh, hence the coefficients ij by solving the following system
(11)
or in matrix for
(12)
where A is a banded matrix because the non-zero elements from a band along the main diagonal. The elements of A are expressed in the form
1
(13)
1
where are noted x1 = x, x2 = y and the coefficients Cstare given by the equation (8) which was written in the form
(14)
The vector with , have
(15)
but is the unknown vector with
For example, when k = 1, l = 1,
m = 1, n = 2
1
(16)
1
which correspond to D11,3, D11,4, D12,1, D12,6.
Solving the equation (12) we find for the values (i,j)D,
. It should be noted that the solution of the problem (4), (2) is an element in the space L2(D1). This space is a set of functions defined in a measurable set D1 with the properties
with the scalar product and .
1
As it had been mentioned in [4] an approximate solution uhof (the class of the functions with a continuous third derivative) verifies the following inequality
1
(17)
1
where M is an upper bound for the derivatives of the third order of u.
is a lower bound for the angles of the triangles which form the D1. If then
1
1
.(18)
2 Numerical results
1
Let us consider the transport equation
1
(19)
1
with the conditions (2). For the domain D [0,1][0,1] we use the network D mentioned before and for D1 we have: (20)
1
O
Fig. 3
1
We begin the construction of an approximate solution by considering the sum
1
(21)
1
with
(22)
Substituting the sum (21) into (4) we get for a point (i,i)
1
(23)
1
and in accordance with (19)
1
(24)
1
We shall now determine the matrices A and g of the equation (12) in the form
1
with
1
and because 0 only for we obtain
1
(25)
1
The calculations were made for the corresponding values (i,i). The approximate solution u is given in the table and compared with the exact solution .
1
Table1
1
(,) / NOD /u
/ / (,) / NOD / u / = 0
= 1/2 / (11) / 0.44 / 0.43 / = 1/4
= -1/4 / (11) / -0.335 / -0.216
(12) / -0.45 / -0.43 / (12) / 0.26 / 0.216
(21) / 0.45 / 0.43 / (21) / -0.26 / -0.216
(22) / -0.44 / -0.43 / (22) / 0.335 / 0.216
= 1/2
= 1/2 / (11) / 0.33 / 0.43 / = 1/4
= 1/4 / (11) / 0.14 / 0.216
(12) / -0.37 / -0.43 / (12) / -0.17 / -0.216
(21) / 0.37 / 0.43 / (21) / 0.17 / 0.216
(22) / -0.33 / -0.43 / (22) / -0.14 / -0.216
= 1/2
= -1/2 / (11) / -0.546 / -0.433 / = -1/4
= 1/4 / (11) / 0.335 / 0.216
(12) / 0.45 / 0.433 / (12) / -0.26 / -0.216
(21) / -0.45 / -0.433 / (21) / 0.26 / 0.216
(22) / 0.546 / 0.433 / (22) / -0.335 / -0.216
1
It should be observed that an increase in the accuracy of the solution is connected with a decrease of the step h in accordance with (17) and (18).
1
1
1
References:
[1] FORRAY M.J., Variational Calculus in Science and
Engineering, Bucharest,Editura Tehnică,1975.
[2] LASCAUX P.,THÉODOR R., Analyse numérique
matricielle appliquée à l’art de l’ingénieur, Paris,
Masson, 1994.
[3] MARCIUK G.I., Numerical Methods,
Romanian Academia Publishing House,
Bucharest, 1983.
[4] MARCIUK G.I., LEBEDEV V.I., Numerical methods
in the neutron transport theory , Moscova, Atomizdat,
1971.
[5] MARTIN O., Une méthode de résolution de l’équation du
transfert des neutrons, Rev. Roum. Sci. Techn.-
Méc.Apl., 37,6(1992).
[6] ROMBALDI J.E., Problèmes corrigés d’analyse
numérique, Paris, Masson, 1996.
1