Struct Milstein

pub struct Milstein<T: Real, V: State<T>, D: CallBackData> {
    pub h: T,
    /* private fields */
}

Milstein Method for solving stochastic differential equations.

The Milstein method is a higher-order method for SDEs that includes an additional term from the Itô-Taylor expansion to achieve better accuracy.

For an SDE of the form: dY = a(t,Y)dt + b(t,Y)dW

The Milstein update is: Y_{n+1} = Y_n + a(t_n, Y_n)Δt + b(t_n, Y_n)ΔW_n + 0.5 * b(t_n, Y_n)^2 * [(ΔW_n)² - Δt]

where ΔW_n is a Wiener process increment.

This implementation is based on the standard form of the Milstein method for scalar SDEs. For geometric Brownian motion, the b(t,Y) term is σY, so the correction term becomes 0.5 * (σY)^2 * [(ΔW)² - Δt].

The Milstein method has strong order of convergence 1.0 (compared to 0.5 for Euler-Maruyama), making it more accurate for SDEs with significant diffusion effects.

Example

use differential_equations::prelude::*;
use nalgebra::SVector;
use rand::SeedableRng;
use rand_distr::{Distribution, Normal};

struct GBM {
    mu: f64,     // Drift rate
    sigma: f64,  // Volatility
    rng: rand::rngs::StdRng,
}

impl GBM {
    fn new(mu: f64, sigma: f64, seed: u64) -> Self {
        Self {
            mu,
            sigma,
            rng: rand::rngs::StdRng::seed_from_u64(seed),
        }
    }
}

impl SDE<f64, SVector<f64, 1>> for GBM {
    fn drift(&self, _t: f64, y: &SVector<f64, 1>, dydt: &mut SVector<f64, 1>) {
        dydt[0] = self.mu * y[0];
    }
     
    fn diffusion(&self, _t: f64, y: &SVector<f64, 1>, dydw: &mut SVector<f64, 1>) {
        dydw[0] = self.sigma * y[0];
    }
     
    fn noise(&self, dt: f64, dw: &mut SVector<f64, 1>) {
        let normal = Normal::new(0.0, dt.sqrt()).unwrap();
        dw[0] = normal.sample(&mut self.rng.clone());
    }
}

let t0 = 0.0;
let tf = 1.0;
let y0 = SVector::<f64, 1>::new(100.0);
let mut solver = Milstein::new(0.01);
let gbm = GBM::new(0.1, 0.2, 42);
let gbm_problem = SDEProblem::new(gbm, t0, tf, y0);

// Solve the SDE
let result = gbm_problem.solve(&mut solver);

Fields

h: T

Implementations

impl<T: Real, V: State<T>, D: CallBackData> Milstein<T, V, D>

pub fn new(h: T) -> Self

Create a new Milstein solver with the specified step size

Arguments

• h - Step size

Returns

• A new solver instance

Trait Implementations

impl<T: Real, V: State<T>, D: CallBackData> Default for Milstein<T, V, D>

fn default() -> Self

Returns the “default value” for a type.

impl<T: Real, V: State<T>, D: CallBackData> Interpolation<T, V> for Milstein<T, V, D>

fn interpolate(&mut self, t_interp: T) -> Result<V, Error<T, V>>

Interpolate between previous and current step

impl<T: Real, V: State<T>, D: CallBackData> SDENumericalMethod<T, V, D> for Milstein<T, V, D>

fn init<F>( &mut self, sde: &F, t0: T, tf: T, y0: &V, ) -> Result<Evals, Error<T, V>>

where F: SDE<T, V, D>,

Initialize SDENumericalMethod before solving SDE

fn step<F>(&mut self, sde: &F) -> Result<Evals, Error<T, V>>

where F: SDE<T, V, D>,

Step through solving the SDE by one step

fn t(&self) -> T

Access time of last accepted step

fn y(&self) -> &V

Access solution of last accepted step

fn t_prev(&self) -> T

Access time of previous accepted step

fn y_prev(&self) -> &V

Access solution of previous accepted step

fn h(&self) -> T

Access step size of next step

fn set_h(&mut self, h: T)

Set step size of next step

fn status(&self) -> &Status<T, V, D>

Status of solver

fn set_status(&mut self, status: Status<T, V, D>)

Set status of solver