Error analysis of an asymptotic preserving dynamical low-rank integrator for the multi-scale radiative transfer equation

Dynamical low-rank algorithm are a class of numerical methods that compute low-rank approximations of dynamical systems. This is accomplished by projecting the dynamics onto a low-dimensional manifold and writing the solution directly in terms of the low-rank factors. The approach has been successfully applied to many types of differential equations. Recently, efficient dynamical low-rank algorithms have been applied to treat kinetic equations, including the Vlasov--Poisson and the Boltzmann equation, where it was demonstrated that the methods are able to capture the low rank structure of the solution and significantly reduce the numerical effort, while often maintaining good accuracy. However, no numerical analysis is currently available. In this paper, we perform an error analysis for a dynamical low-rank algorithm applied to a classical model in kinetic theory, namely the radiative transfer equation. The model used here includes a small parameter, the Knudsen number. This setting is particularly interesting since the solution is known to be rank one in certain regimes. We will prove that the scheme dynamically and automatically captures the low-rank structure of the solution, and preserves the fluid limit on the numerical level. This work thus serves as the first mathematical error analysis for a dynamical low rank approximation applied to a kinetic problem.



page 1

page 2

page 3

page 4


An efficient dynamical low-rank algorithm for the Boltzmann-BGK equation close to the compressible viscous flow regime

It has recently been demonstrated that dynamical low-rank algorithms can...

An asymptotic-preserving dynamical low-rank method for the multi-scale multi-dimensional linear transport equation

We introduce a dynamical low-rank method to reduce the computational com...

Existence of dynamical low-rank approximations to parabolic problems

The existence of weak solutions of dynamical low-rank evolution for para...

Efficient 6D Vlasov simulation using the dynamical low-rank framework Ensign

Running kinetic simulations using grid-based methods is extremely expens...

Efficient dynamical low-rank approximation for the Vlasov-Ampère-Fokker-Planck system

Kinetic equations are difficult to solve numerically due to their high d...

A low-rank power iteration scheme for neutron transport criticality problems

Computing effective eigenvalues for neutron transport often requires a f...

Low-rank lottery tickets: finding efficient low-rank neural networks via matrix differential equations

Neural networks have achieved tremendous success in a large variety of a...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1. Introduction

Kinetic equations are a class of model equations used to describe the statistical behavior of a large number of particles that follow the same physics laws. They have been widely used in many aspects of physics and engineering: the classical Boltzmann equation describes the dynamics of rarefied gases, the radiative transfer equation characterizes the behavior of photons, the neutron transport equation describes the dynamics of neutrons in nuclear reactors, and the Vlasov–Poisson system has been used to describe the motion of plasmas.

These equations, despite there being significant differences for various kinds of particles, share a similar structure. In all situations the dynamics is described in phase space and the solution is thus a distribution function that counts the number of particles at a particular time , location , and velocity . One major reason for kinetic equations being challenging is that they are posed in a higher dimensional space; this is different from most other physical models that describe the corresponding dynamics only on the physical domain, e.g. for a fluid model the solution only depends on .

During the past decade, numerous methods have been proposed to numerically solve kinetic equations and many of them succeed in reducing computational cost without sacrificing accuracy. The main body of work concentrates on developing fast solvers for the collision operators and overcoming the stability requirement via relaxing the time-discretization [24, 22], through either finding fast methods to treat the transport term [17, 44, 5, 12], or implementing the resulting discrete system efficiently [18, 13]

. However, reducing the complexity due to the high dimensionality is largely left unaddressed. The difficulty here is rather clear. To numerically solve kinetic equations, one needs to sample a certain amount of discrete points in each dimension, resulting in a large number of degrees of freedom. The numerical cost, meanwhile, is typically determined mainly by the total number of degrees of freedom.

This viewpoint was challenged in [42] (based on earlier work on very high dimensional systems in quantum mechanics). There the authors abandon the traditional approach. They propose a dynamical low-rank approximation that relies more on the intrinsic dimensionality of the manifold the solution lives on. This method is designed to follow the flow of the dynamics and to identify the important features in the evolution. Since the equation is projected onto a solution manifold of lower dimensionality, only the core information is preserved and some redundant information is thrown away. The numerical cost thus depends mainly on the intrinsic dimensionality of the dynamics, rather than the total number of discrete points. The method was first utilized to deal with coupled ODE systems [36] and it has been demonstrated that, for a range of problems, it can preserve the accuracy of the solution.

The success in dealing with ODE systems inspired the generalization to PDEs. More recently, a class of efficient numerical schemes were proposed for kinetic equations. In the past few years, the method was applied to the Vlasov–Poisson system [28, 15], the Vlasov-Maxwell system  [16] and the classical Boltzmann equation [11] in both collisionless and strong-collisional regimes, and it is observed that both the sophisticated Landau damping phenomenon for the VP system [45, 41], and the incompressible Navier-Stokes limit for the Boltzmann equation are captured rather accurately with low cost. Different equations may require slightly modified dynamical low-rank algorithms to suit the specific structures of the equation, but the general approach is quite similar: one looks for the main features in the evolution and follows the flow of the equation projected onto the solution manifold of a certain low rank.

However, despite the strong intuition and the many promising numerical experiments, a mathematical error analysis is largely absent from all of those works (especially in the PDE setting). There are two major difficulties here. First, we often do not have results to show that the true solution is indeed approximately of low-rank. This point is more subtle than one might think. In fact, a Fourier expansion is also a low-rank approximation. But, a Fourier expansion uses a fixed set of basis functions, leading to unnecessary high ranks. In addition, showing low-rank by performing a Fourier decomposition relies on assuming smoothness of the solution, which for kinetic equations is certainly problematic. It should be noted, however, that for elliptic equations some results can be obtained [6].

The second difficulty is that we have to show that the numerical method actually captures the low-rank structure of the solution. For the ODE case an analysis has been performed in [26] and related works, but this analysis only applies to the non-stiff case and is thus not applicable to PDEs. The only mathematical analysis of a dynamical low-rank scheme in the PDE setting we are aware of is found in [43]. This paper considers a special situation where the stiffness originates only from the linear differential operator whose corresponding flow is exactly computable within the low-rank formulation. This separation suggests a splitting procedure that decouples the stiff and non-stiff dynamics, which permits the convergence of the resulting method. This is not the situation for kinetic equations, and thus the analysis does not apply.

The goal of the current paper is therefore to fill this void. In particular, we will consider a multi-scale radiative transfer equation, essentially a linearized Boltzmann equation. For this equation we will propose a dynamical low-rank algorithm and theoretically justify the validity and efficiency of this method in the strong-collisional regime (small Knudsen number). More specifically, we will show that in the small Knudsen number limit, the method is asymptotic preserving (that is, it correctly approximates the diffusion limit of the radiative transfer equation) and that it automatically captures the corresponding low-rank structure of the solution. It implies that the required degrees of freedom only scale as the number of grid points in physical space, in contrast to the full phase space as is required by other methods. This is, to our knowledge, the first numerical analysis result of this kind. Furthermore, our hope is that the present work will serve as a stepping stone for future studies for more complicated kinetic equations.

1.1. Multi-scale radiative transfer equation

We first give a quick overview of the equation. In dimension, the equation writes:


where is the scattering cross-section, , is called the Knudsen number, which represents the ratio between the mean free path and the typical domain length, and is the angular variable. The density is defined as follows


where is the area of , the unit sphere in dimensional space, and denotes the average with respect to . The equation is written in diffusion scaling, meaning the transport term and the collision term are multiplied by and , respectively.

Under certain regularity assumption on , in the limit of , the solution of (1) will asymptotically approach and thus loses its dependence. Then the solution of satisfies the following heat equation (with being a constant depending on ):


This is called the diffusion limit of the radiative transfer equation. This limit is particularly interesting in the present setting as the corresponding solution , indicating the solution is essentially of rank .

There are two major problems that pose difficulties for numerically solving equation (1). First, a small Knudsen number , the regime of interest in this paper, implies that both the transport term and the collision operator are stiff. A standard numerical scheme thus requires time steps that scale with and a very fine mesh. Second, equation (1) is posed in an dimensional phase space. Therefore, if the equation is discretized using points in each spatial direction and points in each velocity direction, we need to store at least floating point numbers. For this is a five-dimensional problem. The large increase in the degrees of freedom for such high-dimensional problems is usually referred to as the curse of dimensionality in the literature. Both of these issues contribute to the fact that solving equation (1) is extremely expensive from a computational point of view.

The second problem has been essentially left open. It is a rather widely-accepted fact that the numerical cost depends on the number of grid points (the degrees of freedom) used in the whole domain. However, there has been a large body of work addressing the first problem; namely, can we design a numerical solver that advances in time with a time-step size decoupled from the stiffness of the equation. A typical solution is to include some form of implicit solver in the numerical treatment in order to enlarge the stability region, and then devise a numerical method that efficiently solves the implicit part of the scheme. This property was later termed asymptotic-preserving (AP) in a ground-breaking paper [21], although some earlier schemes exists that were designed to satisfy this property [29]. A systematic investigation of AP was then conducted in a series of works for many (mainly nonlinear) Boltzmann-like equations, see [8, 30, 3, 9, 31, 7]. Also see reviews [20, 23, 10].

1.2. Dynamical low rank approximation

Dynamical low rank approximation is a systematic approach to tackle the curse of dimensionality for time-dependent problems. Under the assumption that the solution in fact lives in a low-dimensional manifold, the method looks for the low-rank approximation to the solution at every time step and advances the evolution projected on a low-dimensional tangential space. In particular, for kinetic equations, it is rather straightforward to separate the physical space coordinate and the velocity space:


so that the low rank factors and depend only on one type of coordinate. Assuming , this approximation only requires degrees of freedom.

Historically, dynamical low-rank approximations have been considered extensively in quantum mechanics, in which the dimensionality is high. Finding a low rank approximation there makes the computation tractable, see [40, 39] and [33, 34, 4] for a mathematical treatment. The application of the method in a general setting is studied in [26, 27, 37, 1]

, where the authors study both the matrix case (the approach we will be using here), and a general tensor formats, as is required for very high dimensional problems from quantum mechanics. One significant disadvantage of these algorithm is that they are not robust with respect to over-approximation, i.e. choosing a rank larger than is necessary for the problem at hand. Initially this problem was solved by regularization, and then a major improvement was suggested in

[36], in which the so-called projector splitting integrator

was introduced to make the dynamical low-rank approximation robust with respect to small singular values

[25]. This approach was later extended to various tensor formats [34, 35, 19, 38], and it is the approach to be utilized in the present paper.

The application of dynamical low-rank approximation for kinetic equations is relatively recent. Specifically, numerical methods have been proposed for the Vlasov–Poisson [15, 14], the Vlasov–Maxwell [16], and the Boltzmann equation [11]. These are all nonlinear kinetic models and the numerical results are very promising. However, the currently available results all mainly focus on computational aspects and on conserving the physical structure of the underlying equation. A mathematical analysis is lacking.

To analyze the numerical error for a dynamical low-rank algorithm applied to kinetic equations, we have to answer two questions

  1. Do the solutions have low-rank structures?

  2. How can we capture the structure dynamically in numerics?

Intuitively, to answer the first question, we utilize the fluid limit obtained on the theoretical level going back to equations (1) and (3). Since one can show that as we get , the velocity direction completely degenerates, and the rank of the representation in equation (4) is simply . It is then reasonable to expect that for small the solution is only slightly different from its rank- approximation.

To answer the second question, however, requires us to design an appropriate algorithm. As mentioned above, we will use the projector-splitting approach. This leads to a set of three evolution equations for , , and , respectively. To complete one time step we have to advance all three equations. As will be shown in Section 3, however, the order in which the sub-flows in the splitting are solved is important for preserving the rank. In addition, we will demonstrate that the chosen time integrator plays a crucial role. One can, for example, show that the implicit Euler method, due to the lack of the symmetry, fails to satisfy the correct limit. On the other hand, the Crank–Nicolson method converges to the correct limit and is thus asymptotic preserving.

Before proceeding, let us put the present work in the proper context. There are two basic methods to compute a low-rank approximation of an evolutionary partial differential equations. The first option is to use a so called

projection method (not to be confused with the projector-splitting integrator for a dynamical low-rank approximation). In this case we first discretize the partial differential equations. We then start with an initial value of fixed rank and compute one time step using the numerical integrator chosen. This, in general, takes us outside of the approximation space (i.e. the rank increases). The so obtained result is then projected back to a manifold of functions with fixed rank. For kinetic equations such an approach was suggested in [28]. The disadvantage of this approach, however, is that an object that does not lie in the approximation space has to be constructed, which implies both a memory as well as a performance penalty. In addition, the derived algorithm is very closely tied to the specific space and time discretization that has been chosen. An alternative is to directly formulate the dynamics of the low-rank approximation by projecting the original equation onto the manifold of functions with a fixed rank. This can be done completely on the continuous level and results in a new set of partial differential equations formulated directly in terms the low-rank factors. As mentioned above, we will follow the latter approach which is usually referred to as a dynamical low-rank approximation.

The rest of the paper is organized as follows. In Section 2 we present the numerical method that is applied to the radiative transfer equation. In Section 3 we state and prove the main results. Some parts of the proof in Section 3 are rather tedious. We leave those to the appendix and keep only the core analysis in the main text. Numerical results are then presented in section 4.

2. Numerical scheme

In this section we follow the framework and notations in [15] and present a numerical method that uses a function of the form given in equation (4) as the approximation space. Our goal is to obtain a low rank approximation to the solution of  (1). To define the rank, we first need to equip the spaces with proper measures. For that, we simply use the standard vector space in both the spatial and the velocity domain. That is, we have

where is the spatial measure in , and is a normalized measure in . Then, any function that can be expanded by a set of orthonormal basis functions in and , is called rank-. We collect all these functions together and denote the collection by .

Definition 1 (Rank- function in ).

The collection of all rank- functions is denoted by

where we call a function rank- if there is a set of orthonormal basis functions and a set of orthonormal basis functions such that

Here are orthonormal in physical space and are orthonormal in velocity space with appropriate measures and , namely

We note that in this definition, only the rank, , is fixed. The basis functions and can be arbitrary, as long as the orthogonality condition is satisfied. We further emphasize that is not a function space; it is easily seen that the summation of two rank- functions may not be rank-.

It is unlikely that the solution is of rank-, i.e. in , for all time. However, numerically one can argue that the solution is approximately of low rank. Thus for the numerical solution we seek a rank- approximation in at every time step. The algorithm is consequently looking for a trajectory on the manifold that resembles the evolution guided by the equation. In some sense, we need to project the equation in to the manifold and find the equation that governs the dynamics of this trajectory in . Let us denote by the analytic solution and by the numerical rank- approximation. Using the argument above, the governing equation for is


where stands for the projection of onto the tangential plane of at . This is to ensure that

where denotes the tangential plane of at . This condition is sufficient to guarantee that lies in the manifold for all time.

To explicitly express the tangential plane for the manifold and the projection operator, we first notice that at any function , the tangential plane is different and that the form highly depends on . Denote

then the tangential plane is given by

This set collects all functions whose infinitesimal, when added to , still yields a rank- function. In this definition, we notice that we are allowed to choose arbitrarily an matrix , a function list and a function list , as long as the gauge conditions, , are satisfied. We note that the gauge conditions are imposed to guarantee the uniqueness of the low-rank factors. Interested readers are referred to [36] for details.

With this definition, one has:

  • If , the orthogonality and gauge condition quickly allows us to write

  • if , its projection onto can also be easily formulated as follows


    where the spatial and velocity projection are

Inserting (6) into (5) gives us the governing equation for :


The numerical method will then be developed upon this formulation. We discuss the semi-discrete (in time) and the fully-discrete schemes in details in the following subsections. This will give us a formulation for one time step. To devise an asymptotic preserving scheme for integrating from to the final time we need a thorough understanding of the error analysis. In fact, the implicit Euler and Crank-Nicolson behave slightly different. To optimally suit our purpose we run the Euler for one step before shifting to Crank-Nicolson scheme. We defer a detailed discussion to the end of Section 3.

Remark 1.

We wrap up this introduction with a comment on rank- approximations. In fact, all numerical methods to solve PDEs are rank- approximations. Let us suppose we have a numerical method with grid points in both and . Then

where are basis functions we use to approximate the solution. For example, in a finite difference method we set as the Kronecker delta function or hat functions that peak at (similar argument holds for ) and is the numerical solution evaluated at and . To find a numerical solution on this mesh is equivalent to finding a rank- approximation with fixed and . Then plays the role of above. This is typically not formulated in this way as the rank is very large. When we consider a low-rank approximation, as in this paper, we will always assume that . To achieve this, the functions and need to evolve in time according to the dynamics of the equation.

2.1. Semi-discrete low rank splitting method

In this section, we develop a projector-splitting method to solve equation (7). In order to be concise, we simply denote the low rank solution by and we decompose the solution using its low-rank representation:


where and collects the basis functions

This formulation will also be quite useful later when we introduce a space discretization. We will also be using quantities and :


Computing (7) at discrete times amounts to finding the governing equations that provide the updates of

respectively, for all with .

To specify the initial data into the format, we project the initial condition onto

using the singular value decomposition (SVD):


For updating , the Lie-Trotter splitting is used. From time step to , we split the three operators on the right hand side of (7) into three sub-steps, namely


This splitting (12)-(14) takes place for each time step. That is, given the numerical solution at time we update the solution to by solving the three equations one after another. All equations are advanced for a full time step . Different from directly solving (7), each sub-step only changes one part of the decomposition. In particular, the first splitting step (12) preserves , the last (14) preserves , and the middle step (13) merely updates . This allows us to update the three components separately without disturbing others. Below we detail the evolution of each sub-step:

  1. [-,topsep=0pt]

  2. Updating (12):

    Starting with , in this sub-step, we run (12) for a full time step , and we denote the result by . Since the step preserves :

    To update and , we plug the low rank formulation (8) into (12) and obtain


    with . By using equation (9) we can simplify (15) to


    where , and and , both , are matrix versions of the differential operator and the scattering operator respectively:


    We denote the solution to (16) by , and and

    are obtained through the Gram-Schmidt process (QR factorization) which ensures the orthogonality of


  3. Updating (13):

    In this step, equation (13) is ran for a full time step with initial condition . We denote the result by . Since the step only changes , we immediately obtain

    Plugging the low rank representation (8) into (13), we obtain, for ,

    In matrix form this can be written as


    where and , both , are the matrix versions of the multiplication operator associated to () and the density term, respectively:


    We denote the computed update of (18) by .

  4. Updating (14):

    In this step, (14) is ran for a full time step with initial value . We denote the result by . This step preserves . Thus,

    To update and , the low rank formulation (8) is plugging into (12). Then for all we have

    This can be written in matrix form as follows


    with (since )

    Solving this equation we obtain , and the orthogonality of is ensured through the Gram-Schmidt process:

With these three steps completed, one finally arrives at numerical solution at time

2.2. Fully-discrete low rank splitting method

The beforementioned method will now be discretized. Due to the stiffness of the equations, an implicit scheme will be applied. We denote

the sets of discrete points in and . The discrete solution can then be represented as follows

where and


with and denote the -th mode evaluated at the discrete points in and .

We now discuss the implicit time integration of the three equations (16), (18), and (20

). The asymptotic analysis will be performed based on this fully-discrete formulation.

  • For updating (16), is preserved and thus . The direct application of the implicit Euler scheme gives


    and if Crank–Nicolson is used, the scheme is written as


    We have used


    with . is an matrix with evaluations of at the grid points in assigned as diagonal entries and is the discrete approximation of . The specific form of depends on the spatial discretization. We have used the simple rectangle rule for the integration in . This determines the form of . Other numerical integral rules could also be applied and the specific form of will change accordingly. Obviously and are the discrete versions of (17).

  • For updating (18), we note that and are preserved in this step. That is,

    The direct application of the implicit Euler scheme gives


    Similarly, if Crank–Nicolson is used, the scheme can be written as


    Here and , both , are discrete versions of (19), defined by

  • Finally for updating (16), we note , , are preserved. Defining


    and applying the implicit Euler method we obtain


    where we used the same definition of . Crank–Nicolson method will not be used in the last step and thus we do not specify it.

We finalize the update for

3. Properties of the numerical scheme

We investigate the properties of the numerical method in this section. In particular we will discuss the computational complexity and prove the method preserves the asymptotic limit in the strong collisional regime.

3.1. Computational complexity

To analyze the computational complexity is rather straightforward. Denote the rank, and the number of grid points in physical space and velocity space respectively. The matrices in the updating formula, , , are computed by solving (22), (25), (29). The following cost incurs:

  • Preparation: Calculation of , needs floating point operations (flops).

  • Upadate : Calculation of needs flops.

    Calculation of , need flops.

    Solving needs flops.

    Using QR decomposition to obtain , needs .

  • Update : Calculation of and needs

    Calculation of , needs flops.

    Solving needs flops.

  • Update : Calculation of needs needs flops.

    Calculation of , needs flops.

    Solving needs flops.

    Using QR decomposition to obtain , needs .

Because , in conclusion, we need flops per time step, instead of the .

3.2. Intuition of the error analysis

Before describing and proving our results in detail, in this section we first give a relatively vague justification on why the method works. As described in the introduction, there are two points we need to make:

  • Why the true solution has an approximate low rank structure? This question is a rather fundamental, and is independent of the method chosen: the rank structure of the solution purely depends on the governing PDE we are studying here.

  • Why the method keeps track of the low rank structure? This question concerns the behavior of the specific method (dynamic low-rank approximation) we choose to use here.

The two questions will be addressed respectively in the following two subsections.

3.2.1. Foundation: the low rank structure in the solution

To justify the low rank structure of the linearized Boltzmann equation, we simply will cite the following result [2]

Theorem 1 (Theorem 2 of [2]).

Denote the solution to the linear Boltzmann equation


where and . Then in the zero limit of , the solution converges to the solution of the diffusion equation:

in the sense that

We note that periodic boundary conditions are imposed in the original theorem to avoid complications that may come from the boundary layers. Essentially this theorem states that in the zero limit of , loses its velocity dependence, and the dynamics will purely be reflected in the physical space. In some sense:

can be seen as the rank- approximation with the threshold set at any value bigger than .

The proof for the theorem follows the Hilbert expansion. Formally, writing , we plug it back in the original equation and get:


in which the last equation gives us the diffusion limit. As seen in the expansion, the essence of the proof mainly lies in showing that

and then tracing the dynamics of in the null space. Showing Theorem 1 rigorously then amounts to bounding uniformly in , and we omit it from here.

3.2.2. Capture the structure along the dynamics

It is not straightforward to show that the dynamical low-rank approximation method captures the rank structure of the solution. In fact, the direct calculation would suggest otherwise: recall the dynamic low-rank approximation solution governed by the following equation


and compare it with the original one (30), it is rather easy to see that one essentially need to show that

in some norm, where and are basis functions of . This does not seem to be an easy task. On the contrary, concerning that contains stiff terms such as , a brute-force calculation would suggest that the error is of order .

This straightforward approach, however, overlooks the information hidden in the equation. To show the method keeps small along the evolution, some delicacy from the equation needs to be employed. As will be presented in further details in Section 3.3, the proof largely relies on the separation of scales and one needs to perform the order by order matching of the scales to derive a clearer view of the solution structure. For that we perform the same analysis as done in (3.2.1), and asymptotically expand . Plugging it back into equation (32), one has, in the leading order:

Noting that from the definition of ,


which further suggests that the leading order of the numerical solution

This at least implies the reduced order equation (32) has its leading order lying in the correct space. Whether it preserves the correct equilibrium state is up to more delicate derivations in the higher order expansions, and thus is left to the following section, where we directly tackle the problem in the fully discrete setting, and show that the numerical scheme indeed captures the diffusion phenomenon.

3.3. Asymptotic preserving

Following the intuition from the previous section, we give the rigorous proof that shows that the method captures the diffusion limit in the limit.

It is a rather challenging task to design a numerical method for a stiff equation with the discretization being independent of the smallest scale. For an equation with small dependence, usually the solution shows variations at fine scale, and in order to preserve these variations, the discretization has to be small, leading to a large number of degrees of freedom, driving up the numerical cost. If a method has its discretization relaxed from requirements at the finest scale, while still preserves the asymptotic limit of the equation, we call it an asymptotic preserving (AP) method. It is an attractive property for a numerical method to have.

In this section, we will prove that in our method, by injecting the low-rank structure into numerical solutions, one can automatically capture the diffusion limit (3), with , , and all independent of , and thus is an AP method. Furthermore, can be as small as . We note that the analysis for Crank-Nicolson and implicit Euler shares some similarities but Crank-Nicolson enjoys a symmetry that allows us to pass to the asymptotic limit.

3.3.1. Asymptotic analysis of the implicit Euler method

To unify the notation, throughout this section, we denote a matrix of size with being the numerical solution at . Denote a column vector of length , then is a vector of length representing the discrete version of at time for spatial grid points . The hope is to show that density solves equation (3) in the limit .

Our first result concerns the behavior of the implicit Euler method in the limit .

Theorem 2.

We employ the implicit Euler method to compute from , using equations (22), (25) and (