Solving the 1D diffusion equation (a PDE) using Stan

[Ben Lambert]

Hello all,

I have been playing around with Stan's ``integrate_ode'' functionality - really nice by the way! I was wondering - how easy would it be to do inference for a PDE model in Stan?

As an example, I have been trying to infer the diffusion parameter for a 1D diffusion equation: 

du/dt = D d^2 u / dx^2

where u=u(x,t) measures a population density of some species, and D is the diffusion coefficient. Further, suppose that we know the initial conditions: u(0,t) = f(x), as well as (Dirichlet) boundary conditions: u(x,0) = 0, and u(x,L) = 0. Finally if we assume a normal sampling model about the mean u(x,t), then we have an inference problem.

Typically we would solve this equation by some sort of finite difference method, where we discretise in both x and t. However, if we have a decent integrator, then we should be able to just discretise in distance, and then use the time integrator for each grid point in x. If we use a centralised finite difference approximation for x, then this amounts to integrating:

d v(i,t) / dt = D [v(i+1,t) - 2 v(i,t) + v(i-1,t)] / 2h

where v(i,t) is the approximate solution at grid point i at time t, and h is the step size.

Ideally I was hoping we could use Stan's integrate_ode function for each of the grid points i, to calculate the mean v(i,t), and then use that to increment the log probability by: ~ normal(v(i,t), sigma). At the moment however, I am a little unclear as to how to do this.

My thoughts were that we could essentially do the following at each time step:

vector[T] yHat;

for(i in 1:N){

    real y_hat_individual[T,1];

    y_hat_individual <- integrate_ode()

}

However, I'm not sure how to write the derivative function (that takes three values of v, and returns a single derivative as per the above finite difference approximation.) 

Does anyone have an idea as to how I would do this?

Best,

Ben

[Michael Betancourt]

ODEs, while by no means trivial, are relatively straightforward

and hence admit general solvers like CVODE that we can plug

into the Stan Math Library.  PDEs, however, are much, much

more complicated and don’t admit such generic solvers, hence

we cannot hope to have anything like a generic “integrate_pde”

solver.  Unfortunately PDEs require expert-motivated, bespoke

solvers and the methods are evolving so quickly that there aren’t

the generic libraries we would need to easily be able to integrate

even these bespoke methods into Stan.

That means you’ll have to make approximations and argue their

validity on a case by case basis.  Discretizing over space is one

such approximation which can be okay in some cases and 

disastrous in others.  You can also discretize over space and

time, giving you something equivalent to a (hopefully sparse) 

spatial model — see any of the INLA doc for much more on this 

approach.

--
You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
To post to this group, send email to stan-...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Ben Lambert]

Hi Mike,

For clarity - I wasn't advocating building a PDE solver in Stan. I deal with PDEs all the time, and I know how the concept of a one-size-fits-all solution is ridiculous. My apologies if my question didn't make this clear.

What I was asking was - for the problem at hand - a 1D diffusion equation with measurement error, whether it would be possible to implement such a solver in Stan. As I said, I could use some sort of finite difference (in distance) scheme for this, then use Stan's integrator to integrate over time. For the problem at hand - a simple elliptic PDE with Dirichlet boundary conditions - I was wondering how to do this.

Best,

Ben

[Michael Betancourt]

What I was asking was - for the problem at hand - a 1D diffusion equation with measurement error, whether it would be possible to implement such a solver in Stan. As I said, I could use some sort of finite difference (in distance) scheme for this, then use Stan's integrator to integrate over time. For the problem at hand - a simple elliptic PDE with Dirichlet boundary conditions - I was wondering how to do this.

I know, and I second my initial response.  A good discretization, both in

terms of design and implementation, requires expert knowledge of the 

system at hand.  Once the system of PDEs has been reduced to an

approximate system of ODEs that approximate system can be written

up and fit in Stan, but until then there’s nothing Stan can do.  Again,

see the INLA literature for example discretizations and their

implementations, including ways of dealing with various derivatives.