# Stan demo: stochastic approach of PDE solutions

This =Stan= demo is a translation of one of my previosu R projects, implementing the Walk-on-sphere(WOS) method to solve the Laplace equation on a rectangle. This is the so-called probabilistic mesh-free method for PDE solution. Essentially it’s based on the probabilistic interpretations of Laplace equation, which says the solution is the mean of exit-points’ function values of Brownian motion.

Solution of the 2D Laplace equation by the src below on a rectangle domain. More details here

Source code

``````/* This Stan demo of probalistic appraoch of PDE solution.
It uses Walk-on-sphere method to calculate Laplace equ
solution on a rectangle.
*/

/* run with */
/* .harmonic sample num_samples=1 algorithm=fixed_param data file=./harmonic.data.R */
functions {
real[] rect_boundary_x() {
real xb = { 0.0, 1.0 };  /* rectangle boundary */
return xb;
}
real[] rect_boundary_y() {
real yb = { 0.0, 1.0 };  /* rectangle boundary */
return yb;
}
real bc(real x, real y) {
real xb = rect_boundary_x();
real yb = rect_boundary_y();
real val;
if (x == xb) {           /* left boundary */
val = 0.0;
} else if ( x == xb) {   /* right boundary */
val = 0.0;
} else if ( y == yb) {   /* lower boundary */
val = 0.0;
} else if ( y == yb) {   /* upper boundary */
if (x <= 2.0/3.0)
val = 75*x;
else
val = 150 * (1-x);
}
return val;
}

real rectangle_wos_rng(real x, real y, real tol) {
real xb = rect_boundary_x();
real yb = rect_boundary_y();
real res = {x, y, 0.0};
real dist = { xb - res, res - xb, yb - res, res - yb};
real r = min(dist);
real val;
while (r > tol) {
real theta = uniform_rng(0, 2*pi());
res = res + r * cos(theta);
res = res + r * sin(theta);
dist = { xb - res, res - xb, yb - res, res - yb};
r = min(dist);
res = r;
}
if (dist < tol ) {       /* right boundary */
res = xb;
} else if (dist < tol ) { /* left boundary */
res = xb;
} else if (dist < tol ) { /* upper boundary */
res = yb;
} else if (dist < tol ) { /* lower boundary */
res = yb;
}
val = bc(res, res);
res = val;
return val;
}
}

data{
real tolerance;
int m;
int n;
int N;
}

transformed data {
vector[N] bcsample;
real x;
real y;
real hm = 1.0/m;
real hn = 1.0/n;
real sol;
for ( i in 1:m+1 ) {
for ( j in 1:n+1 ) {
x = (i-1)*hm;
y = (j-1)*hn;
for ( k in 1:N ) {
bcsample[k] = rectangle_wos_rng(x, y, tolerance);
}
sol = mean(bcsample);
}
}
}
``````
4 Likes

Why put this in Stan given that it’s just computing a function using transformed data? Is there a plan to then embed it in a larger application? Or to use it for parameters and model inference?

I’d like to use this model to test two things. One, using multiple MPI communicators for embarrassingly parallel evaluations. Two, testing inference on PDE boundary conditions.