Решение диффузионной конвекции SDE представление с использованием метода Milstein с Boost (C ++) - не сходится к детерми ⇐ C++

[Anonymous]

Я работаю над исследовательским проектом по моделированию уравнения диффузионного адвекции, как здесь, предполагая нулевую дивергенцию. Я решил PDE, используя RK4 и конечные различия и сравниваю результат с решением стохастического представления от Feynman Kac ($ dx = v dt+d dw $), используя метод Милштейна (слайды для стохастической формы Слайд 64).
Но когда я построю обе обе версии, стохастическое решение не соответствует детерминированному. Ниже красным я сделал гистограмму распределения частиц, которые я смоделировал, и построил центр гистограммы. Они не соответствуют решению синей функции.
< Br /> ниже мой код C ++ и мой сценарий Python, который создал сюжет. В качестве примера я использую пользовательский шаг Boost (пользовательский шаговый документ). sin^2 (x), v_0 = 5.
#include // for std::transform and std::accumulate
#include
#include
#include
#include
#include
#include
#include

using namespace std;
const static size_t N = 1;
typedef vector ParticlePositions;
const static double v = 0.1; // Convection coefficient
const static double D = 0.009; // Diffusion coefficient
const static double f = 0.0; // Source Term
typedef boost::array< double , N > state_type;

// Function to calculate initial Gaussian distribution
vector initialize_gaussian(double amp, const vector& x)
{
vector u(x.size());
double sigma = 0.01;
for (size_t i = 0; i < x.size(); ++i) {
u = amp * (1 / sqrt(0.001 * 2 * M_PI)) * exp(-0.5 * (x * x / 0.001));
}
return u;
}

double D_eval(double x) {
return D + (D * x * x);
}
double v_eval(double x) {
return 50 * v * sin(x) * sin(x);
}

double dD_dx(double x) {
return 2 * D * x;
}

// ========================
// === Stochastic Class ===
// === dx=v*dt+D*dW ===
// ========================

// Euler Maruyama: xn+1=xn+dt(w+alpha*dW)
template
class stochastic_euler {
public:
typedef boost::array state_type;
typedef boost::array deriv_type;
typedef double value_type;
typedef double time_type;
typedef unsigned short order_type;
typedef boost::numeric::odeint::stepper_tag stepper_category;

static order_type order(void) { return 1; }

template
void do_step(System system, state_type& x, time_type t, time_type dt) const
{
deriv_type det, stoch;
system.first(x, det);
system.second(x, stoch);
for (size_t i = 0; i < x.size(); ++i) {
double milCorrect = 4 * D_eval(x) * abs(dD_dx(x));
x += dt * (det - milCorrect) + \
+ sqrt(milCorrect * dt * dt) * stoch * stoch + \
sqrt(2 * D_eval(x) * dt) * stoch[i]; // sqrt(dt) is standard dev
}
}
};

// Pair of functions for the deterministic and stochastic parts of the SDE
struct sde_det {
void operator()(const state_type& x, state_type& dxdt) const
{
for (size_t i = 0; i < x.size(); ++i)
dxdt[i] = -2*v_eval(x[i]);
}
};

struct sde_stoch {
mt19937& m_rng;
normal_distribution m_dist;
// Constructor
sde_stoch(mt19937& rng, double sigma) : m_rng(rng), m_dist(0.0, sigma) {}

void operator()(const state_type& x, state_type& dxdt)
{
for (size_t i = 0; i < x.size(); ++i)
dxdt[i] = m_dist(m_rng);
}
};

// Used to store particle positions for each time step into vector argument
struct particle_observer {
ParticlePositions& positions;
int particle_index;
int step;

// Constructor
particle_observer(ParticlePositions& positions, int particle_index)
: positions(positions), particle_index(particle_index), step(0) {}

void operator()(const state_type& x, double t)
{
positions[particle_index][step] = x[0];
step++;
}
};

// Compute PDE derivative at index i for the current solution u_curr
double compute_rhs(const vector& u_curr, int i, double dx) {
double x_i = (i * dx) - 2.0; // physical position
double Dval = D_eval(x_i);
double dDval = dD_dx(x_i);
double w = v_eval(x_i);

double diff_term = (u_curr[i+1] - 2.0 * u_curr[i] + u_curr[i-1]) / (dx * dx);
double convection_term = (u_curr[i+1] - u_curr[i-1]) / (2.0 * dx);

// du/dt
return Dval * diff_term + dDval * convection_term + w * convection_term;
}

// Perform one RK4 time step, updating sol[j] from sol[j-1]
void update_solution_rk4(vector& sol, double dx, double dt, int nx, int j) {
const vector& u_old = sol[j-1];
vector k1(nx, 0.0), k2(nx, 0.0), k3(nx, 0.0), k4(nx, 0.0);
vector u_temp1(u_old), u_temp2(u_old), u_temp3(u_old);
vector u_new(u_old);

// Compute k1 for interior points
for(int i = 1; i < nx - 1; ++i) {
k1[i] = compute_rhs(u_old, i, dx);
u_temp1[i] = u_old[i] + 0.5 * dt * k1[i];
k2[i] = compute_rhs(u_temp1, i, dx);
u_temp2[i] = u_old[i] + 0.5 * dt * k2[i];
k3[i] = compute_rhs(u_temp2, i, dx);
u_temp3[i] = u_old[i] + dt * k3[i];
k4[i] = compute_rhs(u_temp3, i, dx);
u_new[i] = u_old[i] + (dt / 6.0) * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]);
}

// Apply or maintain boundary conditions if needed (e.g., Dirichlet = 0)
u_new[0] = 0; //u_old[0];
u_new[nx-1] = 0; //u_old[nx-1];

// Store result in sol[j]
sol[j] = u_new;
}

int main()
{
cout 0.5) {
cerr

Подробнее здесь: https://stackoverflow.com/questions/794 ... with-boost