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)