Monte Carlo simulation of a second order ODE?

General Usage
1

I’m interested in solving a second order ODE problem using a stochastic approach.
My ODE contains several coefficients, which up to this point, have been constant parameters (µ, Y, η, and φ). Now, I’d like to see how fluctuations in those parameters affect the solution.

All that I’ve managed has been to use the Parallel Ensemble interface to introduce fluctuations in the initial condition (a₀). I have not found an example, which does the same to arbitrary parameters in the function to be solved. There was an example using an SDE which introduced randomness in the “p” vector of parameters as given to the SDE, but I’m not sure if this is the same case.

Help appreciated!

This is what I have up to now:

using Revise
using DifferentialEquations
using Plots
using QuadGK
using Printf
using SparseArrays
using Interpolations
    
µ = 4.5e6                                          
Y = 6.0e6                                           
ρ = 1670.0                                        
η = 1.0e2                                       
φ = 0.04                         
tᵣ = 1.0e-6                         
cur_P = 0.2e9
da₀ = [0.0]                                            
a₀  = [10.0e-6]                                      
b₀ = a₀./cbrt(φ)                                     
tspan = (0.0, 10.0e-6)

slope = cur_P/tᵣ 
P = t-> t<tᵣ ? slope*t : cur_P 

α₀ = b₀.^3 ./(b₀.^3 - a₀.^3)
α₁ = (2*µ.*α₀ .+ Y)./(2*µ.+Y)
α₂ = 2*µ.*α₀/(2*µ.+Y)
     
    
function α(a, b)
    b.^3 ./(b.^3 .- a.^3)
end

function b_radius(a)
    cbrt.(b₀.^3 .- a₀.^3 .+ a.^3)
end

function c_radius(α)
    cbrt.((2*µ.*(a₀.^3)./Y).*(α₀.-α)./(α₀.-1))
end

function Pᵥ(a, c, da, α)
    if (α[1] <= α₀[1]) && (α[1] >= α₁[1])
        0
    elseif (α[1] <= α₁[1]) && (α[1] >= α₂[1])
        12*a[1]^2*da[1]*η*(-(1/3)*((1/c[1]^3 - 1/a[1]^3)))       
    elseif (α[1] <= α₂[1]) && (α[1] >= 1)
        12*a[1]^2*da[1]*η*(-(1/3)*((1/b_radius(a)[1]^3) - 1/a[1]^3))    
    else
        0
    end
end

function Pᵧ(α)
    if (α[1] <= α₀[1]) && (α[1] >= α₁[1])
        4*µ*(α₀[1]-α[1])/(3*α[1]*(α[1]-1.0))
    elseif (α[1] <= α₁[1]) && (α[1] >= α₂[1])
        2*Y/3*(1-(2*µ*(α₀[1]-α[1]))/(Y*α[1]) + log( abs((2*µ*(α₀[1]-α[1]))/(Y*(α[1]-1))) ))
    elseif (α[1] <= α₂[1]) && (α[1] >= 1)
        2*Y/3*log(α[1]/(α[1]-1))
    else
        0
    end
end

function rpnnp(dda, da, a, p, t)
    Pg0 = 101325 
    k = 1.4
    Pg = Pg0*(a₀[1]/abs.(a))^(3*k)    
    @. dda = ((Pg - p(t) + Pᵧ(α(a, b_radius(a))) - Pᵥ(a, c_radius(α(a, b_radius(a))), da, α(a, b_radius(a))))/(ρ) - 0.5*da*(4*da*(1-cbrt(φ))-da*(1-cbrt(φ)^4)))/(a*(1-cbrt(φ)))
end
    
function prob_func(prob, i, repeat)
    @. prob.u0 = [prob.u0[1], rand() * prob.u0[2]]
    prob
end


problem = SecondOrderODEProblem(rpnnp, da₀, a₀, tspan, P) 
ensemble_prob = EnsembleProblem(problem, prob_func = prob_func)
sim = solve(ensemble_prob, KenCarp4(), abstol=1e-12, reltol=1e-12, trajectories = 3)
plot(sim, idxs=(2), linealpha = 0.4)
2

It should be possible to just change prob.p instead of prob.u0 here, which will modify the parameters instead of the initial conditions for each problem in your ensemble.

1 Like
3

Thanks Sevi!

You’re correct. I also had to use a tuple for all the parameters to pass them together into the function, and then what you suggested worked well.

1 Like