{'fraction_non_alphanumeric': 0.03699328704370625, 'fraction_numerical': 0.010919358091561266, 'mean_word_length': 4.569640564826701, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Solving the stationary Liouville equation via a boundary element method

Intensity distributions of linear wave fields are, in the high frequency limit, often approximated in terms of flow or transport equations in phase space. Common techniques for solving the flow equations for both time dependent and stationary problems are ray tracing or level set methods. In the context of predicting the vibro-acoustic response of complex engineering structures, reduced ray tracing methods such as Statistical Energy Analysis or variants thereof have found widespread applications. Starting directly from the stationary Liouville equation, we develop a boundary element method for solving the transport equations for complex multi-component structures. The method, which is an improved version of the Dynamical Energy Analysis technique introduced recently by the authors, interpolates between standard statistical energy analysis and full ray tracing, containing both of these methods as limiting cases. We demonstrate that the method can be used to efficiently deal with complex large scale problems giving good approximations of the energy distribution when compared to exact solutions of the underlying wave equation.


Introduction
Many phenomena in physics and engineering can be accurately described or modelled in terms of linear wave equations. A large range of numerical methods has been developed for solving wave problems, often making heavy use of the linearity of the underlying PDEs leading in general to large, finite matrix equations. Popular tools include finite element or finite volume methods, boundary element methods and various spectral methods. There are, however, basic limitations when approximating the solutions of the wave equations directly: the size of the associated linear system increases with decreasing wavelength and numerical schemes become inefficient when the local wavelengths are orders of magnitude smaller than typical dimensions of the physical system.
In order to overcome this limitation, a range of high-frequency/ short wave length methods have been devised such as the WKB or eikonal approach or semiclassical methods developed in the context of quantum chaos, see [1,2,3] for overviews. In particular, phase information is obtained by solving the Hamilton-Jakobi equations for the action fields, whereas amplitude information is calculated using transport equations. The actual methods can be quite intricate demanding an intimate knowledge of the classical phase space dynamics given by the Hamilton equations of motion for the underlying ray or classical dynamics. If one is prepared to omit phase information altogether, the approximations amount to solving the transport or Liouville equation for phase space densities.
In this paper, we will focus on efficient numerical methods for solving the Liouville equation (LE) for stationary problems. We will show how the solutions of the LE are related to both ray tracing (see for example, [4,5,6]) and statistical energy analysis (SEA) [7,8,9], a popular method in the engineering community for solving vibro-acoustic problems in the high-frequency limit.
The Liouville equation is an example of a PDE describing a conservation law [10]. Typical solution methods are based on characteristics or ray tracing. Ray tracing is the method of choice when considering propagation over relatively short times or length intervals. Prime examples include weather forecasting, where data updates lead to regular re-initialisations of the density distributions [11], or room acoustics, where only very few reflections need to be considered [5]. More sophisticated methods related to tracking the time-dynamics of interfaces in phase space such as moment methods and level set methods [12,13,14] have been developed only relatively recently, finding applications in acoustics, seismology and computer imaging. For an excellent recent overview, see Runborg [2]. Numerical methods based on solving the time dependent Liouville equation directly on a mesh using finite volume or transfer operator methods, also referred to as Ulam methods, have been considered in [15] and more recently in [16,17]. Ray based and tracking methods often become inefficient when considering frequency domain wave problems in bounded domains, for example, determining the wave field in a finite size cavity driven by a continuous, monochromatic excitation. Here, multiple reflections of the rays and complicated folding patterns of the associated level-surfaces often lead to an exponential increase in the number of branches that need to be considered.
In the engineering literature, these problems have been circumvented by subdividing the structure into a set of substructures. Assuming ergodicity of the underlying ray dynamics and quasiequilibrium conditions in each of the subsystems, that is, assuming that the density in each subsystem is approximately constant, greatly simplified SEA equations based only on coupling constants between subsystems can be set up. The disadvantage of SEA is that the underlying assumptions are often hard to verify a priori or are only justified when an additional averaging over 'equivalent' subsystems is considered. These shortcomings have been addressed by Langley [18,19] and more recently in a series of papers by Le Bot [20,21,22]. A computational tool based on a Frobenius-Perron operator approach systematically interpolating between SEA and full ray tracing has been introduced first in [23] and further analysed in [24]. Its name, Dynamical Energy Analysis (DEA), points at the similarities with SEA but stresses at the same time the importance of dynamical correlations present in the ray-densities. Thus, in DEA, the solution in a subsystem of the structure is typically represented by a vector; this is in contrast to SEA where the solution in each subsystem is represented by a single number, the mean ray density in each subsystem.
The implementation of DEA as presented in [23,24] corresponds to a spectral boundary integral method; the integral equations are expanded in orthogonal basis sets and the resulting matrix equations are solved for the coefficients of the basis expansion of the solution. It has been shown in [24] that the method works even in situations where a naive application of SEA fails. However, spectral methods can lead to slow convergence for the matrix elements in the DEA matrix, in particular if the boundary is non-smooth. In this paper, we show that using a boundary element method (BEM) for the spatial variable and a basis function expansion in the momentum component leads to large efficiency gains. We will furthermore present a more in depth derivation of the actual boundary integral formulation.
The paper is structured as follows: in Sec. 2, we sketch the derivation of the Liouville transport equation from the underlying wave equation; we will then derive the boundary integral formulation in Sec. 3 and give a description of the numerical implementation of DEA using the BEM/ basis function approach in Sec. 4. Numerical results will then be presented in the Sec. 5.
2. From wave dynamics to phase space densities

A time dependent formulation
In the following, we sketch the connection between linear wave equations and the Liouville equation via the Eikonal ansatz. We develop the theory in the time domain first and move to a stationary formulation in the next section. We will restrict the discussion to the Helmholtz equation of acoustics; other wave problems such as electromagnetic waves, linear elasticity or Schrödinger's equation can be treated along the same line with the obvious modifications. For a more in depth discussion of the Eikonal approach, see for example [2]. A nice derivation of the connection between wave transport equations and the Liouville equation based on Wigner transformation techniques has been given in [25].
We start from the wave equation for linear acoustics, that is, seeking solutions u(r, t) in some domain Ω with boundary Γ (u may express acoustic pressure in appropriate units). Furthermore, c = c(r) ≥ 0 is the local wave speed. We assume in what follows that the wave equation is not explicitly time dependent and that the walls act as passive elements, that is, we have Dirichlet or Neumann boundary conditions, u(s) = 0 or ∂u/∂n(s) = 0 for s ∈ Γ. Solutions are then determined uniquely in terms of initial conditions u(r, 0) = u 0 (r) and ∂ t u(r, 0) = u 1 (r) at t = 0. We now make the usual Eikonal ansatz with real valued functions A = A(r, t), φ = φ(r, t) and ω acting as the large parameter here. Plugging (2) into (1) and collecting terms of the order ω 2 , we arrive at the Eikonal equation Terms of order ω yield a transport equation we can write Eq. (4) in the form of a conservation law, that is, with "velocity field" v(r, t) = −c 2 ∇φ/ ∂φ ∂t (6) determined through the phase function φ(r, t) alone. Solutions of the Eikonal equation (3) can readily be found using characteristics. Without restricting the discussion we choose the "+" sign in (3); see also comments further down in the text. Setting H(r, p) = c|p| with momentum p = ∇φ , one obtains solutions φ(r, t) along trajectories in phase space X(t) = (r(t), p(t)) by solving Hamilton's equations of motion, that is,ṙ Note that in deriving the Eikonal Eq. (3) and transport Eq. (4) we omitted terms of order O(ω 0 ); neglecting the coupling terms between Eqs. (3), (4) leads to first order differential equations which can formally be solved by using only one of the initial conditions, that is, either u 0 (r) or u 1 (r). The ODEs (8) are thus solved by using, for example, initial conditions X(0) = (r 0 , p 0 = ∇φ 0 (r 0 )). A solution for the phase φ is then obtained by integration, that is, where the integral is taken along X(t) = ϕ t (X 0 ). Here ϕ t (X 0 ) is the phase space flow map obtained from the pair of ODEs (8). In particular we integrate over trajectories starting at t = 0 at some point r 0 with momentum p 0 = ∇φ 0 (r 0 ) and reaching the final point r at time t. Here, φ 0 is the initial phase corresponding to the initial conditions u 0 (r 0 ) = A 0 (r 0 ) exp(iωφ 0 (r 0 )) at t = 0. In general, there will be several such initial points from which a trajectory reaches the point r at time t. That is, the solution to the eikonal equation (3) at a point r has a multi-sheet structure with a set of N = N (r, t) branches, each giving rise to a phase The integrals are taken over paths X j , j = 1, ..., N , reaching r in time t with X j (0) = (r j , ∇φ 0 (r j )). It is this multi-sheet structure which makes solving Eq. (3) difficult in practice [2]. In particular, the number of branches N (t) can increase rapidly with t and for chaotic ray dynamics N (t) actually increases exponentially [1]. We would like to remark that choosing the "−" sign in (3) leads to a sign-change in (8) and gives time-reversed ray trajectories obtained through the time transformation t → −t. The complete solution to the Eikonal equation is thus given by the set of characteristics travelling both in positive and negative time.
We turn now to the solution of the transport Eq. (5); this equation is driven by the solution φ(r, t) of the Eikonal Eq. (3). It thus suffers in principle from the same multi-valuedness problem as the original equation. Note, however, that the solutions of the ray equations (8) in phase space are determined uniquely for a given initial condition; it is the projection onto r -space which leads to the multi-valuedness of the phase solution. It is thus advantageous to solve the transport equation in phase space. The corresponding equation is the Liouville equation (LE) where {· , ·} are the Poison brackets, and H(r, p) is the classical Hamilton function (7). The phase space densityρ indeed solves Eq. (10) if the branches φ j (r, t) andρ j (r, t) obey the Eikonal and transport equations (3), (5), respectively. (This result can be obtained by inserting (11) into (10), see also [2].) A formal solution of the LE is obtained in terms of a Frobenius-Perron operator where the last relation is valid for Hamiltonian flows; this relation is used extensively in the ray tracing approach. The solution of the wave equation (1) is now formally obtained by summing over the solution branches in (9), that is, where ν j are Maslov indices [1] due to caustics along the ray X j (t). As we are not interested in phases in what follows, we will not discuss this issue any further here. An asymptotic treatment of the wave equation using the Eikonal ansatz needs to keep track of the link between phases and amplitudes along the different solution branches. Ray tracing methods are often the only way forward when finding solutions using the Eikonal ansatz. This may be an efficient method for short times and in cases where only a few reflections occur. However, the rapid increase in the number of paths makes ray tracing methods quite cumbersome for large times. The advantages behind the original idea of the Eikonal ansatz, namely replacing a wave equation with highly oscillatory solutions by differential equations for slowly varying quantities such as the amplitude and phase, are lost.
In many applications it is, however, sufficient to find the solution of the transport equations for a given set of initial conditions. The densityρ(X, t) contains a great deal of information about the wave solution, and this may often be sufficient for the required purpose. Classical ray tracing, that is ray tracing without tracking phase information, is indeed used extensively in room acoustics.

Stationary formulation
In what follows, we will consider wave problems with time-harmonic forcing at a frequency ω. We can then separate out time by setting Let us consider problems with a forced excitation by a point source. The corresponding inhomogeneous wave equation then has the form with the source point at r = r 0 . The solution G(r, r 0 ) is the Green function and arbitrary force distributions can be recovered.
For an Eikonal approach, we proceed by splitting the solution into a homogeneous and an inhomogeneous (or direct) solution in the form where G 0 is the free Green function and G h solves the homogeneous Helmholtz Eq.
with appropriate boundary conditions. (For Dirichlet BC, that is G(s) = 0, we obtain G h (s) = −G 0 (s) for s ∈ Γ). One finds for two-dimensional problems with c = const, where H 1 0 is the 0th order Hankel function of the first kind. Proceeding with an Eikonal ansatz, we set G h = A exp(iωψ); we have then φ(r, t) = ψ(r) + t and ∂A ∂t = 0.
The Eikonal and transport Eqs. (3), (5) reduce to c|∇ψ| = 1 (18) and We write the boundary conditions in the form In order to obtain an approximate solution of (18) at a point r ∈ Ω one needs to find those rays starting from a boundary point s ∈ Γ with initial conditions X 0 = (s, p = ∇ψ 0 (s)) which pass through the point r following the equations of motion (8). The phase function ψ is then obtained by integrating along this trajectory with an initial condition ψ(s) = ψ 0 (s); the phase function ψ again has a multi-sheet structure and there are in general infinitely many sheets intersecting with the plane r = const in phase space. (A phase space description of the sheets can be given in terms of Lagrangian manifolds, see for example [26].) The various branches are calculated according to integrating along a trajectory with initial condition X j = (s j , p j = ∇ψ 0 (s j )). The transport equation (19) is solved with initial conditions where j 0 (s) is the incoming flux normal to the boundary. The total solution is then approximatively obtained as where µ j are again Maslov phases and |A j | is the amplitude transported along the trajectory from X j to r; for alternative derivations of Eqs. (13), (21), see [1]. In applications, one is often interested in wave intensities (such as energy densities); the contributions arising from reflections from the boundary have the form Neglecting the oscillatory terms, which is justified after averaging over a small ω range or over an ensemble of similar systems, one arrives at the ray tracing approximation where ρ(r, p) is a solution of the stationary Liouville equation with and boundary conditions given by (20) together with Snell's law of reflection. Note that the total intensity |G(r)| 2 for r ∈ Ω may be estimated simply by adding |G 0 (r)| 2 to ρ(r). It is clear from the description above that a ray tracing approach is even more cumbersome in the stationary case than in the time-dependent case. The stationary solution is constructed from infinitely many branches which originate from reflections at boundaries. If we neglect phase information, however, we can reconstruct the phase space densities using the Frobenius-Perron operator, Eq. (12) [23,24]. In the next section we will introduce an efficient algorithm for determining solutions of the stationary Liouville equation even in the presence of multiple reflections, thus overcoming the limitations of the ray tracing approaches sketched above.

Boundary integral formulation of the Liouville equation
In order to solve the stationary Liouville equation (24), the ray tracing integral formula, Eq. (12), is rewritten in a boundary integral form using a boundary mapping technique. The boundary mapping procedure involves first mapping the ray density emanating continuously from the source onto the boundary Γ and then constructing a Frobenius-Perron operator for the boundary map.
We now give an explicit derivation of the initial density on the boundary arising from Eq. (14) for Dirichlet boundary conditions. The ray density emanating from a source point r 0 ∈ Ω is given by the square of the amplitude of the free space Green function. For Ω ⊂ R 2 , it is given by (16) and we have the following high frequency asymptotic form Lifting the ray density to the full phase space, Eq. (25) corresponds to rays moving out from the source in all directions according to where p 0 = |p|(r − r 0 )/|r − r 0 |. Note that from Eq. (18), the ray dynamics take place on the 'energy' manifold c|p| = 1. To obtain the source density ρ Γ 0 on the boundary Γ, we set where s parameterizes Γ and p s ∈ B d−1 |p| denotes the momentum component tangential to Γ at s for fixed H(X) = c|p| = 1; here B d−1 |p| is an open ball in R d−1 of radius |p| and centre s. Using local coordinates on the boundary, that is, . Combining these identities with (26) and (27) yields where θ(r(s), r 0 ) is the angle between the normal to Γ at r(s) and the ray from r 0 to r(s). The resulting boundary layer density ρ Γ 0 is equivalent to a source density on the boundary producing the same ray field in the interior as the original source field after one reflection.
The operator equivalent to the Frobenius-Perron operator (12), but now acting on boundary densities ρ Γ , is described in terms of a boundary integral operator B, where X s = (s, p s ), Y s = (s , p s ) and γ is the invertible boundary map. Also, ∂P = Γ × B d−1 |p| is the phase space on the boundary at fixed "energy" H(X) = c|p|. The weight w contains factors due to reflection/ transmission at interfaces within Ω and a damping term of the form exp (−µL(s, s )), where L(s, s ) is the length of the trajectory between s and s . Note that, depending on the modelling assumptions, the factors in w may depend on ω. A similar weight term also needs to be applied in Eq. (28) since rays emanating from the source will also be subject to dissipation and possible reflection/ transmission at interfaces. We assume uniform damping in the interior, that is, the damping coefficient µ does not depend on r. The damping term does then not alter the underlying ray dynamics [27]. Note that convexity is assumed to ensure γ is well defined; nonconvex regions can be handled by introducing a cut-off function in the shadow zone as in Ref. [21] or by subdividing the regions further.
The stationary density on the boundary induced by the initial boundary distribution ρ Γ 0 (X s , ω) can then be obtained using 8 where B n contains trajectories undergoing n reflections at the boundary. The resulting density distribution on the boundary ρ Γ (X s , ω) is then mapped back into the interior region. Using Eq. (23), one can relate ρ Γ to a wave energy density after projecting onto r space, that is for Ω ⊂ R 2 , ρ(r, ω) = ρ(r, p)dp, where (|p|, Θ) is a polar coordinate system for p. One finally obtains Now performing a change of variables from Θ to s yields where r s ∈ Γ is the point on Γ in Cartesian coordinates corresponding to s. For systems with damping we must again apply a factor of the form exp (−µL(s, s )) to the final density evaluation (33). The total intensity |G(r, r 0 ; ω)| 2 may be estimated simply by adding |G 0 (r, r 0 ; ω)| 2 to the value of ρ(r, ω) given by (33). We remark in passing that our approach is akin to the infinitesimal generator technique presented in [17] which aims at estimating the long term behaviour of the dynamics by numerically calculating the spectral properties of the generator of the Frobenius-Perron operator and extrapolating to long times.

Numerical implementation: basis approximations and boundary element techniques
The long term dynamics is contained in the operator (I − B) −1 and standard properties of Frobenius-Perron operators ensure that the sum over n in Eq. (30) converges for dissipative or open systems. In order to evaluate (I − B) −1 , a finite dimensional approximation of the operator B must be constructed. In Refs. [23,24], basis expansions have been applied in both position and momentum coordinates, which is straightforward to implement for Ω ⊂ R 2 using univariate expansions in each argument. In particular, in Ref. [23], a Fourier basis was employed globally on each subsystem in both position and momentum. It is well known that the convergence rate of a Fourier basis expansion depends on the smoothness of the approximant and also requires periodic boundary conditions [28]. The periodicity requirement will be fulfilled due to periodicity of the boundary curve and that the ray density corresponding to p s = ±|p| is zero giving periodicity in the momentum coordinate. However, the smoothness criterion often fails if, for example, corners are present in Γ; the initial density due to a point source is not then differentiable at the corners.
In Ref. [24] it was suggested to split up the spatial approximation at corners in Γ and apply a basis approximation with no requirement for periodic boundary conditions in order to attain its optimal convergence properties. An approximation using a Chebyshev basis was implemented since Chebyshev polynomials may be computed simply in terms of trigonometric functions and have similar convergence properties to Fourier basis approximations without the periodicity requirement. Such an approach also has a natural and optimal (in terms of polynomial degrees that can be integrated exactly) quadrature choice for the associated integrals, namely, Gauss-Chebyshev quadrature. The improvement in computational efficiency brought about by this approach opened up the method to larger scale multi-component problems.
In this work we go further along this route by subdividing the boundary using a mesh. This introduces greater flexibility into the numerical method making it possible to both increase the polynomial order of the basis functions and to refine the mesh for the spatial coordinates. The advantage of refining the mesh is that additional degrees of freedom are introduced to the model whilst maintaining relatively low quadrature costs due to the lower order of the spatial basis approximation. As such it is preferable to choose a basis which is orthogonal in the standard L 2 inner product. (From previous experience using a Chebyshev basis, where the associated orthogonal inner product is weighted, we found that low order approximations were very poor). Due to the mesh discretization procedure adopted here, the boundary conditions are not periodic and using a Fourier basis is not advisable; hence an alternative L 2 −orthogonal basis choice is proposed, namely, Legendre polynomials.
The overall form of the proposed numerical scheme in 2D is therefore based on a Legendre polynomial approximation in both position and momentum augmented with a boundary mesh through which the spatial approximation may be localised. For an extension of the method to 3D cavities, see [29]. Recall that p s ∈ (−|p|, |p|) and denote byp s a re-scaling of p s to (−1, 1). Let us also denoteP m (p s ) = 1/|p|P m (p s ). An important special case of the scheme described above is when we fix the spatial basis approximation at zero (constant) order and only refine the mesh, which essentially results in a scaled piecewise constant boundary element approximation. This type of approximation is also often referred to as Ulam's method [15,16], although this would involve performing such an approximation in full phase space, rather than just in its spatial component. The implementation of separate meshes and basis approximations in both position and momentum would, however, complicate the procedure. In addition, our final computations will involve integration over spatial variables only. A spatial mesh therefore has the additional benefit of localising the integration regions, permitting the use of lower order quadrature rules.
The overall approximation is then of the form where N s and N p are the orders of the position and momentum basis expansions, respectively, and n is the number of elements in the mesh. Alsô wheres α parameterizes the α th element and is scaled to take values in (−1, 1), χ α is the characteristic function for element α and A α is the length of element α, α = 1, .., n. An analogous approximation is also made for ρ Γ 0 , the initial density, and the values of ρ α,l,m are to be determined by solving the resulting linear system. Note that the special case of piecewise constant boundary elements is reached by setting N s = 0 above. The other limiting case, N p = 0, corresponds to a 'diffusive' reflection approximation, also denoted as Lambert's law [5], which is equivalent to the assumption that incoming rays are scattered uniformly in all directions after reflecting from a boundary. Equation (30) then becomes equivalent to the integral equations discussed in [21,30]. Note that the case N p = N s = 0 and n = 1 corresponds to an SEA approximation for a single subsystem (n equals the number of substructures for multi-component systems ) [23].
The matrix approximation B of the linear operator B (29) is given by taking the inner product < ·, · > in L 2 (∂P): Here we write γ = (γ s , γ p ), to denote the splitting of the position and momentum parts of the boundary map. The boundary map can be very difficult to obtain in general, and hence we write the operator in terms of trajectories with fixed start and end points s and s as follows B m+1+Np(l+Ns(α−1)),b+1+Np(a+Ns(β−1)) = (2m + 1)(2l + 1) 4 Γ Γ w(Y s , ω)P m (p s (s, s ))P α,l (s)P b (p s (s, s ))P β,a (s ) ∂p s ∂s ds ds.
The resulting boundary integral formulation containing a pair of integrals over boundary coordinates bears a resemblance to standard variational Galerkin boundary integral formulations, see for example [31].

Numerical results
In this section we consider a two-dimensional polygonal example whose boundary is meshed by subdividing each side into equidistant sections governed by a mesh parameter ∆x. The number of elements on any given side is computed using the integer part of the side length divided by ∆x. The Jacobian in Eq. (36) is written in the form where n and n are the inward unit normal vectors to Γ at r and r , respectively. In addition r and r are the Cartesian coordinate vectors corresponding to the boundary parameterised coordinates s and s . In order to treat the corner singularities in (37), Gaussian quadrature is employed where end-points are not included as quadrature points. The convergence of the quadrature rules is still slow due to the peak in the integrand at corners. Telles' transformation techniques are employed to speed up the convergence [32]. A five cavity system is considered as shown in Fig. 1 with Dirichlet boundary conditions on the outer boundary. The application of the method to multiple sub-domains is carried out using the same techniques as detailed in [24]. In general the configuration features irregular shaped, well separated pentagonal subsystems. In such a case one expects the wave energy to be fairly evenly distributed within subsystems and the energy levels in each subsystem to be distinct and well separated. The one feature of the configuration where such trends may not hold is due to the direct channel running to the right of the source point leading to a corridor type effect and energy being localised along this channel. Such dynamical features are likely to require higher order applications of our numerical algorithm to adequately resolve the behaviour.
Energy distributions have been studied as a function of the frequency with a hysteretic damping factor η = 0.01, where µ = ωη/(2c) throughout Ω. Here and in the remainder of this section the subsystems (or cavities) are numbered 1, 2, ..., 5 from left to right. The wave speed is set to unity for simplicity. Note that the method allows c (and thus µ) to vary between subsystems [24], but we do not make use of this in the examples considered here. Therefore the subsystem interface reflection and transmission coefficients appearing in the weight term in (29) are simply 0 and 1, respectively. Fig. 2 shows the energy distribution throughout the five subsystems compared with the results of a boundary element computation for the full wave problem. Explicitly we compute the square of the energy norm The lines show the approximate energy norms given by the Liouville equation model computed for f = (ω/(2π)) = 10, 13, 16, 19, ..., 70 with N p = 10, N s = 0 and n = 627, and hence the spatial approximation has been improved via mesh refinement rather than using higher order Legendre polynomial basis expansions. This has the benefit that the quadrature costs are kept moderate and leads to considerable efficiency gains. The computations here took approximately 2.5 hours per frequency compared with typical times of approximately 53 hours for a full tenth order Chebyshev basis computation using the methods reported in [24]. Here the computations were performed using a single core of a desktop PC with a 2.83 Ghz dual core processor. The circles each show the mean of 22 boundary element method computations, where the average is taken over a small range of frequencies centered around 10, 20,..., 70. The plus signs show the maxima and minima of the boundary element computations. The plot shows a good agreement between the solution of the Liouville equation and the mean of the full wave problem computations. In addition, the Liouville solution always lies within the range of the data for the full wave problem. Fig. 3 shows the convergence of our numerical algorithm in subsystem 5 as the precision of the approximation is increased by refining the mesh and simultaneously increasing the order of the Legendre basis approximation in direction space. The plot shows the relative error in the energy norm G 5 2 compared against the value G e 5 2 computed with N p = 10, N s = 0 and n = 627. The relative error plotted is thus given by ( G e The trend of the energy increasing as the precision of the numerical method is increased is shown by the error values all being positive and is expected since a low order implementation of our numerical scheme will not capture the direct energy channelled to the right of the source and into subsystem 5. Higher order basis functions in direction space are required to capture such directivity patterns in the wave field. In addition the plot shows the approximated solutions converging since the errors are getting closer to zero as higher order approximations are employed. For the case N p = 2, n = 70 one sees the error increasing as the frequency and thus damping are increased. For high damping, the direct channel

15
of energy (before reflection) from the source point will form a more significant part of the overall solution. For the case N p = 4, n = 70, the higher order basis in direction space has been able to model the channelling effect more accurately and the error has become almost uniform across the frequency range. Fig. 4 shows the energy density |G(r, r 0 ; ω)| 2 plotted for r ∈ Ω and f = 15. The upper subplot shows the energy density approximated using a boundary element method for the full wave problem. The central subplot shows the average of 22 full wave solutions at different frequencies 13 < f < 17. One sees that the averaging reduces the influence of small scale fluctuations in the solution and produces a generally smoother looking plot, showing more clearly the influence of the source point. The lower subplot shows the energy density approximated using the Liouville equation model with N p = 10, N s = 0 and n = 627. Many similarities with the central subplot are apparent. For example, the direct energy channel to the right of the source and into subsystem 5, the location of the source and the dip in energy in the second subsystem close to the interface are features present in both plots. Also the energy channelled along the upper boundaries of subsystems 1 and 4 is a common feature. This demonstrates that high frequency wave models based on the Liouville equation are indeed a powerful tool to obtain averaged solutions. Such an average solution can be obtained through averaging over an ensemble of similar structures or over a small frequency range. Applications are abundant in high frequency vibro-acoustics, where small changes in the domain due to, for example, manufacturing tolerances lead to large variations in the response and thus a simulation of any one particular problem will only be of statistical significance. Indeed, replacing highly oscillatory wave problems with smooth ray dynamical problems is often the only practical option in the high frequency regime.

Conclusion
The suitability of the Liouville equation as a model for high frequency wave problems has been discussed. A boundary integral approach has been proposed for its solution using the Frobenius-Perron operator. A numerical solution method has been implemented based on Legendre polynomial approximations in phase space together with a spatial boundary mesh. The proposed method has been verified by comparison with a BEM code for the full wave problem and numerical evidence of convergence has been shown.