- Email: [email protected]
On the numerical solution of Landau's kinetic equation
121
Solution of Landau’s kinetic equation
b) the use of analytic a priori integration of the coordinate functions and the existence of a priori information about the coefficients of the approximating polynomial; the laborious problem of numerically integrating the coordinate functions with the operator of the equation is thereby simplified; c) simplicity and the minimal information needed when specifying the computational model of a multi-zone reactor; d) the possibility of computing neutron diffusion in reactors with virtually arbitrary zone geometry without introducing an error due to the description of the zone boundaries. Translated by D. E. Brown REFERENCES
1.
MIKHLIN, S. G., Vartitional methods in mathematical physics (Variatsionnye
metody v matematicheskoi
fmike), Nauka, Moscow, 1970. 2.
KANG, C. M., and HANSEN, K. F., Finite element methods for reactor analysis, Nucl. Sci. Engng., 51, 456-495,1973.
3.
SLESAREV, I. S., and SIROTKIN, A. M., Use of variational functionals with moving boundaries for constructing approximate methods for solving the neutron transport equation, Zh. vjchisl. Mat. mat. Fiz., 11, No. 2,529-534, 1971.
4.
GONCHAROV,
5.
KANTOROVICH, L. V., and KRYLOV, V. I.,Approximate methods of higher analysisfPriblizhennye metody vysshego anal&a), Gostekhizdat, Moscow, 1962.
6;
KHARRIK, I. Yu., On the approximation of functions, vanishing on the boundary of a domain, by functions of special type,Matem. sb., 37 (79), No. 2,353-384, 1955.
7.
VLADIMIROV,
V. L., Theory of interpolation and approximation of functions (Teoriya interpolirovaniya i priblizheniya funktsii), Gostekhizdat, Moscow, 1954.
V. S., Equations of mathematical physics (Uravneniya matematicheskoi
fiziki), Nauka,
Moscow, 1971. 8.
SLESAREV, I. S., SIROTKIN, A. M., and KHROMOV, V. V., Variational composite methods for computing neutron distribution in reactors, Atomnaya energiya, 32, No. 1,53,1972.
ON THE NUMERICAL
SOLUTION
OF LANDAU’S
KINETIC
EQUATION*
A. V. BOBYLEV and V. A. CHUYANOV Moscow (Received 14 June 1974)
A COMPLETELY conservative difference scheme is constructed for Landau’s kinetic equation. The problem of the spatially homogeneous and velocity-isotropic relaxation of an electron plasma is computed. The problem of the numerical solution of the Landau (Fokker-Planck) equation has recently become very topical because of the need for quantitative results on charged particle dynamics in magnetic thermonuclear traps. Computations are of fundamental importance in such problems, since in essence, they have to decide whether the device concerned is of actual advantage [ 1] . This fact makes specially serious demands on the quality of the computational algorithms. In the present paper, taking the example of the problem of the spatially homogeneous and velocity-isotropic *Zh. vfchisl Mat. mat. F’iz., 16,2,407-416,
1976.
A. V. Bobylev und V. A. Chuyunov
122
relaxation of a single-component (electron) plasma, we reveal the value of the principle of complete conservatism [2,3] for the numerical solution of such problems. We show that, if this principle is not observed, false results are obtained for the later stage of the relaxation. We construct a completely conservative difference scheme, which enables computations to be performed for long time intervals without accumulation of errors. We carry out a computation of the relaxation whereby the results of [4,5] (see also Section 7.8 of [6]), in which non-conservative schemes were used, can be corrected and augmented.
1. The problem Landau’s equation for a system of similarly charged particles (electrons) in the spatially homogeneous and velocity-isotropic case has the form dg(u,t)
=
at
Ifh2e4L mV
a dv 11s
g(s)24dx 0
(1.1) +u3 jg(l)sdz] t)
&Z+g(v,
t) jg(5)Z2dr}) 0
where g(v, t) is the velocity distribution function of the particles, and L is the Coulomb logarithm. We consider for Eq. (1 .l) the Cauchy problem on the half-line 0 < Y< 00for t > 0,assuming that g(v, t) is bounded at Y= 0, and the initial condition
g(u, O)=goW.
(1.2)
It is assumed below that the (classical) solution of this problem, possessing finite moments of second and fourth orders (signifying the number of particles and the energy of the system), exists and is unique. We put m
m
n=4n
s
go (0) U2au,
E=~nT=2nmSg,(u)u4dv.
0
We introduce the dimensionless variables x, T, and the fimctionflx,
(1.3)
0
7):
(l-4)
The equation for Ax, r) how has the form
(1.5)
123
Solution of Landau’s kinetic equation
Further, we shah consider this equation with the initial condition
fw~=foo=~( $)
-=/2
(1.6)
go(v).
Here we denote the variable r by c and call it time, while we call x the velocity, and the dimensionless function fix, t) the distribution function. The moments of the distribution function
J0
n = ;(X,
E=Sf(X,t)X”dX 0
t)z” dz,
willbe called respectively the “number of particles” and the “energy.” We introduce the notation
(1.7)
Equation (1.5) can be written in the equivalent forms
af at
+
=
$;
a
I
-=-_-x2
x4 ax
1[
(1.9)
(
(1.10)
x2
D(fJg+F(f.S)f]
JOY YG(x, Y> -g-+$)fkt)f(Y,t),
G (x, y) =mind(Y,
y”).
Equations (1.8) and (1.9) respectively express the laws of conservation of the number of particles, and of the energy of the system. Let us note the most important properties of the Cauchy problem (1.5), (1.6). Property
A. The integral laws of conservation
J (x, t)x” f
dx =
0
Jf(x, 0
t)x”dx=
J
fo(x)x’dx=rz,
0
Jf,(x)x4dx=E. 0
124
A. V. Bobylev and V. A. Chuyanov
are satisfied. In accordance with Eqs. (1.3) and (1.4), n = E = 1, but it is convenient for the present not to utilize this fact. &operp B. The problem has a unique stationary solution of the “Maxwell” type:
fM(5)=
-$($)
%exp (-gx2).
(1.11)
Here, any function of the type cp(z) =a exp (--by’) , where CIand 6 are arbitrary constants (b > 0), is a stationary solution of Eq. (1.5). The laws of conservation select, from this two-parameter set of functions, the unique equilibrium solution (1.1 I), corresponding to the initial condition (1.6). Z+operry C. We have the H-theorem:
dH(t) d = -=dtJf(Z)t)[lnf(x,t)]x2dxi0. dt 0
(1.12)
The proof is obvious, if we utilize the equation in the form (1.10). The minimum of the functional His reached (for given E and n) on the function&(x).
2. Completely conservative difference scheme The principle of complete conservatism, put forward in [2,3] in connection with problems of gas dynamics and magnetohydrodynamics is of unusual importance for the present problem. In fact, if we use Eqs. (1.8) and (1.9), expressing in explicit form the two different laws of conservation, we can easily construct, by the familiar integro-interpolation method [7], on a symmetric pattern, two difference schemes of order 0 (z+h’) (T and h are the mesh steps in time and velocity space respectively), to which there correspond the difference analogues of the two conservation laws. However, as distinct from the equivalent equations (1.8) and (1.9), the corresponding difference equations are not in general equivalent. It is easy to imagine what this can lead to. If, as the basis of the const~ction of the scheme, we take Eq. (1.8), then the difference analogue of the law of conservation of the number of particles is satisfied, but the energy is conserved only up to quantities sof order o(@). It can easily be shown that the energy in fact decreases monotonically. For long times, when the exact solution of the initial problem is close to the equilibrium solution and changes very slowly, in the case of the difference problem a deftitive role starts to be played by the dissipative properties of the scheme. As a result, all the information about the time-asymptotic behaviour of the exact solution, which can be extracted from the results of numerical computation of relaxation (and which is of greatest interest for the problem in question), proves to be false. Similarly, if, as the basis of the construction of the scheme we take Eq. (1.9), we arrive at a scheme for which the law of conservation of energy is satisfied, but the law of conservation of the number of particles is violated. Here, as distinct, for example, from the equations of gas dynamics [ZJ, the unbalance of the energy (or the number of particles), turns out to be ~dependent of the time step. The reason is that the conserved quantities mand E are here the scalar products of the solution &x, t) with
Solution of Landau’s kinetic equation
125
functions which are time-independent (x2 and x4 respectively). Hence we shall try to arrange for complete ~o~e~t~rn by utilizing only mesh approbation with respect to x of the operator on the ~~t-h~d side of Eqs. (1.8) and (1.9). For this, we write Eq. (1 S) in the form 1 d 1 aW(f,s) af -=:--x2 ax x ax at
’
(2.1)
where
(2.2) P (f,
x) = si(Y, t> y 4/. =
This way of writing the equation takes the most explicit account of the existence of two laws of conservation, and hence is convenient for the const~ction of our scheme. in the domain of space considered x, t (OS&S& t&O) we introduce the difference mesh {x~, t3, ~R,,=x,,+hk+~, k=O, 1,. . . , N-f; x0=0, xN=L; Ei+‘=tj+++‘, j=O> 1,. . . ; t,=O. For simplicity, we shall assume that the mesh is uniform (hh=h=const, rj=r=const). We define on the mesh the mesh functions gkj, Wk’(y) , PA’(#), corresponding to the functions f@V t), W(f, a!), P(f, 2). As usual, yd is the value of the mesh function at the base-point (xk, d), The superscript will as a rule be omitted, and we shall also use the notation with neither subscript nor superscript [7] : ykj=g,
y:‘-‘=
Y-Y -
y’,
= yt,
z
f Y+rYki ---_ YX? h
Yn’-Yk-i h
= Yjr,
y=o.5 (Y~j+Y~-*), E=&b-I,*=0.5 (&:Rf&--1).
Applying the integro-interpolation method to Eq. (2.1), we obtain the (implicit) difference scheme
(2.3) The difference equation (2.3) can be reduced by simple transformations to a form expressing the difference analogue of the law of conservation of energy: I
Yi
(
x--g- -
i?
L
wS-2’i;ij:
1x
(2.4)
The scheme constructed proves to be completely conservative, since Eqs. (2.3) and (2.4), expressing the two different laws of ~onse~ation, are ~gebr~c~y equivalent. It still needs to be
126
A. K Bobyiev and K A. Chuyanov
proved that the integral laws of conservation of the number of particles and the energy are satisfied ~ou~out the mesh domain. It follows from Eqs. (2.3) and (2.4) that ~tisfa~ion of the integral laws of con~rvation is connected with the boundary conditions. For, the initial Cauchy problem is posed on the half-line 0 6 x < m, and the corresponding difference problem in the interval OGTS&=L. The order of decrease of the distribution function for large x can be est~t~d from the Maxwell ~st~bution (1 .I 1). The boundary value x = L has to be suf~ciently large, in order to prevent the neglect of quantities of the order of the ~t~bution function for x = L by comparison with unity from destroying the required accuracy. For instance, with E = n = 1 and L = 5, we have, in accordance with (1 .I l), fM(L) -1W’. Hence the bo~&~ condition at x = L is not important. For brevity, we shall here only consider the simplest method of stating the bound monition. In fact, we put ye = 0. Then, on appro~mat~g the integrals in (2.2) according to the trapezoid formula, we get N--1
h
pk(y)=pk=--xkyk+ 2
E
hkyi, k=O,1,..,, N-l, pN=O,
i=k+I
According to the last expression, Wa =Wf=& which enables us to dispense with a boundary condition at x = 0, since yo does not appear in the difference equations, The difference analogues of the energy and number of particles of the system are
E
=
z k-i
N--l
hxklyk,
n=
x hx,,2yk. k=l
Using (2.3) and (2_4), we can easily show that the variation of E and 12over a time step is proportional to v#_~ I i.e., ex~onent~y small. Actually, this variation is small compared with the rounding errors in computer evaluations. Notice that, by stating a more complicated boundary condition at x = L, we can arrange for the laws of conservation to be satisfied exactly throughout the mesh domain, For the difference equation (2.3), with such a boundary condition and the initial condition yk’=fo (a), k=& 2, N1, it is possible to prove that the difference analogues of Properties A-and B of the initial . ..( problem are satisfied, and it is possible to find an equilibrium solution, approximating the Maxwell ~t~bution (1 .l 1).However, in view of the exponential decrease of the distribution function for large x, the results of computations t~ou~out the mesh domain, except for a small nei~bourh~d of the boundary x = L, prove to be virtually the same for any sensible choice of boundary condition at x = L. The scheme is constructed on a symmetric pattern and has the order of approximation 0 (r+h*) , if no account is taken of the error arising on the passage from the half-line to the interval; this error is exponentially small for a rapidly decreasing distribution function. The scheme is implicit and hence unconditionally stable. Its diverence properties are unaltered if y is replaced on the right-hand side of (2.3) by y (“)=oy+ (1-o) y, O=GcrG1. In particular, on putting u = 0.5, we can obtain a scheme of order 0 (za+h2), The extension of the proposed scheme to the case of a nonuniform mesh is obvious.
Solution of imdau’s kinetic equation
127
Returning to the notation with subscripts and superscripts, we write the scheme in a form convenient for computations:
(2.5)
Here we have used the notation
N-i h pkj(y) = y Xt,ykj+ C
hX,yij,
ic=l,
2, . * . 2 N-1,
ia=&i
Dt,,(y) -R/f (y) -0. The initial and boundary conditions are yk*=fo(xk),
k=f,
2,.. .,N--3,
y,j=o,
j=d,
2,. I..
The difference scheme (2.5) approximates the quasi-linear equation of parabolic type (1.8). Hence the solution of our problem is completely similar to the numerical solution of an ordinary quasi-linear equation of heat condition (see e.g., [7], p. 214) in accordance with an implicit difference scheme. At each time step, the non-linear difference equations are solved by the iterative method, while at each iteration the equations are solved by the pivotal condensation method. It is easily shown that the pivotal condensation is correctly set up provided that h<2E/3nL.
3. Numerical study of rekxation The most typical initial condition for the present problem is the monoenergetic ~stribution
f*(x)=6(5-I).
(3.1)
A problem with similar initial conditions was considered in [4,5] (see also Section 7.8 of [6], where the results of [4] are quoted). In these papers, the computations are performed according to nonconservative schemes, and hence the results are only valid for small times (of the order of unity), For the numerical solution of the problem according to our scheme, we took the interval OGzG5, L-=5. The velocity step was fixed at h = 0.05, while the time step had two values: r = 0.01 for t Q 5, and r = 0.05 for t > 5. The function (3.1) was approximated on the mesh in the usual way, i.e., y&) = l//z for x = 1, and y&) = 0 for x # 1. To monitor the operation of the scheme, we used the &theorem (1.12). Since the equably solution of the ~fference problem differs somewhat from the Maxwell ~tribution (1 .I 1), it cannot be expected that the H-theorem will be satisfied exactly. In fact, during the computations the U-function was decreased until tz25, then slightly
128
A. V. Bobylev and V. A. Chuyanov
increased (in such a way that the fust five places after the decimal point remained unchanged) and was stabilized at a value corresponding to the equilibrium solution of the difference problem. For comparison, computations were also performed using nonconservative schemes. In this case the H-function increases monotonically, starting from iFi+2. Due to its exponential decrease for large x, the function fix, t) is unsuitable for graphical representation. We therefore give below curves of the function g(x, t), defined in such a way that
f(? t>=fna(x)[l+&?(Tt) I,
(3.2)
where, in accordance with the initial condition (3.l),r
The calculation shows that the relaxation proceeds, roughly speaking, in two stages. In the first (fast) stage, damping of the initial conditions occurs, the distribution peak moves to the point x = 0, and the function j(x, t) takes a shape close to the Maxwell distribution. The duration of this stage is of the order of unity. Hence the characteristic time
taken as unity in (1.4) actually characterizes the duration of the first stage of relaxation. Our results, for times of the order of unity, are the same as the results of [4,.5] , so that we shall not dwell on them.
FIG. 1. Profile of the function g(x, f) a) for t = 0.01, b) for t=0.2,c)fort=1.6,d)fort=lO,e)fort=16,f)fort=24
Sa~utionof Landau‘s kinetic equation
129
The second (slow) stage of relaxation is more interesting; in this, the system has almost “forgotten” about the initial conditions, and passes to equilibrium. Here, the relaxation proceeds differently in the domain of thermal velocities (0 < x < 2) and in the high-energy domain (x > 2). In the thermal velocity range, in a time of the order of unity, a stable profile of the function g(x, t) is established (see Fig. 1); the function then damps in time, while hardly changing its shape. For t > IO, the damping is close to exponential. Notice that, for the chosen initial condition (3.1) the low-energy domain is under-populated with particles for long times (as compared with the equ~brium distribution). The reverse result, obtained in [4] by using a nonconservative scheme, is a result of the dissipative properties of the scheme. computations using nonconservative schemes tend to break off at comparatively short times (t < lo), since it becomes obvious from then on that the numerical solution not only does not approach the Maxwell equilibrium solution, but moves away from it.
FIG. 2 In the course of time, the maximum of gfx, t) moves into the high energy region (Fig. 2). The relaxation of the tail of the distribution takes place fairly slowly. The solution in this region recalls a travelling wave (see Fig. 2). If, for characterizing the relaxation of the tail, we use the speed of the front of this wave (xf is defined by the condition g(s*, t) ==-0.5 ), we obtain the approximate dependence
axa
1
dt_=501_1
x~s~3t
+ const.
The tail relaxation can thus be characterized by the time
tl ,
(3.3)
which depends on the velocity:
?TZ2V3
,,=,,+
12ne"Ln
‘
A similar tail relaxation time was introduced in [4] . The behaviour of the solution for large x is easily explained as follows. As x + 00,the Landau equation (1 S) can be replaced approximately by a linear equation, which, in the light of the initial condition (3.1), has the form df
-=-_ dt
1
d
x2 ds (
;;+i,.
130
A. V. Bobyiev and V. A. Chuyanov
For the function g(x, t) of (3.2) we now obtain
d&T 1 dg dt+--=z_-_. x2 dx
1
d
1 dg
3x
dx
x dx
If the right-hand side of this equation is regarded as a viscous term, then (3.3) defmes the characteristics of the ~o~espond~g hyperbolic equation (with zero viscosity). The computations results of Fig. 2 thus show that, for comparatively low velocities (x < 5) admitting numerical investigation, the viscosity does not seriously affect the velocity of wave propagation. Notice that, though we have here considered, for brevity, only one version of the problem with initial condition (3.1), the regular features of the slow stage of relaxation, indicated above (stability of the profile of the function g(x, t) in the domain of thermal velocities and the wave nature of the propagation of the solution in the high velocity domain) prove to be typical for other finite initial conditions. The authors thank A. A. Sarmaskii for useful discussions. ~Q~s&~e~by D. E. Brown REFERENCES 1.
KILLING, J., and MARKS, K. D., Solution of the Fokker-Planck equation for a plasma in a trap with magnetic probes, in: Computational methods in plasm4 physics (Vychisl. metody v fii. plazmy), Mir, Moscow, 417-482,1974. (Russian translation).
2.
POPOV, Yu. P., and SAMARSKII, A. A., Completely conservative difference schemes, .?h. vj%hisl.Mat. m4t. Fiz., 9, No. 4,953-958,1969.
3.
POPOV, Yu. P., and SAMARSKII, A. A., Completely conservative difference schemes for the ~~etohydr~y~i~, Zh. vj%h&f.b44t. mat. FEZ., 10,No. 4,990-998,197O.
4.
MacDONALD, W. M., ROSENBL~H, M. N., and CHUCK, W., Relaxation of a system of particles with coulomb interactions,Phys. Rev., 107, No. 2,350-353,1957.
5.
DOLINSKY, A., Numerical integration of kinetic equations, Phys. F&ids, 8, No. 3,436-443,
6.
SHKAROFSKY, I., JOHNSTON, T., and BACHINSKY, M., Particle kinetics of pkwws, 1966.
7.
SAMARSKII, A. A., Introduction to the theory of difference schemes (Vvedenie v teoriyu skhem), Nauka, Moscow, 1970.
equationsof
1965.
Addison-Wesley, raznostnykh
Solution of Landau’s kinetic equation
b) the use of analytic a priori integration of the coordinate functions and the existence of a priori information about the coefficients of the approximating polynomial; the laborious problem of numerically integrating the coordinate functions with the operator of the equation is thereby simplified; c) simplicity and the minimal information needed when specifying the computational model of a multi-zone reactor; d) the possibility of computing neutron diffusion in reactors with virtually arbitrary zone geometry without introducing an error due to the description of the zone boundaries. Translated by D. E. Brown REFERENCES
1.
MIKHLIN, S. G., Vartitional methods in mathematical physics (Variatsionnye
metody v matematicheskoi
fmike), Nauka, Moscow, 1970. 2.
KANG, C. M., and HANSEN, K. F., Finite element methods for reactor analysis, Nucl. Sci. Engng., 51, 456-495,1973.
3.
SLESAREV, I. S., and SIROTKIN, A. M., Use of variational functionals with moving boundaries for constructing approximate methods for solving the neutron transport equation, Zh. vjchisl. Mat. mat. Fiz., 11, No. 2,529-534, 1971.
4.
GONCHAROV,
5.
KANTOROVICH, L. V., and KRYLOV, V. I.,Approximate methods of higher analysisfPriblizhennye metody vysshego anal&a), Gostekhizdat, Moscow, 1962.
6;
KHARRIK, I. Yu., On the approximation of functions, vanishing on the boundary of a domain, by functions of special type,Matem. sb., 37 (79), No. 2,353-384, 1955.
7.
VLADIMIROV,
V. L., Theory of interpolation and approximation of functions (Teoriya interpolirovaniya i priblizheniya funktsii), Gostekhizdat, Moscow, 1954.
V. S., Equations of mathematical physics (Uravneniya matematicheskoi
fiziki), Nauka,
Moscow, 1971. 8.
SLESAREV, I. S., SIROTKIN, A. M., and KHROMOV, V. V., Variational composite methods for computing neutron distribution in reactors, Atomnaya energiya, 32, No. 1,53,1972.
ON THE NUMERICAL
SOLUTION
OF LANDAU’S
KINETIC
EQUATION*
A. V. BOBYLEV and V. A. CHUYANOV Moscow (Received 14 June 1974)
A COMPLETELY conservative difference scheme is constructed for Landau’s kinetic equation. The problem of the spatially homogeneous and velocity-isotropic relaxation of an electron plasma is computed. The problem of the numerical solution of the Landau (Fokker-Planck) equation has recently become very topical because of the need for quantitative results on charged particle dynamics in magnetic thermonuclear traps. Computations are of fundamental importance in such problems, since in essence, they have to decide whether the device concerned is of actual advantage [ 1] . This fact makes specially serious demands on the quality of the computational algorithms. In the present paper, taking the example of the problem of the spatially homogeneous and velocity-isotropic *Zh. vfchisl Mat. mat. F’iz., 16,2,407-416,
1976.
A. V. Bobylev und V. A. Chuyunov
122
relaxation of a single-component (electron) plasma, we reveal the value of the principle of complete conservatism [2,3] for the numerical solution of such problems. We show that, if this principle is not observed, false results are obtained for the later stage of the relaxation. We construct a completely conservative difference scheme, which enables computations to be performed for long time intervals without accumulation of errors. We carry out a computation of the relaxation whereby the results of [4,5] (see also Section 7.8 of [6]), in which non-conservative schemes were used, can be corrected and augmented.
1. The problem Landau’s equation for a system of similarly charged particles (electrons) in the spatially homogeneous and velocity-isotropic case has the form dg(u,t)
=
at
Ifh2e4L mV
a dv 11s
g(s)24dx 0
(1.1) +u3 jg(l)sdz] t)
&Z+g(v,
t) jg(5)Z2dr}) 0
where g(v, t) is the velocity distribution function of the particles, and L is the Coulomb logarithm. We consider for Eq. (1 .l) the Cauchy problem on the half-line 0 < Y< 00for t > 0,assuming that g(v, t) is bounded at Y= 0, and the initial condition
g(u, O)=goW.
(1.2)
It is assumed below that the (classical) solution of this problem, possessing finite moments of second and fourth orders (signifying the number of particles and the energy of the system), exists and is unique. We put m
m
n=4n
s
go (0) U2au,
E=~nT=2nmSg,(u)u4dv.
0
We introduce the dimensionless variables x, T, and the fimctionflx,
(1.3)
0
7):
(l-4)
The equation for Ax, r) how has the form
(1.5)
123
Solution of Landau’s kinetic equation
Further, we shah consider this equation with the initial condition
fw~=foo=~( $)
-=/2
(1.6)
go(v).
Here we denote the variable r by c and call it time, while we call x the velocity, and the dimensionless function fix, t) the distribution function. The moments of the distribution function
J0
n = ;(X,
E=Sf(X,t)X”dX 0
t)z” dz,
willbe called respectively the “number of particles” and the “energy.” We introduce the notation
(1.7)
Equation (1.5) can be written in the equivalent forms
af at
+
=
$;
a
I
-=-_-x2
x4 ax
1[
(1.9)
(
(1.10)
x2
D(fJg+F(f.S)f]
JOY YG(x, Y> -g-+$)fkt)f(Y,t),
G (x, y) =mind(Y,
y”).
Equations (1.8) and (1.9) respectively express the laws of conservation of the number of particles, and of the energy of the system. Let us note the most important properties of the Cauchy problem (1.5), (1.6). Property
A. The integral laws of conservation
J (x, t)x” f
dx =
0
Jf(x, 0
t)x”dx=
J
fo(x)x’dx=rz,
0
Jf,(x)x4dx=E. 0
124
A. V. Bobylev and V. A. Chuyanov
are satisfied. In accordance with Eqs. (1.3) and (1.4), n = E = 1, but it is convenient for the present not to utilize this fact. &operp B. The problem has a unique stationary solution of the “Maxwell” type:
fM(5)=
-$($)
%exp (-gx2).
(1.11)
Here, any function of the type cp(z) =a exp (--by’) , where CIand 6 are arbitrary constants (b > 0), is a stationary solution of Eq. (1.5). The laws of conservation select, from this two-parameter set of functions, the unique equilibrium solution (1.1 I), corresponding to the initial condition (1.6). Z+operry C. We have the H-theorem:
dH(t) d = -=dtJf(Z)t)[lnf(x,t)]x2dxi0. dt 0
(1.12)
The proof is obvious, if we utilize the equation in the form (1.10). The minimum of the functional His reached (for given E and n) on the function&(x).
2. Completely conservative difference scheme The principle of complete conservatism, put forward in [2,3] in connection with problems of gas dynamics and magnetohydrodynamics is of unusual importance for the present problem. In fact, if we use Eqs. (1.8) and (1.9), expressing in explicit form the two different laws of conservation, we can easily construct, by the familiar integro-interpolation method [7], on a symmetric pattern, two difference schemes of order 0 (z+h’) (T and h are the mesh steps in time and velocity space respectively), to which there correspond the difference analogues of the two conservation laws. However, as distinct from the equivalent equations (1.8) and (1.9), the corresponding difference equations are not in general equivalent. It is easy to imagine what this can lead to. If, as the basis of the const~ction of the scheme, we take Eq. (1.8), then the difference analogue of the law of conservation of the number of particles is satisfied, but the energy is conserved only up to quantities sof order o(@). It can easily be shown that the energy in fact decreases monotonically. For long times, when the exact solution of the initial problem is close to the equilibrium solution and changes very slowly, in the case of the difference problem a deftitive role starts to be played by the dissipative properties of the scheme. As a result, all the information about the time-asymptotic behaviour of the exact solution, which can be extracted from the results of numerical computation of relaxation (and which is of greatest interest for the problem in question), proves to be false. Similarly, if, as the basis of the construction of the scheme we take Eq. (1.9), we arrive at a scheme for which the law of conservation of energy is satisfied, but the law of conservation of the number of particles is violated. Here, as distinct, for example, from the equations of gas dynamics [ZJ, the unbalance of the energy (or the number of particles), turns out to be ~dependent of the time step. The reason is that the conserved quantities mand E are here the scalar products of the solution &x, t) with
Solution of Landau’s kinetic equation
125
functions which are time-independent (x2 and x4 respectively). Hence we shall try to arrange for complete ~o~e~t~rn by utilizing only mesh approbation with respect to x of the operator on the ~~t-h~d side of Eqs. (1.8) and (1.9). For this, we write Eq. (1 S) in the form 1 d 1 aW(f,s) af -=:--x2 ax x ax at
’
(2.1)
where
(2.2) P (f,
x) = si(Y, t> y 4/. =
This way of writing the equation takes the most explicit account of the existence of two laws of conservation, and hence is convenient for the const~ction of our scheme. in the domain of space considered x, t (OS&S& t&O) we introduce the difference mesh {x~, t3, ~R,,=x,,+hk+~, k=O, 1,. . . , N-f; x0=0, xN=L; Ei+‘=tj+++‘, j=O> 1,. . . ; t,=O. For simplicity, we shall assume that the mesh is uniform (hh=h=const, rj=r=const). We define on the mesh the mesh functions gkj, Wk’(y) , PA’(#), corresponding to the functions f@V t), W(f, a!), P(f, 2). As usual, yd is the value of the mesh function at the base-point (xk, d), The superscript will as a rule be omitted, and we shall also use the notation with neither subscript nor superscript [7] : ykj=g,
y:‘-‘=
Y-Y -
y’,
= yt,
z
f Y+rYki ---_ YX? h
Yn’-Yk-i h
= Yjr,
y=o.5 (Y~j+Y~-*), E=&b-I,*=0.5 (&:Rf&--1).
Applying the integro-interpolation method to Eq. (2.1), we obtain the (implicit) difference scheme
(2.3) The difference equation (2.3) can be reduced by simple transformations to a form expressing the difference analogue of the law of conservation of energy: I
Yi
(
x--g- -
i?
L
wS-2’i;ij:
1x
(2.4)
The scheme constructed proves to be completely conservative, since Eqs. (2.3) and (2.4), expressing the two different laws of ~onse~ation, are ~gebr~c~y equivalent. It still needs to be
126
A. K Bobyiev and K A. Chuyanov
proved that the integral laws of conservation of the number of particles and the energy are satisfied ~ou~out the mesh domain. It follows from Eqs. (2.3) and (2.4) that ~tisfa~ion of the integral laws of con~rvation is connected with the boundary conditions. For, the initial Cauchy problem is posed on the half-line 0 6 x < m, and the corresponding difference problem in the interval OGTS&=L. The order of decrease of the distribution function for large x can be est~t~d from the Maxwell ~st~bution (1 .I 1). The boundary value x = L has to be suf~ciently large, in order to prevent the neglect of quantities of the order of the ~t~bution function for x = L by comparison with unity from destroying the required accuracy. For instance, with E = n = 1 and L = 5, we have, in accordance with (1 .I l), fM(L) -1W’. Hence the bo~&~ condition at x = L is not important. For brevity, we shall here only consider the simplest method of stating the bound monition. In fact, we put ye = 0. Then, on appro~mat~g the integrals in (2.2) according to the trapezoid formula, we get N--1
h
pk(y)=pk=--xkyk+ 2
E
hkyi, k=O,1,..,, N-l, pN=O,
i=k+I
According to the last expression, Wa =Wf=& which enables us to dispense with a boundary condition at x = 0, since yo does not appear in the difference equations, The difference analogues of the energy and number of particles of the system are
E
=
z k-i
N--l
hxklyk,
n=
x hx,,2yk. k=l
Using (2.3) and (2_4), we can easily show that the variation of E and 12over a time step is proportional to v#_~ I i.e., ex~onent~y small. Actually, this variation is small compared with the rounding errors in computer evaluations. Notice that, by stating a more complicated boundary condition at x = L, we can arrange for the laws of conservation to be satisfied exactly throughout the mesh domain, For the difference equation (2.3), with such a boundary condition and the initial condition yk’=fo (a), k=& 2, N1, it is possible to prove that the difference analogues of Properties A-and B of the initial . ..( problem are satisfied, and it is possible to find an equilibrium solution, approximating the Maxwell ~t~bution (1 .l 1).However, in view of the exponential decrease of the distribution function for large x, the results of computations t~ou~out the mesh domain, except for a small nei~bourh~d of the boundary x = L, prove to be virtually the same for any sensible choice of boundary condition at x = L. The scheme is constructed on a symmetric pattern and has the order of approximation 0 (r+h*) , if no account is taken of the error arising on the passage from the half-line to the interval; this error is exponentially small for a rapidly decreasing distribution function. The scheme is implicit and hence unconditionally stable. Its diverence properties are unaltered if y is replaced on the right-hand side of (2.3) by y (“)=oy+ (1-o) y, O=GcrG1. In particular, on putting u = 0.5, we can obtain a scheme of order 0 (za+h2), The extension of the proposed scheme to the case of a nonuniform mesh is obvious.
Solution of imdau’s kinetic equation
127
Returning to the notation with subscripts and superscripts, we write the scheme in a form convenient for computations:
(2.5)
Here we have used the notation
N-i h pkj(y) = y Xt,ykj+ C
hX,yij,
ic=l,
2, . * . 2 N-1,
ia=&i
Dt,,(y) -R/f (y) -0. The initial and boundary conditions are yk*=fo(xk),
k=f,
2,.. .,N--3,
y,j=o,
j=d,
2,. I..
The difference scheme (2.5) approximates the quasi-linear equation of parabolic type (1.8). Hence the solution of our problem is completely similar to the numerical solution of an ordinary quasi-linear equation of heat condition (see e.g., [7], p. 214) in accordance with an implicit difference scheme. At each time step, the non-linear difference equations are solved by the iterative method, while at each iteration the equations are solved by the pivotal condensation method. It is easily shown that the pivotal condensation is correctly set up provided that h<2E/3nL.
3. Numerical study of rekxation The most typical initial condition for the present problem is the monoenergetic ~stribution
f*(x)=6(5-I).
(3.1)
A problem with similar initial conditions was considered in [4,5] (see also Section 7.8 of [6], where the results of [4] are quoted). In these papers, the computations are performed according to nonconservative schemes, and hence the results are only valid for small times (of the order of unity), For the numerical solution of the problem according to our scheme, we took the interval OGzG5, L-=5. The velocity step was fixed at h = 0.05, while the time step had two values: r = 0.01 for t Q 5, and r = 0.05 for t > 5. The function (3.1) was approximated on the mesh in the usual way, i.e., y&) = l//z for x = 1, and y&) = 0 for x # 1. To monitor the operation of the scheme, we used the &theorem (1.12). Since the equably solution of the ~fference problem differs somewhat from the Maxwell ~tribution (1 .I 1), it cannot be expected that the H-theorem will be satisfied exactly. In fact, during the computations the U-function was decreased until tz25, then slightly
128
A. V. Bobylev and V. A. Chuyanov
increased (in such a way that the fust five places after the decimal point remained unchanged) and was stabilized at a value corresponding to the equilibrium solution of the difference problem. For comparison, computations were also performed using nonconservative schemes. In this case the H-function increases monotonically, starting from iFi+2. Due to its exponential decrease for large x, the function fix, t) is unsuitable for graphical representation. We therefore give below curves of the function g(x, t), defined in such a way that
f(? t>=fna(x)[l+&?(Tt) I,
(3.2)
where, in accordance with the initial condition (3.l),r
The calculation shows that the relaxation proceeds, roughly speaking, in two stages. In the first (fast) stage, damping of the initial conditions occurs, the distribution peak moves to the point x = 0, and the function j(x, t) takes a shape close to the Maxwell distribution. The duration of this stage is of the order of unity. Hence the characteristic time
taken as unity in (1.4) actually characterizes the duration of the first stage of relaxation. Our results, for times of the order of unity, are the same as the results of [4,.5] , so that we shall not dwell on them.
FIG. 1. Profile of the function g(x, f) a) for t = 0.01, b) for t=0.2,c)fort=1.6,d)fort=lO,e)fort=16,f)fort=24
Sa~utionof Landau‘s kinetic equation
129
The second (slow) stage of relaxation is more interesting; in this, the system has almost “forgotten” about the initial conditions, and passes to equilibrium. Here, the relaxation proceeds differently in the domain of thermal velocities (0 < x < 2) and in the high-energy domain (x > 2). In the thermal velocity range, in a time of the order of unity, a stable profile of the function g(x, t) is established (see Fig. 1); the function then damps in time, while hardly changing its shape. For t > IO, the damping is close to exponential. Notice that, for the chosen initial condition (3.1) the low-energy domain is under-populated with particles for long times (as compared with the equ~brium distribution). The reverse result, obtained in [4] by using a nonconservative scheme, is a result of the dissipative properties of the scheme. computations using nonconservative schemes tend to break off at comparatively short times (t < lo), since it becomes obvious from then on that the numerical solution not only does not approach the Maxwell equilibrium solution, but moves away from it.
FIG. 2 In the course of time, the maximum of gfx, t) moves into the high energy region (Fig. 2). The relaxation of the tail of the distribution takes place fairly slowly. The solution in this region recalls a travelling wave (see Fig. 2). If, for characterizing the relaxation of the tail, we use the speed of the front of this wave (xf is defined by the condition g(s*, t) ==-0.5 ), we obtain the approximate dependence
axa
1
dt_=501_1
x~s~3t
+ const.
The tail relaxation can thus be characterized by the time
tl ,
(3.3)
which depends on the velocity:
?TZ2V3
,,=,,+
12ne"Ln
‘
A similar tail relaxation time was introduced in [4] . The behaviour of the solution for large x is easily explained as follows. As x + 00,the Landau equation (1 S) can be replaced approximately by a linear equation, which, in the light of the initial condition (3.1), has the form df
-=-_ dt
1
d
x2 ds (
;;+i,.
130
A. V. Bobyiev and V. A. Chuyanov
For the function g(x, t) of (3.2) we now obtain
d&T 1 dg dt+--=z_-_. x2 dx
1
d
1 dg
3x
dx
x dx
If the right-hand side of this equation is regarded as a viscous term, then (3.3) defmes the characteristics of the ~o~espond~g hyperbolic equation (with zero viscosity). The computations results of Fig. 2 thus show that, for comparatively low velocities (x < 5) admitting numerical investigation, the viscosity does not seriously affect the velocity of wave propagation. Notice that, though we have here considered, for brevity, only one version of the problem with initial condition (3.1), the regular features of the slow stage of relaxation, indicated above (stability of the profile of the function g(x, t) in the domain of thermal velocities and the wave nature of the propagation of the solution in the high velocity domain) prove to be typical for other finite initial conditions. The authors thank A. A. Sarmaskii for useful discussions. ~Q~s&~e~by D. E. Brown REFERENCES 1.
KILLING, J., and MARKS, K. D., Solution of the Fokker-Planck equation for a plasma in a trap with magnetic probes, in: Computational methods in plasm4 physics (Vychisl. metody v fii. plazmy), Mir, Moscow, 417-482,1974. (Russian translation).
2.
POPOV, Yu. P., and SAMARSKII, A. A., Completely conservative difference schemes, .?h. vj%hisl.Mat. m4t. Fiz., 9, No. 4,953-958,1969.
3.
POPOV, Yu. P., and SAMARSKII, A. A., Completely conservative difference schemes for the ~~etohydr~y~i~, Zh. vj%h&f.b44t. mat. FEZ., 10,No. 4,990-998,197O.
4.
MacDONALD, W. M., ROSENBL~H, M. N., and CHUCK, W., Relaxation of a system of particles with coulomb interactions,Phys. Rev., 107, No. 2,350-353,1957.
5.
DOLINSKY, A., Numerical integration of kinetic equations, Phys. F&ids, 8, No. 3,436-443,
6.
SHKAROFSKY, I., JOHNSTON, T., and BACHINSKY, M., Particle kinetics of pkwws, 1966.
7.
SAMARSKII, A. A., Introduction to the theory of difference schemes (Vvedenie v teoriyu skhem), Nauka, Moscow, 1970.
equationsof
1965.
Addison-Wesley, raznostnykh