Homogenous linear RODE with the sine of a Geometric Brownian motion coefficient

This time we take, as the coefficient of a homogeneous linear equation, the sine of a Geometric Brownian motion process. This is a multiplicative noise.

The equation

We consider the RODE

\[ \begin{cases} \displaystyle \frac{\mathrm{d}X_t}{\mathrm{d} t} = \sin(Y_t) X_t, \qquad 0 \leq t \leq T, \\ \left. X_t \right|_{t = 0} = X_0, \end{cases}\]

where $\{Y_t\}_{t\geq 0}$ is a geometric Brownian motion process. The explicit solution is

\[ X_t = e^{\int_0^t \sin(Y_s) \;\mathrm{d}s} X_0.\]

Computing a higher order approximation of the solution

As in the previous examples, the integral $\int_0^{t_j} \sin(Y_s)\;\mathrm{d}s$ and, hence, the exact solution, is not uniquely defined from the values $W_{t_j}$ of the noise on the mesh points. This time, an exact distribution for the collection of exact solutions conditioned on the mesh points is not available in closed form. Hence, we consider an approximation of an exact solution by solving the equation numerically, with the Euler method itself, but with a much higher mesh resolution.

Indeed, the convergence will be estimated from a set of discretizations with mesh points with time step $\Delta t_N = (t\_f - t\_0) / 2^N$, for $N = N_1 < N_2 < \ldots N_n$, for some $n\in \mathbb{N}$, by comparing the error of such solutions to an approximated solutions computed in a finer mesh with $\Delta t_{\textrm{fine}} = \Delta t_{N_n}^2$, hence with $N_\textrm{fine} = N_n^2$.

Numerical approximation

Setting up the problem

First we load the necessary packages

using Plots
using Random
using Distributions
using RODEConvergence

Then we set up some parameters, with a Distributions.Normal random variable as the initial condition.

rng = Xoshiro(123)

f(t, x, y, p) = sin(y) * x

params = nothing

t0, tf = 0.0, 1.0
x0law = Normal()

Distributions.Normal{Float64}(μ=0.0, σ=1.0)

The geometric Brownian motion noise is defined via GeometricBrownianMotionProcess, with initial value $y_0$, drift $\mu$, and dissipation $\sigma$ as given by

μ = 1.0
σ = 0.2
y0 = 1.0
noise = GeometricBrownianMotionProcess(t0, tf, y0, μ, σ)

GeometricBrownianMotionProcess{Float64}(0.0, 1.0, 1.0, 1.0, 0.2)

The mesh parameters are

ntgt = 2^18
ns = 2 .^ (4:9)

6-element Vector{Int64}:
  16
  32
  64
 128
 256
 512

and

nsample = ns[[1, 2, 3, 4]]

4-element Vector{Int64}:
  16
  32
  64
 128

The number of simulations for the Monte Carlo estimate is set to

m = 400

We add some information about the simulation, for the caption of the convergence figure.

info = (
    equation = "\$\\mathrm{d}X_t/\\mathrm{d}t = \\sin(Y_t) X_t\$",
    noise = "a geometric Brownian motion process noise \$\\{Y_t\\}_t\$ (ic=$y0, drift=$μ; diffusion=$σ)",
    ic = "\$X_0 \\sim \\mathcal{N}(0, 1)\$"
)

We define the target solution as the approximation obtained by the Euler method in the much higher resolution ntgt of mesh points. The approximations are also obtained via the Euler method, in the coarser meshes defined by ns.

target = RandomEuler()
method = RandomEuler()

RandomEuler{Float64, Distributions.Univariate}(Float64[])

Order of convergence

With all the parameters set up, we build the ConvergenceSuite:

suite = ConvergenceSuite(t0, tf, x0law, f, noise, params, target, method, ntgt, ns, m)

ConvergenceSuite{Float64, Distributions.Normal{Float64}, GeometricBrownianMotionProcess{Float64}, Nothing, typeof(Main.var"Main".f), 1, 1, RandomEuler{Float64, Distributions.Univariate}, RandomEuler{Float64, Distributions.Univariate}}(0.0, 1.0, Distributions.Normal{Float64}(μ=0.0, σ=1.0), Main.var"Main".f, GeometricBrownianMotionProcess{Float64}(0.0, 1.0, 1.0, 1.0, 0.2), nothing, RandomEuler{Float64, Distributions.Univariate}(Float64[]), RandomEuler{Float64, Distributions.Univariate}(Float64[]), 262144, [16, 32, 64, 128, 256, 512], 400, [1, 1, 1, 1, 1, 1], [-10.021430432325799, 1.3901747e-317, 6.9452245652387e-310, 2.547989973e-315, 0.0, 0.0, -6.907421002359041, -7.3081940879107865, -8.398285958468641, -8.606521037231587  …  -4.4610120006407685, -5.653656550000688, -5.918937121052622, -6.547402630423461, -7.3081940879107865, -8.398285958468641, -8.606521037231587, -9.249108707956637, -10.447822800086016, -10.758392253849077], [NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN  …  -8.398285958468641, -8.606521037231587, -9.249108707956637, -10.021430432325799, -10.758392253849077, -4.149887265017013, -4.909382477812359, -5.222434485176767, -5.918937121052622, -6.907421002359041], [-0.06389100157913419, -0.06389832123980099, -0.06399176896251398, -0.06401059154825825, -0.06410301059555304, -0.0642056926639871, -0.0644023558544875, -0.0646604258177009, -0.06500341799733947, -0.06541953538693149  …  -0.2998268101145127, -0.3002077746653318, -0.3006036347566172, -0.30100593902869144, -0.3013634499594975, -0.30165137532922903, -0.3019093966416785, -0.3023484223138704, -0.30271709411887887, -0.30289678603658343])

Then we are ready to compute the errors via solve:

@time result = solve(rng, suite)

  2.325944 seconds (807.65 k allocations: 39.917 MiB, 22.88% compilation time)

The computed strong error for each resolution in ns is stored in result.errors, and a raw LaTeX table can be displayed for inclusion in the article:

table = generate_error_table(result, suite, info)

    \begin{center}
        \begin{tabular}[htb]{|r|l|l|l|}
            \hline N & dt & error & std err \\
            \hline \hline
            16 & 0.0625 & 0.0285 & 0.00157 \\
            32 & 0.0312 & 0.0145 & 0.000784 \\
            64 & 0.0156 & 0.00734 & 0.000404 \\
            128 & 0.00781 & 0.00368 & 0.000201 \\
            256 & 0.00391 & 0.00186 & 0.0001 \\
            512 & 0.00195 & 0.000918 & 4.97e-5 \\
            \hline
        \end{tabular}
    \end{center}

    \bigskip

    \caption{Mesh points (N), time steps (dt), strong error (error), and standard error (std err) of the Euler method for $\mathrm{d}X_t/\mathrm{d}t = \sin(Y_t) X_t$ for each mesh resolution $N$, with initial condition $X_0 \sim \mathcal{N}(0, 1)$ and a geometric Brownian motion process noise $\{Y_t\}_t$ (ic=1.0, drift=1.0; diffusion=0.2), on the time interval $I = [0.0, 1.0]$, based on $M = 400$ sample paths for each fixed time step, with the target solution calculated with $262144$ points. The order of strong convergence is estimated to be $p = 0.99$, with the 95\% confidence interval $[0.8873, 1.0971]$.}

The calculated order of convergence is given by result.p:

println("Order of convergence `C Δtᵖ` with p = $(round(result.p, sigdigits=2)) and 95% confidence interval ($(round(result.pmin, sigdigits=3)), $(round(result.pmax, sigdigits=3)))")

Order of convergence `C Δtᵖ` with p = 0.99 and 95% confidence interval (0.887, 1.1)

Plots

We illustrate the rate of convergence with the help of a plot recipe for ConvergenceResult:

plt = plot(result)

For the sake of illustration, we plot some approximations of a sample target solution:

plot(suite, ns=nsample)

We can also visualize the noise associated with this sample solution,

plot(suite, xshow=false, yshow=true, label="gBm noise")

and the sine of the noise, which is the coefficient of the equation

plot(suite, xshow=false, yshow=sin, label="sin of gBm noise")

This page was generated using Literate.jl.