Difference Scheme

for the Vlasov-Manev System

R. Illner and S. Rjasanow

Received: January 22, 1999

Communicated by Alfred K. Louis

Abstract. We develop and test a finite difference scheme for the Vlasov-Manev Equation in one space and one velocity dimension. The Manev correction to the Newtonian potential produces visible qualitative differences in the behaviour of stellar systems; the most notable effect observed in this paper is a stabilisation of the separate identities of two Maxwellian concentrations at different locations.

1991 Mathematics Subject Classification: 65N06,65N12,65N22,85A05 Keywords and Phrases: Vlasov-Manev system, Finite difference methods, stability

1 Introduction

We are concerned with the numerical solution of the Vlasov equation with a Manev-type correction to the potential in 1 + 1 (one space, one velocity) dimensions. In the three-dimensional case, the Newtonian potential is changed to the “Manev” potential

U(|x − y|) = −\frac{γ}{|x − y|} − \frac{δ}{|x − y|^{2}} .  (1)

The correction −δ/|x − y|^{2} was introduced by Manev in a series of papers [7],[8],[9], [10] in the 1920s in an attempt to find a semiclassical approximation to the relativistic central force problem. For γ the universal gravitational constant and δ = 3γ^{2}/c^{2}, where c is the speed of light, this correction gives a qualitatively accurate prediction of the precession of the perihelion of Mercury. Manev’s work was the main motivation for the recent paper [2], in which the authors discuss the properties of the corresponding stellar dynamic equation. Notably, it is shown that the Cauchy problem for this equation does not, in general, admit global solutions (the corresponding result for the classical stellar dynamic equation holds only if the number of space dimensions is larger than or equal to 4, see, e.g., [4]). This means that a stellar system driven by Manev forces will typically develop features where the spatial density ρ loses smoothness such that the Manev force term, i.e., the Riesz transform of ρ

E_{2}[ρ](t, x) := −δ ∫ \frac{x − y}{|x − y|^{4}} ρ(t, y) dy

diverges. Possible reasons for this are the local formation of singularities in ρ (“concentrations”) or in ∇_{x}ρ (e.g., “shock waves”). As it is well known that such singularities do not occur in solutions of the classical stellar dynamic equation (see [4]), the Manev correction may be of physical relevance in modelling the evolution of large stellar systems like galaxies, globular clusters or interstellar dust clouds. We remark that an equation with γ = 0 and δ > 0 (referred to as the “pure” Manev case in [2]) possesses an interesting “projective” invariance in addition to the standard translation, scaling and Galilei invariances. Specifically, as shown in [2], if f(t, x, v) is a solution of the pure Manev equation

∂_{t}f + v · ∇_{x}f + E_{2}[ρ] · ∇_{v}f = 0

and if for some a > 0, τ = t/(1 + at), y = x/(1 + at) and w = (1 + at)v − ax, then

F (τ, y, w) := f(t, x, v)

solves

∂_{τ} F + w · ∇_{y}F + E_{2}[ρ̃] · ∇_{w}F = 0

with ρ̃ = ∫ F dw. This “projective” invariance, described in the context of the corresponding N-body problem by Bobylev and Ibragimov [1], may be of significance (and use) in regions where ρ or ∇ρ are large and the Manev correction dominates the Newtonian forces. As also discussed in [2], Boltzmann collision terms are dimensionally of the same order as the Manev force term and should therefore be included in a proper model. However, the present study aims at the identification of effects which can be attributed to the Manev correction alone; we therefore omit Boltzmann collision terms and any other conceivable correction (such as, e.g., Smoluchowski type coagulation terms). It is tempting to try a particle or particle-in-cell scheme for this equation, as is common for the Vlasov-Poisson (VP) system. The most advanced scheme of this type for VP is due to Greengard and Rokhlin [5]. We experimented with particle schemes for the Vlasov-Manev equation and found that they performed poorly due to the strong singularity of the Manev correction at very short range. Specifically, particles could be accelerated to extreme velocities over one time step, an effect which can also happen for VP but which is rare enough to cause no difficulties. It should also be an insignificant effect for the Vlasov-Manev equation while the densities are smooth (this follows from the existence and uniqueness proof shown in [6]), so a good numerical method should reflect this; but particle methods do not. Of course, one could mollify the Manev potential in order to avoid the difficulty, but unless the mollification parameter is chosen with great care, this may obfuscate the effects of the Manev correction relative to the Newtonian forces. It is for these reasons that we decided to avoid particle methods altogether. The main objective of this paper is therefore the construction and testing of a difference scheme for the Vlasov-Manev (VM) equation in one-dimensional geometry; effective generalisation for multidimensional cases is planned for future work. Our paper is organised as follows. In Section 2, we describe the VM equation and summarise its properties. In Section 3, a difference scheme for VM is derived, and its properties are formulated and proved. Section 4 contains a few informative numerical examples.

2 The Vlasov-Manev equation and its properties

The one-dimensional Vlasov-Manev equation is associated with the potential

U(|x − y|) = −γ|x − y| − δ ln |x − y|.

This potential arises from (1) by assuming homogeneity of the stellar system in the y− and z− directions, and the Vlasov-Manev equation can be written in the form

f_{t} + v f_{x} + E f_{v} = 0,  (2)

where f = f(t, x, v) : ℝ_{+}×ℝ_{x}×ℝ_{v} → ℝ_{+} is a non-negative distribution density function, t ≥ 0 is the time variable, x is the space variable and v is the velocity variable. The force function E = E(t, x) is defined as follows:

E(t, x) = −γ ∫_{ℝ_{y}} \frac{x − y}{|x − y|}ρ(t, y)dy − δ ∫_{ℝ_{y}} \frac{x − y}{|x − y|^{2}} ρ(t, y)dy,  (3)
ρ(t, x) = ∫_{ℝ_{v}} f(t, x, v)dv.

or E = −∇U, with

U(t, x) = −γ ∫ |x − y|ρ(t, y) dy − δ ∫ ln |x − y|ρ(t, y) dy.

182 R. Illner and S. Rjasanow Here, ρ denotes the spatial density. The non-negative constants γ and δ are given, and typically γ >> δ. The equation (2) is complemented with the initial condition

f(0, x, v) = f_{0}(x, v) ≥ 0.  (4)

We summarise the main properties of this equation.

1. Conservation of non-negativity

if f_{0}(x, v) ≥ 0, then for t > 0 f(t, x, v) ≥ 0.

2. Conservation of mass

m(t) = ∫_{ℝ_{x}} ∫_{ℝ_{v}} f(t, x, v)dvdx = m(0) = ∫_{ℝ_{x}} ∫_{ℝ_{v}} f_{0}(x, v)dvdx.

3. The continuity equation reads

ρ_{t}(t, x) + j_{x}(t, x) = 0,  (5)
j(t, x) = ∫_{ℝ_{v}} v f(t, x, v)dv.

4. Conservation of energy

e(t) = \frac{1}{2} ∫_{ℝ_{v}} ∫_{ℝ_{x}} v^{2}f(t, x, v)dxdv
− \frac{1}{2} ∫_{ℝ_{x}} ∫_{ℝ_{y}}(γ|x − y| + δ ln |x − y|)ρ(t, x)ρ(t, y)dxdy = e(0).

5. Second derivative of the moment of inertia

\frac{d^{2}}{dt^{2}} ∫_{ℝ_{x}} x^{2}ρ(t, x)dx = 2 ∫_{ℝ_{x}} ∫_{ℝ_{v}} v^{2}f(t, x, v) dv dx
− γ ∫_{ℝ_{x}} ∫_{ℝ_{y}} |x − y|ρ(t, x)ρ(t, y)dy dx − δ m^{2}(0).

Documenta Mathematica 4 (1999) 179–201

In the three-dimensional case and for γ = 0, the last identity is

\frac{d^{2}}{dt^{2}} ∫ x^{2}ρ(t, x) dx = 4e(0),

and this can be used to show non-global existence (for γ = 0) whenever the initial energy is negative. This argument is not applicable in the one-dimensional case under consideration; it may well be that solutions in this situation always exist globally.

3 A difference scheme

We begin our numerical study of the initial value problem (2),(4) with the discretisation of the physical and velocity spaces. First, we restrict the whole space ℝ_{x} × ℝ_{v} to a rectangle

Q_{L} = {(x, v) ∈ ℝ_{x} × ℝ_{v}, −L_{x} ≤ x ≤ L_{x}, −L_{v} ≤ v ≤ L_{v}} ,

and assume that f(t, x, v) has its support with respect to x and v in the box Q_{L}. We can then compute the force field by integrating over Q_{L} alone; later, we shall act as if f is extended periodically in x and v. The next step is the discretisation of the rectangle Q_{L} using the nodes

(x_{i}, v_{j}) = (−L_{x} + i h_{x}, −L_{v} + j h_{v}), (i, j) ∈ Q_{n}
h_{x} = \frac{2L_{x}}{n_{x}} , n_{x} ∈ ℕ, and n_{x} is even,
h_{v} = \frac{2L_{v}}{n_{v}} , n_{v} ∈ ℕ, and n_{v} is even,
Q_{n} = {(i, j) ∈ ℤ^{2}, 0 ≤ i ≤ n_{x}, 0 ≤ j ≤ n_{v}}.

Furthermore, we introduce the index set Q̃_{n} as a subset of all vectors in Q_{n} excluding those which have the form (0, j) or (i, 0). Let τ > 0 be the time discretisation parameter, and t_{k} = k τ, k = 0, 1, …. The function f(t_{k}, x, v) will now be represented by a vector f^{k} ∈ ℝ^{n}, n = n_{x} n_{v} with components

f_{l}^{k} = f_{i,j}^{k} ≈ f(t_{k}, x_{i}, v_{j}), (i, j) ∈ Q̃_{n}.  (6)

Here, l denotes the global index of the vector f^{0}, given by

l = (j − 1) n_{x} + i, l = 1, …, n.

Note that (6) defines the vector f^{k} for all (i, j) ∈ Q_{n} by the assumed periodic extension of f. It is also convenient to use the matrix form of the unknown function:

F^{k} ∈ ℝ^{n_{x}×n_{v}}.

The numerical density ρ_{i}^{k} ≈ ρ(t_{k}, x_{i}) can be computed using the midpoint integration rule

ρ_{i}^{k} = h_{v} ∑_{j=1}^{n_{v}} f_{i,j}^{k}, i = 1, …, n_{x},
ρ_{0}^{k} = ρ_{n_{x}}^{k}, k = 0, 1, …,

or in the matrix form

ρ^{k} = h_{v}F^{k}e_{n_{v}},  (7)

where ρ^{k} ∈ ℝ^{n_{x}} and e_{n_{v}} = (1, …, 1)^{T} ∈ ℝ^{n_{v}}. The total mass of the system is

m^{k} = h_{x} ∑_{i=1}^{n_{x}} ρ_{i}^{k}

and can be computed as follows:

m^{k} = h_{x} (ρ^{k}, e_{n_{x}}) = h_{x} h_{v} (f^{k}, e_{n}) = h_{x} h_{v} (F^{k}e_{n_{v}}, e_{n_{x}}) .  (8)

The next and most involved step is the numerical computation of the force due to (3). We take advantage of the fact that it is sufficient to integrate over one spatial period in (3), because what is really done is treat the case where the support of ρ stays inside such a period. Using the notation

P (x − y) = −γ \frac{x − y}{|x − y|} − δ \frac{x − y}{|x − y|^{2}}

and the piecewise representation of the density ρ we compute

E(t_{k}, x_{i}) ≈ E_{i}^{k} = ∑_{j=1}^{n_{x}} ρ_{j}^{k}G_{ij},  (9)
E^{k} = G ρ^{k}, E^{k} ∈ ℝ^{n_{x}}, G ∈ ℝ^{n_{x}×n_{x}}.  (10)

The elements of the matrix G are defined by

G_{ij} = ∫_{x_{j}−h_{x}/2}^{x_{j}+h_{x}/2} P (x_{i} − y)dy.  (11)

Direct computation of the force via (9) will require O(n_{x}^{2}) arithmetical operations in each time step and is therefore an “expensive” step. The combination of the special form of the matrix G and uniform discretisation leads to a special, Toeplitz form of the matrix G.

Lemma 1 The matrix G defined in (11) is a skew-symmetric Toeplitz matrix.

Proof: A matrix G is Toeplitz if

G_{i+1,j+1} = G_{ij}, i, j = 1, …, n_{x} − 1.

The analytical integration in (11) leads to

G_{ij} = γ h_{x} − δ ln \frac{j − i − 1/2}{j − i + 1/2}, j > i.

For j < i we get

G_{ij} = −G_{ji}

because of (11). The matrix G is therefore uniquely defined by its first row. The element G_{11} is a strongly singular integral which should be considered as a Cauchy integral

G_{11} = lim_{ε→0} (∫_{x_{1}−h_{x}/2}^{x_{1}−ε} P (x_{1} − y)dy + ∫_{x_{1}+ε}^{x_{1}+h_{x}/2} P (x_{1} − y)dy) .

Using the substitution y = −y′ + 2x_{1} in the second integral and the obvious property

P (x − y) = −P (y − x)

we obtain

∫_{x_{1}+ε}^{x_{1}+h_{x}/2} P (x_{1} − y)dy = − ∫_{x_{1}−ε}^{x_{1}−h_{x}/2} P (x_{1} + y′ − 2x_{1})dy′
= ∫_{x_{1}−h_{x}/2}^{x_{1}−ε} P (y′ − x_{1})dy′

= − ∫_{x_{1}−h_{x}/2}^{x_{1}−ε} P (x_{1} − y′)dy′,

and therefore

G_{11} = 0.

Remark 1 The multiplication of a Toeplitz matrix with a vector can be realised efficiently using the following well known trick. The matrix G can be considered as a left-upper block of the circulant matrix G̃ of the dimension m which is a power of two:

G̃ = ( \begin{matrix}G & G_{12} \\ G_{21} & G_{22} \end{matrix} ) ,

The matrix G̃ has the following additional property

G̃_{i,m} = G̃_{i+1,1}, i = 1, …, m − 1

and its first row is defined as

(G_{11}, …, G_{1,n_{x}}, 0, …, 0, −G_{n_{x},1}, −G_{n_{x}−1,1}, …, −G_{2,1}) ∈ ℝ^{1×m}.  (12)

The dimension m of the matrix G̃ is the next power of two for the number 2 n_{x} − 1 and the number of zeros in (12) is equal to m − 2n_{x} + 1. It is obviously true that 2 n_{x} − 1 < m < 4 n_{x} − 4. Each circulant matrix C of the dimension m can be represented as

C = m^{−1}F_{m} Λ F_{m}^{∗},

where F_{m} denotes the matrix of the Discrete Fourier Transform (DFT) of the dimension m. The diagonal matrix Λ contains the eigenvalues of C and can be computed as

Λ = diag(λ_{1}, …, λ_{m}) = diag(F_{m}C^{T} e_{1}),  (13)

e.g. as the DFT of the first row of the matrix C.

The computation of the force due to (10) therefore reduces to one DFT in a preparatory step for the computation of the eigenvalues of the matrix G̃ using (13) (the matrix G̃ remains fixed during the time steps), and then to two DFTs and one multiplication of the diagonal matrix n_{x}^{−1}Λ with a vector. If m is a power of two then the computation of the force only requires O(n_{x} log_{2}(n_{x})) arithmetical operations using the Fast Fourier Transform (FFT) [3],[13] and is therefore much “cheaper” than the computation of the density via (7) which requires O(n) = O(n_{x} n_{v}) arithmetical operations. The number of arithmetical operations for the multiplication

G̃ρ̃^{k} = ( \begin{matrix}G & G_{12} \\ G_{21} & G_{22} \end{matrix} ) ( \begin{matrix}ρ_{k} \\ 0 \end{matrix} ) = ( \begin{matrix}G ρ_{k} \\ G_{21} ρ_{k} \end{matrix} )  (14)

is then of the same order O(n_{x} log_{2}(n_{x})).

Remark 2 In a two- (d = 2) or three-dimensional case (d = 3) we will obtain a circulant-block matrix [12]. Such matrices can again be efficiently multiplied with a vector using the FFT. The amount of arithmetical work would be O(n_{x}^{d} log_{2}(n_{x})) in this case.

The next step in the numerical procedure we are describing is the discretisation of the equation (2) using a semi-implicit difference scheme. At the given time level k we compute the density ρ^{k} via (7) and the force E^{k} via (10) or (14), and then we use the following “upwind” approximation for the derivatives in (2):

f_{t}(t_{k}, x_{i}, v_{j}) ≈ f_{t,ij}^{k} = \frac{f_{ij}^{k+1} − f_{ij}^{k}}{τ} ,  (15)

v f_{x}(t_{k}, x_{i}, v_{j}) ≈ v_{j} f_{x,ij}^{k+1} = \{ \begin{matrix}v_{j} \frac{f_{ij}^{k+1} − f_{i−1,j}^{k+1}}{h_{x}} & , & v_{j} ≥ 0 \\ v_{j} \frac{f_{i+1,j}^{k+1} − f_{ij}^{k+1}}{h_{x}} & , & v_{j} < 0 \end{matrix} , (16)

E f_{v}(t_{k}, x_{i}, v_{j}) ≈ E_{i}^{k} f_{v,ij}^{k+1} = \{ \begin{matrix}E_{i}^{k} \frac{f_{ij}^{k+1} − f_{i,j−1}^{k+1}}{h_{v}} & , & E_{i}^{k} ≥ 0 \\ E_{i}^{k} \frac{f_{i,j+1}^{k+1} − f_{ij}^{k+1}}{h_{v}} & , & E_{i}^{k} < 0 \end{matrix} . (17)

The resulting difference scheme can now be written in the form

f_{t,ij}^{k} + v_{j} f_{x,ij}^{k+1} + E_{i}^{k} f_{v,ij}^{k+1} = 0, i, j ∈ Q̃_{n}, k = 0, 1, ….  (18)

The initial values f_{ij}^{0} = f_{0}(x_{i}, v_{j}) are given. After multiplication by τ, (18) is a system of linear equations which can be written in the matrix form

A_{k} f^{k+1} = f^{k}, A_{k} ∈ ℝ^{n×n}, f^{k}, f^{k+1} ∈ ℝ^{n}, n = n_{x} n_{v}.

The elements of the matrix A_{k} using the global numbering with l = n_{x}(j−1)+i are of the form

(A_{k})_{ll} = 1 + \frac{τ}{h_{x}} |v_{j}| + \frac{τ}{h_{v}} |E_{i}^{k}|,  (19)
(A_{k})_{l,l−1} = \{ \begin{matrix}−\frac{τ}{h_{x}} |v_{j}| & , & v_{j} ≥ 0 \\ 0 & , & v_{j} < 0 \end{matrix} ,  (20)

(A_{k})_{l,l+1} = \{ \begin{matrix}0 & , & v_{j} ≥ 0 \\ −\frac{τ}{h_{x}} |v_{j}| & , & v_{j} < 0 \end{matrix} ,  (21)

(A_{k})_{l,l−n_{x}} = \{ \begin{matrix}−\frac{τ}{h_{v}} |E_{i}^{k}| & , & E_{i}^{k} ≥ 0 \\ 0 & , & E_{i}^{k} < 0 \end{matrix} ,  (22)

(A_{k})_{l,l+n_{x}} = \{ \begin{matrix}0 & , & E_{i}^{k} ≥ 0 \\ −\frac{τ}{h_{v}} |E_{i}^{k}| & , & E_{i}^{k} < 0 \end{matrix} ,  (23)
l = 1, …, n .

All other elements of the matrix A_{k} are equal to zero, i.e. the matrix A_{k} is extremely sparse. Exactly three elements of each row of this matrix are unequal to zero. In a d−dimensional case this number would be 2d + 1.

Remark 3 If the indices in (15)-(17) or in (20)-(24) are not from the set Q̃_{n} then we always assume the periodic property (e.g. f_{n_{x}+1,j}^{k+1} ≡ f_{1,j}^{k+1} etc.).

The main properties of the difference scheme (18) correspond to the properties of the matrices A_{k}, k = 0, ….

Lemma 2 The matrix A_{k} has the following properties

1. A_{k} is a regular M-matrix,

2. A_{k} e_{n} = e_{n}, A_{k}^{T} e_{n} = e_{n},

3. ∥A_{k}^{−1}∥_{2} = 1.

Here ∥A_{k}^{−1}∥_{2} denotes the spectral norm of the matrix A_{k}^{−1}, i.e. its biggest singular value.

Proof:

1. The elements of the matrix A_{k} fulfil

(A_{k})_{ll} > 0, (A_{k})_{lm} ≤ 0, l ≠ m,
∑_{m=1}^{n} (A_{k})_{lm} = 1, l = 1, …, n.  (24)

By (24), the matrix A_{k} is strongly diagonal-dominant and therefore a regular M−matrix.

2. The first property is given trivially in (24). This means that the vector e_{n} is an eigenvector of the matrix A_{k} and corresponds to the eigenvalue one. The matrix A_{k}^{T} has the same eigenvector and the same eigenvalue because

(A_{k}^{T} e_{n})_{l} = ∑_{m=1}^{n} (A_{k})_{ml} =
= (A_{k})_{ll} + (A_{k})_{l+1,l} + (A_{k})_{l−1,l} + (A_{k})_{l−n_{x},l} + (A_{k})_{l+n_{x},l} .

Using the representations (l + 1, l) = (l + 1, (l + 1) − 1) and (l − 1, l) = (l − 1, (l − 1) + 1), the property l ± 1 = (n_{x} − 1)j + (i ± 1) and (20),(21) we obtain

(A_{k})_{l+1,l} + (A_{k})_{l−1,l} = −\frac{τ}{h_{x}} |v_{j}|.

By analogy

(A_{k})_{l−n_{x},l} + (A_{k})_{l+n_{x},l} = −\frac{τ}{h_{v}} |E_{i}^{k}|.

Together with (19) we obtain the required result.

3. The matrix A_{k}^{−1} is element-wise non-negative, because it is the inverse of the M−matrix. The spectral norm of the matrix A_{k}^{−1} is equal to its largest singular value or to the square root of the largest eigenvalue of the matrix A_{k}^{−T} A_{k}^{−1}. This matrix only has non-negative elements (as a product of two element-wise non-negative matrices) and the real eigenvector e_{n} only has positive components. Then the corresponding eigenvalue (Perron-Frobenius theorem) is the largest Perron-eigenvalue of this matrix. In our case this eigenvalue is equal to one, because of the properties in 2. Hence the spectral norm of the matrix A_{k}^{−1} is equal to one.

More details concerning M−matrices can be found in [11]. The above lemma enables us to prove some important properties of the difference scheme

1. Initial step
f_{ij}^{0} = f_{0}(x_{i}, v_{j}),
2. Time step for k = 0, 1, …
2.1 ρ^{k} = h_{v} F^{k} e_{n_{v}},
2.2 E^{k} = G ρ^{k},
2.3 A_{k}f^{k+1} = f^{k}.

(25)

Corollary 1 The solution of the difference scheme (25) exists for all k = 0, 1 ….

Proof: This property follows directly from the regularity of the matrix A_{k} for all k = 0, 1 ….

Corollary 2 If the initial function f_{0}(x, v) is non-negative then the vectors f^{k} remain component-wise non-negative for all k = 0, 1, ….

Proof: The initial vector f^{0} is component-wise non-negative because of its definition in step 1 of (25). If f^{k}, k = 0, 1, … is component-wise non-negative, then we obtain from Step 2.3 of (25)

f^{k+1} = A_{k}^{−1}f^{k}.

The matrix A_{k}^{−1} is component-wise non-negative because it is inverse of an M-matrix. The proof is then completed by induction.

Corollary 3 The difference scheme (25) conserves mass.

Proof: The mass of the system can be computed for k = 1, 2, … corresponding to formula (8)

m^{k} = h_{x}h_{v}(f^{k}, e_{n}) = h_{x}h_{v}(A_{k−1}f^{k−1}, e_{n})
= h_{x}h_{v}(f^{k−1}, A_{k−1}^{T}e_{n}) = h_{x}h_{v}(f^{k−1}, e_{n}) = m^{k−1} = … = m^{0}.

Corollary 4 The difference scheme (25) is stable in the discrete maximum norm with respect to the initial data.

Proof: The discrete maximum norm of f^{k}, k = 0, 1, … is defined as

∥f^{k}∥_{∞} = max_{l} |f_{l}^{k}| = max_{l} f_{l}^{k} = f_{l^{∗}}^{k}

because the components of the vector f^{k} are non-negative. Here we have used the global numbering l = n_{x}(j − 1) + i of the components of f^{k}. Using

(A_{k−1})_{l^{∗}l^{∗}} > 0, (A_{k−1})_{ij} ≤ 0, i ≠ j

in the time step k = 1, 2, … we obtain the following estimate for the index l^{∗}

∥f^{k}∥_{∞} = f_{l^{∗}}^{k} = ((A_{k−1})_{l^{∗}l^{∗}} + (A_{k−1})_{l^{∗}l^{∗}−1}
+ (A_{k−1})_{l^{∗}l^{∗}+1} + (A_{k−1})_{l^{∗}l^{∗}−n_{x}} + (A_{k−1})_{l^{∗}l^{∗}+n_{x}}) f_{l^{∗}}^{k}
≤ (A_{k−1})_{l^{∗}l^{∗}}f_{l^{∗}}^{k} + (A_{k−1})_{l^{∗}l^{∗}−1}f_{l^{∗}−1}^{k}
+ (A_{k−1})_{l^{∗}l^{∗}+1}f_{l^{∗}+1}^{k} + (A_{k−1})_{l^{∗}l^{∗}−n_{x}}f_{l^{∗}−n_{x}}^{k} + (A_{k−1})_{l^{∗}l^{∗}+n_{x}}f_{l^{∗}+n_{x}}^{k}=
+ (A_{k−1}f^{k})_{l^{∗}} = f_{l^{∗}}^{k−1} ≤ ∥f^{k−1}∥_{∞} ≤ … ≤ ∥f^{0}∥_{∞}.

Corollary 5 The difference scheme (25) is stable in the discrete L_{2}−norm with respect to the initial data.

Proof: The discrete L_{2}−norm of f^{k}, k = 0, 1, … is defined as

∥f^{k}∥_{2}^{2} = h_{x}h_{v}(f^{k}, f^{k}).

Using this definition and property 3 in Lemma 2 we obtain for k = 1, …

∥f^{k}∥_{2} = ∥A_{k−1}^{−1}f^{k−1}∥_{2} ≤ ∥A_{k−1}^{−1}∥_{2}∥f^{k−1}∥_{2} = ∥f^{k−1}∥_{2} ≤ … ≤ ∥f^{0}∥_{2}.

Corollary 6 If the sequence {f^{k}} converges then it converges to the constant

lim_{k→∞} f^{k} = \frac{m^{0}}{4L_{x}L_{v}} e_{n}

Proof: If the sequence {f^{k}} converges to f^{∞} then this vector fulfils

A_{∞}f^{∞} = f^{∞},

where A_{∞} denotes the limit of the sequence of matrices {A_{k}}. Since the matrix A_{∞} is still a regular M−matrix, its eigenvalue 1 is simple. It means that only a constant vector f^{∞} = αe_{n} can fulfil the equation (26). The constant α can be obtained using the conservation of mass

m^{0} = h_{x}h_{v}(f^{∞}, e_{n}) = h_{x}h_{v}α(e_{n}, e_{n}) = α(h_{x}n_{x})(h_{v}n_{v}) = α(4L_{x}L_{v}).

Next, we obtain the discrete form of the continuity equation (5). We will use the following notations

v = (v_{1}, …, v_{n_{v}})^{T} ∈ ℝ^{n_{v}}  − vector of the velocities,
D_{v} = diag(v) ∈ ℝ^{n_{v}×n_{v}}  − corresponding diagonal matrix,
D_{v}^{+} = diag(0.5(|v| + v)) ∈ ℝ^{n_{v}×n_{v}} − positive part of D_{v},
D_{v}^{−} = diag(0.5(|v| − v)) ∈ ℝ^{n_{v}×n_{v}} − negative part of D_{v},
w = (v_{1}^{2}, …, v_{n_{v}}^{2})^{T} ∈ ℝ^{n_{v}}  − vector of squares of the velocities,
E^{k} = (E_{1}^{k}, …, E_{n_{x}}^{k})^{T} ∈ ℝ^{n_{x}} − vector of the forces, (∗)
D_{E} = diag(E^{k}) ∈ ℝ^{n_{x}×n_{x}}  − corresponding diagonal matrix,
D_{E}^{+} = diag(0.5(|E^{k}| + E^{k})) ∈ ℝ^{n_{x}×n_{x}} − positive part of D_{x},
D_{E}^{−} = diag(0.5(|E^{k}| − E^{k})) ∈ ℝ^{n_{x}×n_{x}} − negative part of D_{v},
J_{m} = circ(0, 1, 0, …, 0) ∈ ℝ^{m×m}  − circulant matrix of the dimension m,
ρ^{k} = h_{v} F^{k} e_{n_{v}} ∈ ℝ^{n_{x}}  − density,
j^{k} = h_{v} F^{k} v ∈ ℝ^{n_{x}}  − numerical flux.

Using (∗) we rewrite the difference scheme (18) in the matrix form

\frac{F^{k+1} − F^{k}}{τ} + \frac{1}{h_{x}} ((I_{n_{x}} − J_{n_{x}})F^{k+1}D_{v}^{−} + (I_{n_{x}} − J_{n_{x}}^{T})F^{k+1}D_{v}^{+})
+ \frac{1}{h_{v}} (D_{E}^{−}F^{k+1}(I_{n_{v}} − J_{n_{v}}^{T}) + D_{E}^{+}F^{k+1}(I_{n_{v}} − J_{n_{v}})) = 0,
k = 0, 1, ….

If we multiply this matrix with the vector h_{v}e_{n_{v}} then we obtain using

(I_{n_{v}} − J_{n_{v}})e_{n_{v}} = (I_{n_{v}} − J_{n_{v}}^{T})e_{n_{v}} = 0,
D_{v}^{−}e_{n_{v}} = v^{−},
D_{v}^{+}e_{n_{v}} = v^{+}

the following equation

\frac{ρ^{k+1} − ρ^{k}}{τ} + \frac{h_{v}}{h_{x}} ((I_{n_{x}} − J_{n_{x}})F^{k+1}v^{−} + (I_{n_{x}} − J_{n_{x}}^{T})F^{k+1}v^{+}) ,

\frac{ρ^{k+1} − ρ^{k}}{τ} + \frac{h_{v}}{h_{x}} (0.5(I_{n_{x}} − J_{n_{x}})F^{k+1}(|v| − v) + 0.5(I_{n_{x}} − J_{n_{x}}^{T})F^{k+1}(|v| + v)),
\frac{ρ^{k+1} − ρ^{k}}{τ} + \frac{1}{2h_{x}} (J_{n_{x}} − J_{n_{x}}^{T})j^{k+1} + \frac{1}{2}h_{x}h_{v} \frac{1}{h_{x}^{2}} (2I_{n_{x}} − J_{n_{x}} − J_{n_{x}}^{T})F^{k+1}|v| = 0

or

\frac{ρ^{k+1} − ρ^{k}}{τ} + \frac{1}{2h_{x}} (J_{n_{x}} − J_{n_{x}}^{T})j^{k+1} = −\frac{1}{2}h_{x}h_{v} \frac{1}{h_{x}^{2}} (2I_{n_{x}} − J_{n_{x}} − J_{n_{x}}^{T})F^{k+1}|v|.

The short form of this equation is

ρ_{t}^{k} + jẋ = −\frac{1}{2}h_{x}(h_{v}F^{k+1}|v|)_{xx},  (26)

where y_{x} denotes the central difference and y_{xx} the second difference of the grid function y. The equation (26) corresponds to the continuous equation (5). While the left hand side of (26) is a possible correct approximation of the derivatives in (5), the right hand side forms an artificial viscosity of our scheme. Because of this term which is of the order O(h_{x}) we are not able to obtain the conservation of the energy of the scheme directly. However, our numerical tests show that the variation of the energy in one time step is small.

4 Numerical examples

In this section we calculate some examples using our difference scheme. The initial distribution f_{0}(x, v) is given by

f_{0}(x, v) = \frac{1}{2π\sqrt{T_{x}T_{v}}} (exp (−\frac{(x − x_{0})^{2}}{2T_{x}} ) + exp (−\frac{(x + x_{0})^{2}}{2T_{x}} )) exp (−\frac{v^{2}}{2T_{v}} ) ,

where T_{x}, T_{v} and x_{0} are some positive parameters. In Figures 1,2 we present the time evolution of the density and of the force in the time interval (0, 1.4) for the following setting of parameters: γ = 4, δ = 0, T_{x} = 2, T_{v} = 0.05, x_{0} = 4 and n_{x} = 60, n_{v} = 90, i.e. for the pure Vlasov case. The time interval (0, 1.4) is sufficient to show the main numerical effects.

194 R. Illner and S. Rjasanow Figure 1: The density profiles for γ = 4, δ = 0 Figure 2: The force profiles for γ = 4, δ = 0

Documenta Mathematica 4 (1999) 179–201

Difference Scheme for the Vlasov-Manev System 195 Figure 3: The solution and its iso-lines for γ = 4, δ = 0

Documenta Mathematica 4 (1999) 179–201

We observe a very clear unification of the two particle “clouds” in space and no remarkable concentration of mass during the time evolution. Figure 3 shows the function f(t, x, v) and its iso-lines for the time t = 1.4. In the second test we consider the pure Manev case with the same initial distribution and the same parameter of discretisation.

Figure 4: The density profiles for γ = 0, δ = 4

Figure 5: The force profiles for γ = 0, δ = 4

Documenta Mathematica 4 (1999) 179–201

Difference Scheme for the Vlasov-Manev System 197 Figure 6: The solution and its iso-lines for γ = 0, δ = 4

Documenta Mathematica 4 (1999) 179–201

There is a very clear difference in the behaviour of the two examples. The pure Manev case leads to a significant concentration of the mass in the two “clouds”, and during the evolution they remain separated. Figure 6 shows the function f(t, x, v) and its iso-lines for the time t = 1.4. Finally, we consider the mixed case γ = 2, δ = 2 in order to illustrate the influence of the two effects: unification and concentration. The results are presented in Figures 7,8,9.

Figure 7: The density profiles for γ = 2, δ = 2

Figure 8: The force profiles for γ = 2, δ = 2 Difference Scheme for the Vlasov-Manev System 199 Figure 9: The solution and its iso-lines for γ = 2, δ = 2

Documenta Mathematica 4 (1999) 179–201

5 Conclusions

Our calculations suggest that the Manev correction will have a “stabilising” effect on isolated one-dimensional matter concentrations; this effect counteracts the tendency of the long-range Newtonian potential to accumulate all matter in one location; while this is only an isolated phenomenon which is observed here as a consequence of the Manev correction, we believe it to be evidence that truly interesting effects may occur in the more relevant three-dimensional case. Numerical experiments to this end are planned.

References

[1] A. Bobyblev and N. Ibragimov. The relationship of the group symmetry properties in dynamics, kinetic theory and hydrodynamics series. Mathem. Modelirovanie, 1 : 100–109, 1989.

[2] A. V. Bobylev, P. Dukes, R. Illner, and H. D. Victory. On Vlasov-Manev equations, I: Foundations, Properties and Global Nonexistence. JSP, 88:885–911, 1997.

[3] J. W. Cooley and J. W. Tukey. An algorithm for the machine calculation of complex Fourier series. Math. Comput., 19 : 297–301, 1965.

[4] R. Glassey. The Cauchy Problem in Kinetic Theory. SIAM Publ., 1996.

[5] L. Greengard and V. Rokhlin. A Fast Algorithm for Particle Simulations. J. Comput. Phys., 73: 325–348, 1987.

[6] R. Illner, P. Dukes, H. D. Victory, and A. V. Bobylev. On Vlasov-Manev equations, II: Local Existence and Uniqueness. JSP, to appear, 1998.

[7] G. Manev. La gravitation et le principe de l’égalité de l’action et de la réaction. Comptes Rendues, 178: 2159–2161, 1924.

[8] G. Manev. Die Gravitation und das Prinzip von Wirkung und Gegenwirkung. Zeitschrift für Physik, 31: 786–802, 1925.

[9] G. Manev. La gravitation et l’eńrgie au zéro. Comptes Rendues, 190: 1374– 1377, 1930.

[10] G. Manev. Le principe de la moindre action et la gravitation. Comptes Rendues, 190: 963–965, 1930.

[11] H. Minc. Nonnegative Matrices. J. Wiley & Sons, 1988.

[12] S. Rjasanow. Effective algorithms with circulant-block matrices. Linear Alg. Appl., 202: 55–69, 1994.

[13] C. Van Loan. Computational Frameworks for the Fast Fourier Transform. SIAM, Philadelphia, 1992.

R. Illner Department of Mathematics and Statistics University of Victoria P.O. Box 3045 Victoria, B.C. V8W 3P4 Canada rillner@Math.UVic.CA S. Rjasanow Fachbereich 9 - Mathematik Universität des Saarlandes Postfach 151150 66041 Saarbrücken Germany rjasanow@num.uni-sb.de 202

Documenta Mathematica 4 (1999)