Kinetic layers and coupling conditions for macroscopic equations on networks I: the wave equation
Raul Borsche Axel Klar





[1]Technische Universität Kaiserslautern, Department of Mathematics, Erwin-Schrödinger-Straße, 67663 Kaiserslautern, Germany 
  ({borsche, klar}@mathematik.uni-kl.de)
[2]Fraunhofer ITWM, Fraunhoferplatz 1, 67663 Kaiserslautern, Germany 
 




Kinetic layers and coupling conditions for  macroscopic equations on networks I: the wave equation
    R. Borsche[1] 
     A. Klar[1] [2]
    
==================================================================================================




We consider kinetic and associated macroscopic equations on networks. 
The general approach will be explained in this paper for a linear kinetic BGK model and the corresponding limit for small Knudsen number, which is  the wave equation.
Coupling conditions for the macroscopic equations are derived from the kinetic conditions via an asymptotic analysis near the nodes of the network. 
This analysis leads to the consideration of a fixpoint problem involving the coupled solutions of kinetic half-space problems. 
A new approximate method for the solution of kinetic half-space problems is derived and used for the determination of the coupling conditions. 
Numerical comparisons between the solutions of the macroscopic equation with different coupling conditions and the kinetic solution are presented for the case of tripod and more complicated networks.


Keywords. 
Kinetic layer, coupling condition, kinetic half-space problem, networks

AMS Classification.  
 82B40, 90B10,65M08






§ INTRODUCTION


There have been many attempts  to define coupling conditions for macroscopic partial differential equations on networks including, for example,  drift-diffusion equations, scalar  hyperbolic equations, or hyperbolic systems like  the wave equation or Euler type models, see for example <cit.>. 
In <cit.> coupling conditions for scalar hyperbolic equations on networks are discussed and investigated.
<cit.> treat the wave equation and general nonlinear hyperbolic systems are considered in <cit.>. 
We finally note, that, for example, for hyperbolic systems  on networks 
there are still many unsolved problems, like finding suitable coupling conditions without restricting to subsonic situations.

On the other hand, coupling conditions for kinetic equations on networks have been discussed in a much smaller number of publications,
see <cit.>. 
In <cit.> a first attempt to derive a coupling condition for a
macroscopic equation from the underlying kinetic model has been presented for the case of a kinetic equations for chemotaxis.
In the present  paper, we will present a more general and more accurate procedure.
It is motivated by the classical procedure to find kinetic slip boundary conditions for macroscopic equations based on the analysis of the 
kinetic layer <cit.>.
In this work, we will derive coupling conditions for macroscopic equations 
on a network from underlying microscopic or kinetic models via an asymptotic analysis of the situation near the nodes. 
To explain the basic approach we concentrate on a simple linear BGK-type kinetic model with the linear wave equation as the associated macroscopic model. 
More complicated problems and, in particular, nonlinear models will be discussed in future work. 

The paper is organized in the following way. 
In section <ref> we discuss the kinetic and macroscopic equations and the boundary and coupling conditions for these equations.
In section <ref>  kinetic boundary layers are discussed, as well as   an 
asymptotic analysis of the kinetic equations near the nodes. 
This leads to an abstract formulation of the coupling conditions for the macroscopic equations at the nodes based on a fix-point problem involving kinetic half-space equations. 
In the following section <ref> an approximate coupling condition is derived based on the so called Maxwell approximation of the half-space problem. 
A  refined method to determine the solution of the half space problems is derived, compared to previous approximate solution methods for half-space problems and applied to the problem of finding accurate coupling conditions for the macroscopic equations in sections <ref>. 
Moreover, the macroscopic equations on the network with the different coupling conditions are numerically compared to each other and to the full solutions of the kinetic equations on the network  in section <ref>.
The results show the very good approximation of the underlying kinetic model by the macroscopic model with the new coupling conditions.




§ EQUATIONS AND BOUNDARY AND COUPLING CONDITIONS



We consider  the  linear kinetic  BGK model in 1D for  
  x ∈ℝ, v  ∈.

    ∂_t f + v ∂_x f = -1/ϵ(f-(ρ + v/ a^2  q ) M (v))

with the Maxwellian

    M (v) = 
    1/√(2 π a^2)exp( - | v|^2/2 a^2),

where

    ρ_ϵ = ∫_-∞^∞ f(v) dv ,      q_ϵ = ∫_-∞^∞ v f(v) dv .

The associated macroscopic equation  for ϵ→ 0 is the wave equation

    ∂_t ρ + ∂_x q =0
    ∂_t q + ∂_x (a^2 ρ ) =0 .


To illustrate the influence of the underlying kinetic model on the coupling conditions for the macroscopic equations, we also consider the following equation with  bounded velocity space  v  ∈ [-1,1] 

    ∂_t f + v ∂_x f = -1/ϵ(  f-  (ρ_ϵ+ v/a^2  q_ϵ ) 1/2)

with a^2 = 1/3 or 

    ∂_t f + v ∂_x f = -1/ϵ(  f-  (ρ_ϵ/2+ 3/2 v q_ϵ ) )

with

    ρ_ϵ = ∫_-1^1 f(v) dv ,    q_ϵ = ∫_-1^1 v f(v) dv  .

The associated macroscopic equation is again (<ref>) 
with
a^2 = 1/3.


 §.§ Boundary conditions

For x ∈ [0,b] we prescribe for the kinetic  equation  

    f(0,v), v > 0,   f(b,v), v<0.

For (<ref>) the boundary conditions are given in characteristic variables <cit.>.
The corresponding Riemann Invariants are

    r_1/2 = q   ∓ a ρ.

As boundary data the value of q + aρ
at the left boundary and q - aρ at the right boundary are prescribed.



 §.§ Coupling conditions


If these equations are considered on a network, it is sufficient to study a single coupling point.
At each node so called coupling conditions are required. 
In the following we consider a node connecting n edges, which are oriented away from the node, as in Figure <ref>.
Each edge i is parametrized by the interval [0,b_i] and the kinetic and macroscopic quantities are denoted by
f^i and ρ^i, q^i respectively.
A possible choice of coupling conditions for the kinetic problem are given by 

    f^i(0,v) = ∑_j=1^n c_ijf^j (0,-v), v >0,

compare <cit.>.
The total mass in the system is conserved, if 

    ∑_i=1^n c_ij = 1

 holds.
In the following we use the vector notation 

    f^+ = C f^-, v >0 ,

where f^+=(f^1(0,v),…,f^n(0,v)) and f^-=(f^1(0,-v),…,f^n(0,-v)).
The coupling conditions for the macroscopic quantities are
conditions on the characteristic variables

    r_2^i (0) = q^i(0) + a ρ^i(0)

using the given values of 

    r_1^i (0) = q^i(0) - a ρ^i(0) .

We refer to <cit.> for systems of macroscopic equations on networks.
In the following we will derive, via asymptotic analysis, 
macroscopic coupling conditions from the kinetic coupling conditions (<ref>).




§ BOUNDARY AND COUPLING CONDITIONS FOR MACROSCOPIC EQUATIONS VIA KINETIC LAYER ANALYSIS




 §.§ Boundary conditions

A kinetic layer analysis, see <cit.>, at the left boundary of the interval [0,b], i.e.  a rescaling of the spatial variable in equation (<ref>) with ϵ, gives   to first order in ϵ  the  following stationary kinetic  half space  problem for 
x ∈ [0,∞]

    v ∂_x φ = - ( φ -(ρ + v/ a^2 q ) M (v) ) ,

where ρ and q are the zeroth and first moments of  φ. 
At x=0 the boundary conditions for the  half space problem are

    φ(0,v) =  k(v) = f(0,v), v > 0 .


On the left side, i.e. at x=∞,  a condition is prescribing 
an arbitrary linear combination of the invariants of the half-space problem
< v φ >
and < v^2φ > .
Here, and in the following we use the notation < φ  > = ∫φ dv and  < φ  >_+ = ∫_v>0φ dv
or < φ  >_- = ∫_v<0φ dv.
We use the values of the first Riemann Invariant (<ref>) r_1 = q -a  ρ of the macroscopic system (<ref>) to fix 

    < (v - v^2/ a )φ> = r_1.



The boundary condition for (<ref>) is obtained by 
determining r_2 from  the asymptotic solution
of the half space problem and setting r_2 = q_∞ +a  ρ_∞.
The values ρ_∞ and q_∞ are the macroscopic quantities associated to the solution of the half-space problem at infinity, which has the form

    φ(∞,v) = (ρ_∞+ v/a^2 q_∞)M(v) .


The solution of the half space problem is also used to determine the outgoing distribution

    A[k](v) = f(0,v)  = φ(0,v), v<0 ,

where the notation A is used for the so called Albedo operator of the  half space problem. 
The structure of a half space problem is illustrated in Figure <ref>.







 §.§ Coupling conditions


We use the corresponding procedure to determine the coupling conditions for the macroscopic equations.


Starting from the kinetic coupling conditions

    f^i(0,v) = ∑_j=1^n c_ijf^j (0,-v), v >0

we determine
the coupling conditions for the macroscopic equations in the following way.
We use the kinetic  coupling conditions to obtain conditions on the in- and outgoing solutions of the half space problems on the different arcs.
That means 

    φ^i(0,v)
    = ∑_j=1^n c_ijφ^j(0,-v) , v>0

or, if we denote the ingoing function of the half-space problem on arc i by
k^i (v),v>0  and the outgoing solution by A^i [k^i] (v), v<0

    k^i(v)
     = ∑_j=1^n c_ij  A^j [ k^j] (-v) , v  >0 .



This is a fix point equation for k^i, i= 1 …,n.

Additional conditions are needed to solve the half-space problems, i.e.


    < ( v- v^2/a )φ^i>= r^i_1

with r^i_1 =  q^i -    a ρ^i.
The coupling conditions for the wave equation  are  conditions on the outgoing characteristic variables at x=0. We define 

    r^i_2 (0) = 
     q^i_∞[k^i]+ a ρ^i_∞[k^i].


The main task is now to find tractable expressions for these 
coupling conditions. In the next section we discuss the 
so called Maxwell method to solve the half-space problem approximately and the resulting approximate coupling conditions.
In further sections a new refined method to determine 
the solution of the half space problems based on half moment equations is derived and applied to the problem of finding accurate coupling conditions. 
Finally, the results are compared to full solutions of the kinetic equations on the network and to other approximate solution methods for the half-space problem.




§ APPROXIMATE SOLUTION  OF THE HALF SPACE PROBLEM VIA HALF-FLUXES AND APPROXIMATE COUPLING CONDITIONS



To solve kinetic  half space problem approximately a variety of different methods can be found in the literature. 
Approaches via a Galerkin method can be found in <cit.>. 
Approximate  methods to determine only the  asymptotic states and outgoing distributions can be found in <cit.>.
For the determination of the macroscopic  boundary conditions we use in this section  a simple approximation, the so called Maxwell approximation, see <cit.>. 



 §.§ Approximate solution of the half space problem



We use  the equality of half-fluxes

    <v k(v) >_+  
    = <v (ρ_∞ + v/ a^2 q_∞) M >_+

to obtain one condition. 
Together with the first Riemann Invariant

    q_∞ -aρ_∞
     = q -a ρ = r_1(0) = C =

this determines ρ_∞ and q_∞.
This can be rewritten as

    ([ a 1/√(2 π)        1/2;         -a          1 ])
    ([ ρ_∞; q_∞ ])
    =
    ([ <vk>_+;      C ])

with the solution

    ([ ρ_∞; q_∞ ])
    =
    1/a( 1/√(2 π)- 1/2)([          <vk>_+- C/2; a <vk>_+ + aC/√(2 π) ])
     .

The outgoing distribution is  approximated by

    φ(0,v)  =  ( ρ_∞ + v/ a^2  q_∞) M (v ) , v <0.



For bounded velocity space and equation (<ref>) we have

    ([ 1/4 1/2;  -a   1 ])
    ([ ρ_∞; q_∞ ])
    =
    ([ <vk>_+;      C ])
     .





 §.§ Approximate coupling conditions


The fix-point problem (<ref>) is in the present case 
approximated by the problem

    k^i(v) = ∑_j=1^n c_ij( ρ_∞^j [k^j]- v/ a^2  q_∞^j [k^j]
    ) M (-v) ,  v >0.

Thus, the asymptotic values are determined by

    ([ a 1/√(2 π)        1/2;         -a          1 ])
     ([ ρ^i_∞; q^i_∞ ])
    =
    ([ ∑_j c_ij( ρ_∞^j a/√(2 π) -  q_∞^j/ 2);                             r^i_1 (0) ]).

These values can be assigned to the states of the wave equation (<ref>) at the coupling point by

    r^i_2 (0) = q^i_∞ +a ρ^i_∞ .

Summing the first equation of (<ref>) for all i=1,…,n and using the conservation property of the kinetic coupling conditions (<ref>) directly yields the conservation of mass in the macroscopic variables

    ∑_i=1^n q_∞^i =0.

In the special case of a  uniform node, i.e. c_ij = 1/n-1 for i ≠ j and c_ij=0 for i =j  we have

    a 1/√(2 π)ρ^i_∞ +   1/2  q^i_∞
    
    
    =1/n-1∑_j i(a/√(2 π)ρ^j_∞  -   q^j_∞/ 2).

Multiplying by √(2 π) and adding √(2 π)/n times (<ref>) we obtain

    a  ρ^i_∞ +   √(2 π)/2n-2/n q^i_∞ 
    =1/n-1∑_j i(a  ρ^j_∞  +√(2 π)/2n-2/n   q^j_∞)     i = 1,…,n.

These n equations are not linearly independent and they can be reformulated as

    a ρ^i_∞ +   √(2 π)/2n-2/n q^i_∞
    
    
    =
    a   ρ^j_∞  +   √(2 π)/2n-2/n q^j_∞     i,j=1,…,n , i≠ j.

Thus we have found a macroscopic invariant at the coupling point, i.e.

    ρ +   √(2 π)/2 a n-2/n q.




For bounded velocity space we have

    1/4ρ^i_∞ +   1/2  q^i_∞
    
    
    =∑_j c_ij(
     1/ 4ρ^j_∞ -   q^j_∞/ 2).

For a uniform node this leads to

    ρ^i_∞ +   2(n-2)/n  q^i_∞
    
    
    =
     ρ^j_∞ +   2(n-2)/n q^j_∞

and to the invariant

    ρ_∞ +   2(n-2)/n  q_∞ .

For example for n=3, we obtain for  the above invariant (<ref>)  a factor 
√(6 π)/6∼ 0.72360 for a^2 = 1/3 in contrast to the value 2/3 for the bounded velocity case. 
Thus although the kinetic coupling conditions and the macroscopic equations (<ref>) are identical we obtain different coupling conditions for a model with bounded velocity, since the half space problems are different.



 Another direct approach to obtain coupling conditions for the wave equation has been used in <cit.>.
 Using a full moment approximation of the distribution function in the case of bounded velocities, i.e. using
 
    f^± (v) = ρ/2±3/2 v q,    v∈ [0,1]

 in the kinetic coupling conditions  (<ref>), one obtains after integration over the positive velocities
 
    ρ^i_∞  + 3/2 q^i_∞ =    ∑_j c_ij(ρ^j_∞ - 3/2 q^j_∞)   .
 
 This leads, for a uniform node, to the invariant
   
    ρ_∞   + 3(n-2)/2n q_∞ .
 
      In the case of unbounded velocities we obtain 
      
    ρ^i_∞  + √(2)/ a √(π) q^i_∞ =    ∑_j c_ij(ρ^j_∞ - √(2)/ a √(π) q^j_∞)   .
 
      and the invariant
         
    ρ_∞   + √(2)/ a √(π)(n-2)/n q_∞ .

    
  
 
 



§ COUPLING CONDITIONS VIA  APPROXIMATION BY HALF-MOMENT EQUATIONS


In this section we develop a refined method to solve the half-space problem based on a half moment approximation and use it for the derivation of refined coupling conditions.
We begin with the case of equation (<ref>).
For further approaches to the approximate solution of half space problems we refer to <cit.>.



 §.§ The case of a bounded velocity domain


Consider  the linear BGK model with bounded velocities (<ref>) and the corresponding limit equation (<ref>).



  §.§.§ Half moment equations
 


We determine a half moment approximation for the half-space solution, compare <cit.>.
We define

    ρ_ϵ^-    =∫_-1^0 f(v)dv ,        ρ_ϵ^+    =∫_0^1f(v)dv ,
    
    
     q_ϵ^-    =∫_-1^0 vf(v)dv  ,
        q_ϵ^+    =∫_0^1vf(v)dv  .

As  closure assumption we use the following approximation of the distribution function f by  affine linear functions in v to determine half-moment equations, see <cit.> and references therein.

    f(v) = a^+ + vb^+, v≥ 0      and      
    f(v) = a^- + vb^-, v≤ 0
     .

One obtains

    ρ_ϵ^-    = a^- -1/2b^-
     ,     ρ_ϵ^+    =a^+ + 1/2b^+ , 
    
    q_ϵ^-    = -1/2a^- + 1/3b^- 
     ,     
    q_ϵ^+    = 1/2a^+ + 1/3b^+

and 

    ∫_0^1 v^2f(v)dv = -1/6ρ_ϵ^+ + q_ϵ^+  ,
          ∫_-1^0 v^2f(v)dv = -1/6ρ_ϵ^- - q_ϵ^- .

Finally, integrating the kinetic equation with respect to the corresponding half-spaces, we get the half-moment approximation of the  kinetic equation  as

    {[                 ∂_tρ_ϵ^+ + ∂_xq_ϵ^+                                   =    -1/ϵ(ρ_ϵ^+ - (ρ_ϵ/2 + 3 q_ϵ/4) );   ∂_tq_ϵ^+ + ∂_x(-1/6ρ_ϵ^+ + q_ϵ^+)                                   =       -1/ϵ(q_ϵ^+ -(ρ_ϵ/4 +  q_ϵ/2));                 ∂_tρ_ϵ^- + ∂_xq_ϵ^-                                   =    -1/ϵ(ρ_ϵ^- -
    (ρ_ϵ/2 - 3 q_ϵ/4)
    ); ∂_tq_ϵ^- + ∂_x( -1/6ρ_ϵ^- - q_ϵ^- )                                   =  -1/ϵ(q_ϵ^- - (- ρ_ϵ/4 +  q_ϵ/2)) .;                                     ].

Introducing the even-odd variables

    ρ_ϵ   = ρ_ϵ^+ + ρ_ϵ^- ,
           ρ̂_ϵ   = ρ_ϵ^+ - ρ_ϵ^- ,
    
    
    q_ϵ   = q_ϵ^+ + q_ϵ^- , 
       q̂_ϵ   = q_ϵ^+ -q_ϵ^- ,

we can rewrite the system as 

    {[               ∂_tρ_ϵ + ∂_xq_ϵ                             =                             0;  ∂_tq_ϵ + ∂_x(-1/6ρ_ϵ + q̂_ϵ)                             =                             0;             ∂_tρ̂_ϵ + ∂_xq̂_ϵ                             =            -1/ϵ(ρ̂_ϵ -3/2q_ϵ); ∂_tq̂_ϵ + ∂_x(-1/6ρ̂_ϵ + q_ϵ)                             =            -1/ϵ(q̂_ϵ-ρ_ϵ/2) . ].

Obviously, the  half-moment model  has again the wave equation (<ref>) as  macroscopic limit as ϵ goes to 0. 




  §.§.§ Half space-half moment problem


Rescaling the spatial variable in the half-moment problem with ϵ one obtains the following half-space problem for 
x ∈^+


    {[                  ∂_xq^+                       = -(ρ^+ - (ρ/2 + 3 q/4) );      ∂_x(-1/6ρ^+ + q^+)                       =    -(q^+ -(ρ/4 +  q/2));                  ∂_xq^-                       = -(ρ^- -
    (ρ/2 - 3 q/4)
    );    ∂_x( -1/6ρ^- - q^- )                       = -(q^- - (- ρ/4 +  q/2));                         ].

or

    {[            ∂_xq               =               0; ∂_x(-1/6ρ + q̂)               =               0;           ∂_xq̂               =     -(ρ̂ -3/2q); ∂_x(-1/6ρ̂ + q)               =     -(q̂-ρ/2) . ].

We  have to provide boundary conditions for ρ^+(0) and q^+ (0), as well as a condition at x=∞

    q_∞ - a ρ_∞ = r_1 (0) = C.

Then, the half space problem can be solved explicitly.
We determine a solution up to 3 constants which will be fixed with the above 3 conditions.
First, we observe, that we have 2 invariants

    q    = C_1 = const 
    
    -ρ/6 + q̂   = C_2

From the last equation in (<ref>) we can deduce that at x = ∞

    q̂_∞ = ρ_∞/2 .

Combining this with (<ref>) gives 
ρ_∞ = 3 C_2 or q̂_∞ = 3 C_2/2.
From the third equation of (<ref>) we obtain 
ρ̂_∞ = 3q/2 = 3 C_1/2.
This simplifies (<ref>) to 

    {[                     q                     =                   C_1;                     ρ                     =          6 q̂ - 6 C_2;                 ∂_xq̂                     =        -(ρ̂ -3/2 C_1);           ∂_x(-1/6ρ̂)                     = -(-2 q̂ +  3 C_2  ) . ].

The ODEs for ρ̂ and  q̂ have the solutions 

    ρ̂= γexp(-2x/a)   + γ̂exp(2x/a) + 3 C_1/2
    q̂ = a/2γexp(-2x/a)   -  a/2γ̂exp(2x/a) +3 C_2/2

Since we are looking only for bounded solutions we are left with

    ρ̂= γexp(-2x/a)    + 3 C_1/2
    q̂ = a/2γexp(-2x/a)   +3 C_2/2
    
    q=C_1
    ρ = 3 a γexp(-2x/a)   + 3  C_2 .

The three  parameters are  fixed with the 3 conditions mentioned above.
At x=0 inflow data is given

    1/2( q(0)+q̂ (0) ) = q_+ (0)
    1/2( ρ (0)+ρ̂(0) ) = ρ_+ (0)

and the Riemann Invariant at x=∞ gives

    q_∞ -a ρ_∞ = C.

Inserting the above determined solution we obtain

    1/2( C_1 +a/2γ   +3 C_2/2)    = q_+ (0)
    1/2( 3 a γ   + 3  C_2 +γ     + 3 C_1/2)    = ρ_+ (0) ,

 which can be rewritten in terms of the asymptotic values

    q_∞/2   +ρ_∞/4 +a/4γ   = q_+ (0)
    3 q_∞/4 +  ρ_∞/2       + (3 a+1/2) γ   = ρ_+ (0)  .
 
  Together with the condition at infinity (<ref>), this determines  the   asymptotic values q_∞,ρ_∞ and γ.
 The   outgoing quantities ρ_-(0), q_-(0) are then determined by 
 
    q_∞/2   - ρ_∞/4 - a/4γ   = q_-(0)
    
      - 3 q_∞/4  +    1 /2ρ_∞     + (3 a-1 /2) γ   = ρ_- (0)  .
 
 
 
With the Maxwell approximation for the half-moment problem in (<ref>) the asymptotic states are determined by 
 
    q_+ (0) = q_+(∞)  = 1/2 (q_∞ + q̂_∞)
     = q_∞/2 + ρ_∞/4
 
 and the condition at infinity (<ref>).
 The outgoing quantities are 
 
    q_- (0)    = q_-(∞) = - ρ_∞/4 +  q_∞/2
    ρ_-(0)    = ρ_-(∞)
      = ρ_∞/2 -  3 q_∞/4 .
 



  §.§.§ The extrapolation length

 
 To estimate the accuracy of our method, we consider the classical problem of determining the so-called extrapolation length <cit.>.
 
 For x ∈ℝ^+, v  ∈ [-1,1] we consider the half space 
   equation 
  
    v ∂_x f = - (  f-  (ρ/2+ 3/2 v q ))

  with
  ∫_-1^1 v f dv = q = 0 and f (0,v) = v , v >0.
  That means, we consider
  
    v ∂_x f = - (  f-  ρ/2) .

  The extrapolation length is the value of λ_∞ = f(∞,v) = ρ_∞/2.
  The Maxwell approximation  gives λ_∞ =   2/3.
  The above half moment approximation gives
  
    ρ_∞/4 +a/4γ  = q_+ (0) = 1/3
    ρ_∞/2       + (3 a+1/2) γ  = ρ_+ (0) = 1/2 .
 
  This leads to 
   
    ρ_∞    =   9 a+ 4/6 a+3
 
   and with a^2 =1/3 we obtain
   
    ρ_∞    =   3 √(3) + 4/2 √(3)+3  .
 
   Thus, the extrapolation length is approximated as λ_∞∼ 0.7113.
   The exact value computed from a spectral method is 0.7104,
   see <cit.>. 
   This yields an error for the above half-moment method of approximately 0.1 %.
   In contrast, the Maxwell approximation gives 0.6666, which is an error of 6.1 %. 
   The variational method in <cit.> gives 0.7083, which is an error of 0.3%.
   
  


  §.§.§ Half-moment coupling conditions


In this subsection we determine the coupling conditions on the basis of the half-moment approximation of the half-space problem.
Multiplying with v and integrating the kinetic coupling conditions (<ref>) with respect to the  positive and negative half moments  gives

    q^i_+ (0) = - ∑_ j=1^n  c_ij q_-^j (0)  .

Inserting the half moment approximations (<ref>) yields

    q^i_∞/2   +ρ^i_∞/4 +a/4γ^i  
    =  ∑_ j=1^n  c_ij( -q^j_∞/2   + ρ^j_∞/4 + a/4γ^j).
 
Again, a summation w.r.t. i=1,…,n directly gives 
the equality of fluxes

    ∑_i q^i_∞ =0 .

For a uniform node  with equal distribution c_ij = 1/n-1, i ≠ j and 0 otherwise, 
one obtains using the equality of fluxes

    n-2/2nq^i_∞   +ρ^i_∞/4 +a/4γ^i  
    = n-2/2nq^j_∞   + ρ^j_∞/4 + a/4γ^j .
 
or the invariance of 

    ρ_∞ +  2(n-2)/nq_∞ + aγ.



Further, integrating the kinetic coupling conditions (<ref>)
with respect to the positive and negative  half moments, we obtain 

    ρ^i_+ (0) =  ∑_j =1^n  c_ijρ_-^j (0) .
 
With the half moment approximations (<ref>) this reads

    3 q^i_∞/4 +  ρ^i_∞/2       + (3 a+1/2) γ^i  
      =  ∑_j=1^n  c_ij( - 3 q^j_∞/4  +    1 /2ρ^j_∞     + (3 a-1 /2) γ^j) .
 
 Summing these conditions for i=1,…,n yields 
 
    ∑_i=1^nγ^i=0 .

 Thus, in the case of a uniform node, we  derive another  coupling invariant



  
    ρ_∞    
          + 3(n-2)/2n q_∞+ (3 a+ n-2/n) γ .
 
  
  Alltogether this yields 2n  conditions at a node, i.e. (<ref>), (<ref>), (<ref>) and (<ref>). 
  In combination with the conditions at infinity we have 3n conditions for 3n quantities γ^i, ρ^i_∞, q^i_∞.
   
Note that the invariants (<ref>) and  (<ref>) can be combined such that γ is eliminated, which gives the invariance of
 
    ρ_∞ + (n-2/n)9a+4(n-2/n)/4a+2(n-2/n) q_∞ .
 


 All the  coupling conditions for the wave equation with uniform nodes derived up to now are given by the conservation of mass (<ref>) and an invariant of the form 
 
    ρ   + C q .
 
 They differ in the factor C, see Table <ref>. 
 

	
 			
 
 
   With a =1 /√(3) and n=3 the half moment approximation  gives  C=0.7313
   compared to 0.6666 for Maxwell (<ref>) or 0.5 for the full moment approximation.
   
   A numerical comparison of the resulting network solutions is presented in the next section.
   
    

	For the solution of the wave equation (<ref>), we consider the mathematical entropy
	
    e = 1/2(ρ^2+1/a^2q^2).

	It evolves according to the conservation law
	
    ∂_t e + ∂_x (ρ q) = 0  .

	Along one edge this entropy  is conserved, but the total entropy in the network can change according to the entropy-fluxes at the nodes.
	Note that for all above models with ρ + Cq= C̃ and C>0 the total entropy decays, since
	
    ∑_i=1^nρ  q = 
    		∑_i=1^n(C̃-Cq)  q = 
    		C̃∑_i=1^nq -C∑_i=1^nq^2 =  -C∑_i=1^nq^2 <0 .
	
   
 In the following we discuss also the case of an underlying kinetic model with unbounded velocities.


 §.§ The case of unbounded velocity domain

We now apply the above procedure to the linear BGK model (<ref>) with an unbounded velocity space, which has (<ref>) with arbitrary a as limit equation.



  §.§.§ Half moment equations

The half moments are defined as

    ρ ^-    =∫_-∞^0 f(v)dv ,        ρ^+    =∫_0^∞ f(v)dv ,
    
    
    q ^-    =∫_-∞^0 vf(v)dv  ,
        q^+    =∫_0^∞ vf(v)dv  .

As closure we consider the following approximations of the distribution function  

    f(v) = (a^+ + vb^+)M(v), v≥ 0      and      
    f(v) = (a^- + vb^-)M(v), v≤ 0
     ,

which leads to

    ρ^-    = 1/2a^- -a/√(2π)b^-
     ,     ρ^+    =1/2a^+ + a/√(2π)b^+ , 
    
    q^-    = -a/√(2π)a^- + a^2/2b^- 
     ,     
    q^+    = a/√(2π)a^+ + a^2/2b^+ .

Inverting these relations


























gives

    ∫_0^∞ v^2f(v)dv    = 
    
     1/π - 2((π - 4)a^2 ρ^++ a√(2π)q^+)
     ,
    
    ∫_-∞^0 v^2f(v)dv    = 
    
     1/π - 2((π - 4)a^2 ρ^– a√(2π)q^-)
     .

This leads to the half-moment system 

    {[                                  ∂_tρ^+ + ∂_xq^+                                                =                   -1/ϵ(ρ^+ - (ρ/2 +  q/√(2π)a) ); ∂_tq^+ + ∂_x(1/π - 2((π - 4)a^2 ρ^++ a√(2π)q^+))                                                =                     -1/ϵ(q^+ -(aρ/√(2π) +  q/2));                                  ∂_tρ^- + ∂_xq^-                                                =                  -1/ϵ(ρ^- -
    (ρ/2 -   q/√(2π)a)
    ); ∂_tq^- + ∂_x( 1/π - 2((π - 4) a^2ρ^– a√(2π)q^-))                                                =                -1/ϵ(q^- - (- aρ/√(2π) +  q/2)) .;                                                  ].

Introducing the even-odd variables as before 












we can rewrite the system as 

    {[                             ∂_tρ + ∂_xq                                       =                                       0; ∂_tq + ∂_x(π-4/π-2a^2ρ +  a√(2π)/π-2q̂)                                       =                                       0;                           ∂_tρ̂ + ∂_xq̂                                       =                    -1/ϵ(ρ̂ -2/√(2π) aq); ∂_tq̂ + ∂_x(π-4/π-2a^2ρ̂ +a√(2π)/π-2 q)                                       =                    -1/ϵ(q̂-2aρ/√(2π)) . ].









  §.§.§ Half space-half moment problem


The corresponding half space problem is for 
x ∈^+ given by

    {[                             ∂_xq                                =                                0; ∂_x(π-4/π-2a^2ρ +  a√(2π)/π-2q̂)                                =                                0;                            ∂_xq̂                                =                -(ρ̂ -2/√(2π) aq);  ∂_x(π-4/π-2a^2ρ̂ +a√(2π)/π-2 q)                                =                -(q̂-2aρ/√(2π)) . ].

As before we need boundary conditions on 
ρ^+(0), q^+ (0)
and a condition at infinity

    q_∞ - a ρ_∞ = r_1 = C.

We now construct the explicit solution of this half space problem. 

The first two equations of (<ref>) state the invariance of 

    q    = C_1 = const 
    π-4/π-2a^2ρ +  a√(2π)/π-2q̂   = C_2

or 





    ρ = π-2/a^2(π-4)C_2 - √(2π)/a(π-4)q̂.

From the last equation of (<ref>) we can deduce that at x = ∞ we have

    q̂_∞ = 2a/√(2π)ρ_∞ ,

which leads with the above invariant to ρ_∞ = 1/a^2C_2.
Moreover, from the third line in (<ref>) we obtain

    ρ̂_∞ = 2/√(2π) aq_∞





and thus we can transform the lower part of (<ref>) to

    {[                                                  ∂_xq̂                                                      =                                    -(ρ̂ -2/√(2π) aq_∞);                      ∂_x(π-4/π-2a^2ρ̂ +a√(2π)/π-2 q_∞)                                                      = -( q̂ - 2aρ/√(2π)(π-2/a^2(π-4)C_2 - √(2π)/a(π-4)q̂)) . ].

Rearranging gives

    {[                       ∂_xρ̂                           =     -(π-2)^2/(π-4)^2a^2q̂ + 2(π-2)^2/√(2π)a^3(π-4)^2C_2;                       ∂_xq̂                           =                       -ρ̂ +              2/√(2π) aq_∞ . ].

By defining λ =  π-2/a(π-4) the solution of ρ̂ and q̂ is

    ρ̂= γλexp(λ x)   + γ̂λexp(-λ x) +  2/√(2π) aq_∞
    q̂ = -γexp(λ x)   + γ̂exp( -λ x) +2a/√(2π)ρ_∞ .

Since we are only interested in bounded solutions, we have γ̂=0 and thus

    ρ̂= γλexp(λ x) +  2/√(2π) aq_∞
    q̂ = - γexp(λ x)  +2a/√(2π)ρ_∞ .

As before, there are three parameters which can be determined with the three conditions

    1/2( ρ (0)+ρ̂(0) ) = ρ_+ (0) 
    1/2( q(0)+q̂ (0) ) = q_+ (0)
    
    q_∞ -a ρ_∞ = r_1 (0) .

Inserting the expressions (<ref>) yields 

    1/2(ρ_∞+(λ +  √(2π)/a(π-4))γ  +  2/√(2π) aq_∞)    = ρ_+ (0) 
    1/2( q_∞ -γ  +2a/√(2π)ρ_∞)    = q_+ (0)

and for the outgoing quantities ρ_-(0) and q_-(0) we obtain

    1/2(  ρ_∞+(-λ + √(2π)/a(π-4))γ  -  2/√(2π) aq_∞)    = ρ_- (0) 
    1/2( q_∞ + γ  -2a/√(2π)ρ_∞)    = q_- (0) .



The Maxwell approximation  is in this case given by the equation

    q_+ (0) = q_+(∞)  = 1/2 (q_∞ + q̂_∞)
    = q_∞/2 + a/√(2π)ρ_∞ .
 
The outgoing quantities are computed as

    q_- (0) = q_-(∞) =   q_∞/2 - a/√(2π)ρ_∞ ,
    ρ_-(0) 
    = ρ_∞/2 -  1/√(2π) aq_∞ .
 

    


  §.§.§ The extrapolation length
  
In order to estimate the quality of the above approximation, we consider

    v ∂_x f = - (f-(ρ + v/ a^2  q ) M (v))

with q =0 and k(v) = v M, in order to determine f(∞,v ) = λ_∞ M.
The Maxwell approximation gives
ρ_∞ = √(2 π)/2 a.
With the half moment method we obtain

    1/2(ρ_∞+(λ +  √(2π)/a(π-4))γ  +  2/√(2π) aq_∞)    = ρ_+ (0)  = 1/√(2 π)
    1/2( q_∞ -γ  +2a/√(2π)ρ_∞)    = q_+ (0) = 1/2 .

With q_∞ =0 this results in 

    ρ_∞    = π√(2 π) +  2π(1+a)  -2 √(2 π) -8 a/a(√(2π)π + 2 π -2 √(2 π) -4)

For a=1 this gives ρ_∞ = 1.443. 
Compared to the very accurate value 1.4371 obtained with a spectral method <cit.> this yields an error of 0.4%.
Maxwell gives  1.2533, which is an error of 12.8%.
The variational method <cit.> gives  1.4245 or an error of 0.9%.




  §.§.§ Coupling conditions


A half space half integration of the coupling conditions gives

    q^i_+ (0)    = -∑_ j=1^n c_ij  q_-^j (0)

which yields

    q^i_∞ -γ^i  +2a/√(2π)ρ^i_∞   = 
    -∑_ j=1^n c_ij( q^j_∞ +γ^j  -2a/√(2π)ρ^j_∞)

Summing these equations gives  ∑_i q^i_∞ =0.
Further, for a uniform node we observe the invariance of

    n-2/nq^i_∞ -γ^i  +2a/√(2π)ρ^i_∞   = 
     n-2/nq^j_∞ -γ^j  +2a/√(2π)ρ^j_∞ .





Moreover, from 

    ρ^i_+ (0) = ∑_ j=1^n c_ijρ_-^j (0)
 
one obtains






    ρ^i_∞+(λ +  √(2π)/a(π-4))γ^i  +  2 q^i_∞/√(2π) a =
      ∑_ j=1^n c_ij( ρ^j_∞+(-λ + √(2π)/a(π-4))γ^j  -  2 q^j_∞/√(2π) a) .
 
Summing these  equations we obtain again ∑_i=1^n γ^i=0.












In the case of a uniform node we can further rearrange, such that






    ρ^i_∞+(n-2/n + √(2π)/π-2)   λγ^i  +  2/√(2π) an-2/nq^i_∞
       =
    ρ^j_∞+(n-2/n +  √(2π)/π-2)λγ^j  +  2/√(2π) an-2/nq^j_∞
 
holds. 
Together with the conditions at infinity we have again 3n conditions for 3 n quantities γ^i, ρ^i_∞, q^i_∞.
Combining (<ref>) and (<ref>) we can eliminate the γ^i
and obtain as invariant









































    ρ
    
    +
    n-2/n n-2/n(π-2)√(2π)+4π-8/√(2π)(π-4)+n-2/n(2π-4)+2√(2π) 1/a 
    q  .
 
or

    ρ
    
    +
    n-2/n1/a 4 + n-2/n√(2 π)/√(2 π)+2 n-2/n q
 


The values of the factor C for the coupling invariants ρ +C q for the present case with unbounded velocities are summarized in Table <ref>. 
 
				

For a=1 and n=3 this  factor is approximately 0.5079.
In contrast, the  Maxwell approximation gives 0.4178.
For a^2 =1/3 the factor is 
6 π/6∼ 0.723 for Maxwell and approximately 0.8797 for the above half moment formula.
This has to be compared to the case of bounded velocities discussed in the last section, where the factor has been determined as 
0.7313 for the half-moment-  and 0.6666 for the Maxwell approximation.
Thus, depending on the underlying kinetic model different coupling conditions are obtained for the same macroscopic  equation.




§ NUMERICAL RESULTS


In this section we compare the numerical results of the different models on networks.
The solutions of the  kinetic equation (<ref>) are compared with the half-moment approximation (<ref>) and the macroscopic wave equation (<ref>) with different coupling conditions  (<ref>) (<ref>) (<ref>)
. 


The networks are composed of coupled edges, each arc is given by an interval x ∈ [0,1], which is discretized with 400 spatial cells if not otherwise stated.
In the kinetic model the velocity domain  [-1,1] is discretized with 400 cells and we choose ϵ = 0.001 if not otherwise stated.

For the advective part of the equations we use an upwind scheme. 
The source term in the kinetic and half moment equations is approximated with the implicit Euler method.
We note that for the wave equation the upwind scheme yields the exact solution by choosing the CFL number equal to 1.

At the nodes we consider in addition to the coupling conditions discussed above for comparison a coupling based on the assumption of equal density at the node, see  <cit.>

    ρ^i   =ρ^j      i≠ j, i,j=1,2,3
    ∑_i=1^3q_i   =0 .


In general, at the outer boundaries  of the network, boundary conditions have to be imposed.
For the kinetic problem the values for the ingoing velocities have to be prescribed as mentioned before. For the half-moment approximation the natural boundary conditions are given by integrating the kinetic conditions and using ρ_+ and q_+ at left boundaries and 
 ρ_- and q_- at right  boundaries.



For the wave equations with full moment,  Maxwell and half-moment conditions
we use the corresponding approximations discussed above to provide boundary values for the macroscopic equations. 

We note that for the boundary values at the  right end of the edges we apply the same procedure as detailed in the previous sections for the left end, reversing  the spatial orientation in the half space problem.
More precisely,  a given kinetic inflow ℓ(v) for v∈[-1,0] at the right boundary leads to the following boundary data.
In the  kinetic model we directly impose f(x=1,v) = ℓ(v) v∈[-1,0], for the  half moment model the quantities ρ^-=<ℓ(v)>_- and q^-=<vℓ(v)>_- are fixed.
In case of the wave equation, the right going Riemann invariant q+aρ=C is given by the inner data.
The remaining information is given by the following  approximations
	
    wave Maxwell          
    		<vℓ(v)>_-= -ρ/4+1/2q
     wave half-moment         {[                <ℓ(v)>_- =1/2ρ - 3/4q +-3a-1/2γ,;               <vℓ(v)>_-     =-1/4ρ + 1/2q +a/4γ ].
    		
     wave full-moment         
    		<ℓ(v)>_-=1/2ρ - 3/4q ,

	which can be obtained by revisiting section <ref> with x∈[-∞,0].

	These are used to prescribe the value of the left going characteristic q-aρ.




 §.§ Tripod network

In a first example, we consider a tripod network with initial conditions
f^1(x,v) = 1/2,f^2(x,v) = 1/3 and f^3(x,v) = 0.
The corresponding macroscopic states are
(ρ^1,q^1) = (1,0),
(ρ^2,q^2) = (2/3,0) and 
(ρ^3,q^3) = (0,0).
We use free boundary conditions at the exterior boundaries.
The computational time is chosen such that the waves generated at the node do not reach the exterior boundaries.

In Figure <ref> we compare the kinetic, the half-moment
and the wave equation with coupling conditions given by the 
Maxwell, the half-moment and the  full-moment approach and the assumption of equal density
at time  T = 1.
		
	We observe first, that the half moment model gives a very accurate approximation of the kinetic equation.
	Considering the wave equation with different coupling conditions one observes that the interior state is very accurately  approximated by the half moment coupling conditions. 
	Also the Maxwell approximation provides a  good approximation in this case.
	Note that a boundary layer is appearing at the node in the kinetic model, see Figure <ref> for a magnification of the situation on edge 1. 
		
			
	Moreover, we investigate the evolution of the total entropy in the network, i.e.
	
    ∑_i=1^31/2∫( (ρ^i)^2+ 1/a^2(q^i)^2 ) dx.

	Initially, the total entropy at t=0 is equal to 0.722222.

	In this case we use a very fine grid with 30000 spatial cells for all models and 400 cells in  velocity space for the kinetic equation.
	In Table <ref> the value of the total entropy at time T=0.1 is shown for the different coupling conditions together with the half moment and the kinetic solution for comparison.
	One observes the very accurate approximation given by the wave equation with half moment coupling conditions.
	Note, that even for the equal density a small amount of  entropy  is lost. 
	This is caused by the numerical diffusion in the very last time step, since we have Δ t<Δ x/a to hit the final time.
						
In Table <ref> the value of the total entropy at time T=0.1 is shown for kinetic and half moment equations for different spatial discretizations. 
					

Finally, we investigate the kinetic equation for different values of ϵ.
For  ϵ→ 0 the kinetic and the half moment model are very well approximated by the solution of the wave equation with half moment coupling conditions, see Figure <ref> and the corresponding magnification at the node in Figure <ref>.

	

	
	


 §.§ Diamond network
	
As a second example we consider a more complicated network, see Figure <ref>, as, for example, studied in <cit.> for the wave equation.

		
As initial conditions for the kinetic equation we choose f^1(x,v) = 1,f^2(x,v) = 5/6 and f^j(x,v) = 1/2 for j=3,…,7,
which corresponds to macroscopic densities ρ^1=2, ρ^2=5/3 and ρ^j=1 for j=3,…,7 and fluxes q^j = 0 j=1,…,7.
These data are also prescribed at the two outer boundaries, i.e. k^1(v)=1, v∈ [0,1] and ℓ^7(v)=1/2, v∈ [-1,0].
Boundary conditions for the wave equation with full moment, Maxwell and half moment conditions are derived as detailed above.
In case of the equal density conditions, we determine the ingoing characteristic using ρ=1, q=0 at the E_1-boundary and ρ=1/2, q=0 at the E_7-boundary.

In Figure <ref> the density ρ^4 on edge 4 is displayed at time t=3 and t=10.
As before, we observe a good agreement of the half moment coupling with the kinetic and half moment model. 
Also the Maxwell approximation is relatively close to the kinetic results.

	
The states of the full moment coupling and the equal density coupling deviate remarkably from the kinetic results.
In Figure <ref> on the left, ρ^4 at time t=50 is shown. 
All models have reached a stationary state except the equal density coupling.
Since the entropy is conserved, a set of waves remains trapped in the network, oscillating back and forth.
The stationary states of the models with entropy losses almost coincide.
	
On the right hand side of Figure <ref> the evolution of the total entropy over time is plotted. 
The total entropy is increasing due to inflow from the boundaries. 
All models saturate at a certain level, but we again observe a deviation of the equal density and the full moment coupling compared to the  accurate  Maxwell coupling. In this situation the results for the  half-moment coupling, the half moment model and the kinetic model coincide.

 


§ CONCLUSION AND OUTLOOK


In this work we have derived coupling conditions for the wave equation on a network based  
on coupling conditions for an underlying  kinetic BGK type model via a layer analysis  of the situation near the nodes.  
The presentation in this work  includes a new half-moment approximation of  the kinetic  half-space problem.
The general approach can be extended to more  complicated problems like linearized Euler equations or kinetic based coupling conditions for nonlinear problems like Burgers and Lighthill-Whitham type equations which will be considered in a forthcoming publication.






§  ACKNOWLEDGMENT
 

The first author is supported by the Deutsche  (DFG) grant BO 4768/1. The second author by DFG KL 1300/26.
Moreover, funding by the DFG
within the RTG 1932 "Stochastic Models for Innovations in the Engineering
Sciences" is gratefully acknowledged.

 






10
siam
ALM10 S. Avdonin, G. Leugering, V. Mikhaylov, On an inverse problem for tree-like networks of elastic strings, Journal of Applied Mathematics and Mechanics (ZAMM), 90, 2, 136–150, 2010

BHK06b
M. Banda, M. Herty, A. Klar, Coupling conditions for gas networks governed by the isothermal Euler equations, NHM 1(2), 295-314, 2006

BHK06a
M. Banda, M. Herty, A. Klar, Gas flow in pipeline networks, NHM 1(1), 41-56, 2006


BSS84 C. Bardos, R. Santos, and R Sentis, Diffusion approximation and computation of the critical size, Trans. Amer. Math. Soc. 284, 2, 617-649, 1984
BLP79 A. Bensoussan, J.L. Lions, and G.C. Papanicolaou, Boundary-layers and homogenization of transport processes, J. Publ. RIMS Kyoto Univ. 15, 53-157, 1979

BCG10
R. Borsche, R. Colombo, M. Garavello. On the coupling of systems of hyperbolic conservation laws with ordinary differential equations. Nonlinearity 23, 11,2749, 2010
BGKS14
R. Borsche, S. Göttlich, A. Klar, and P. Schillen.
The scalar Keller-Segel model on networks.
Math. Models Methods Appl. Sci., 24, 2, 221–247, 2014.
BKKP16
R. Borsche, J. Kall, A. Klar, and T.N.H. Pham.
Kinetic and related macroscopic models for chemotaxis on networks.
Mathematical Models and Methods in Applied Sciences,
26, 6, 1219–1242, 2016.

BNR14
G. Bretti, R. Natalini, and M. Ribot.
A hyperbolic model of chemotaxis on a network: a numerical study.
ESAIM: M2AN, 48, 1, 231–258, 2014.
CC17
F. Camilli and L. Corrias.
Parabolic models for chemotaxis on weighted networks.
to appear in Journal de Mathematiques Pures et Appliquees,
2017.

C69 C. Cercignani, A Variational Principle for Boundary Value Problems,
 J. of Stat. Phys. 1,  2, 1969
C88 C. Cercignani, The Boltzmann Equation and its Applications,
Springer, 1988
CLT14 I.K. Chen, T.P. Liu, and S. Takata, Boundary singularity for thermal transpiration problem of the linearized Boltzmann equation, Arch. Ration. Mech. Anal. 212, 2, 575–595, 2014. 
CGP05
G. M. Coclite, M. Garavello, and B. Piccoli, Traffic flow on a road
	network, SIAM J. Math. Anal., 36 , 1862–1886 2005.
CHS08
R. Colombo, M. Herty, V.Sachers - On 2 × 2 conservation laws at a junction - SIAM J. Math. Anal. 40, 2, 2008.


Coron F. Coron, Computation of the Asymptotic States for Linear Halfspace
Problems, TTSP 19, 2,  89, 1990
CGS88 F. Coron, F. Golse, C. Sulem, A Classification of Well-posed Kinetic
Layer Problems, CPAM,  41,  409, 1988
DZ01 R. Dager, E. Zuazua, Controllability of tree-shaped networks of vibrating strings, C. R. Acad. Sci. Paris, 332, 1087–1092, 2001
DL02 B. Dubroca, A. Klar, Half Moment closure for radiative transfer equations, J. Comp. Phys., 180 (2), 584-596, 2002
EK16H. Egger, Thomas Kugler, Damped wave systems on networks: Exponential stability and uniform approximations
https://arxiv.org/abs/1605.03066 (2016)
FT15
L. Fermo and A. Tosin.
A fully-discrete-state kinetic theory approach to traffic flow on
road networks.
Math. Models Methods Appl. Sci., 25,3, 423–461, 2015.


G10
M. Garavello,
A review of conservation laws on networks NHM 5, 3, 565 - 581, 2010      

G92 F. Golse, Knudsen Layers from a Computational Viewpoint, TTSP 21,3 ,  211, 1992
G08 F. Golse, Analysis of the boundary layer equation in the kinetic theory of gases, Bull. Inst. Math. Acad. Sin. 3, 1, 211-242, 2008
GK95 F. Golse, A. Klar, Numerical Method for Computing Asymptotic States and Outgoing Distributions for a Kinetic Linear Half Space Problem, J. Stat. Phys. 80, 5-6, 1033-1061, 1995
GMP W. Greenberg, C. van der Mee, V. Protopopescu, Boundary Value
Problems in Abstract Kinetic Theory, Birkhäuser, 1987
GMS09
G. Guerra, F. Marcellini, V. Schleper, Balance laws with integrable unbounded sources. SIAM Journal on Mathematical Analysis 41, 3, 1164–1189, 2009


GHKLS11
M. Gugat, M. Herty, A. Klar, G. Leugering, V. Schleper, Well–posedness of networked hyperbolic systems of balance laws, International Series of Numerical Mathematics, 160, 175-198, Springer, 2011





HKP07
M. Herty, A. Klar and B. Piccoli.
Existence of  solutions for supply chain models based on  partial
differential equations.
SIAM J. Math. Anal. 39, 1, 160-173, 2007.
HM09
M. Herty and S. Moutari.
A macro-kinetic hybrid model for traffic flow on road networks.
Comput. Methods Appl. Math., 9, 3,238–252, 2009.

LLS162 Q. Li, J. Lu, and W. Sun, Half-space kinetic equations with general boundary conditions, Math. Comp. 86, 1269-1301,  2017
LLS16 Q. Li,J.Lu,W. Sun, A convergent method for linear half-space kinetic equations, to appear in  ESAIM: M2AN, 2017


LF67  S.K. Loyalka, J.H. Ferziger, Model Dependence of the Slip Coefficient, Phys. Fluids 10,  8, 1967
LF81  S.K. Loyalka, Approximate Method in the Kinetic Theory, Phys. Fluids
11,  14, 1971
M47 R. E. Marshak, Note on the spherical harmonic method as applied to the Milne problem for a sphere, Phys. Rev. 71, 443-446, 1947
M J.C. Maxwell, Phil. Trans. Roy. Soc. I, Appendix, 1879; reprinted in
The Scientific Papers of J.C.Maxwell, Dover, New York, 1965
ST81  C.E. Siewert, J.R. Thomas, Strong Evaporation into a Half
Space I, Z. Angew. Math. Physik, 32, 421, 1981
SO Y. Sone, Y. Onishi, Kinetic Theory of Evaporation and Condensation,
Hydrodynamic Equation and Slip Boundary Condition, J. Phys. Soc. of Japan, 
44,  6,  1981, 1978


T09E.F. Toro Riemann solvers and numerical methods for fluid dynamics, Springer, 2009

UTY03 S. Ukai, T. Yang, and S.-H. Yu, Nonlinear boundary layers of the Boltzmann equation. I. Existence, Comm. Math. Phys. 236,  3, 373-393, 2003


 
VZ09 J. Valein, E. Zuazua, Stabilization of the Wave Equation on 1-d Networks, SIAM Journal on Control and Optimization, 48, 4, 2771-2797, 2009 



  
