First, lets consider the following linear model.
\[ u' = p u \]
f(u,p,t) = p.*u
f (generic function with 1 method)
We then wish to solve this model on the timespan t=0.0
to t=10.0
, with an intial condition u0=10.0
and parameter p=-0.3
. We can then setup the differential equations, solve, and plot as follows
using DifferentialEquations, Plots u0 = [10.0] p = [-0.3] tspan = (0.0,10.0) prob = ODEProblem(f,u0,tspan,p) sol = solve(prob) plot(sol) ylims!(0.0,10.0)
However, what if we wish to consider a random initial condition? Assume u0
is distributed uniformly from -10.0
to 10.0
, i.e.,
using Distributions u0_dist = [Uniform(-10.0,10.0)]
1-element Array{Distributions.Uniform{Float64},1}: Distributions.Uniform{Float64}(a=-10.0, b=10.0)
We can then run a Monte Carlo simulation of 100,000 trajectories by solving an EnsembleProblem
.
prob_func(prob,i,repeat) = remake(prob, u0 = rand.(u0_dist)) ensemble_prob = EnsembleProblem(prob,prob_func=prob_func) ensemblesol = solve(ensemble_prob,Tsit5(),EnsembleThreads(),trajectories=100000)
EnsembleSolution Solution of length 100000 with uType: DiffEqBase.ODESolution{Float64,2,Array{Array{Float64,1},1},Nothing,Nothing, Array{Float64,1},Array{Array{Array{Float64,1},1},1},DiffEqBase.ODEProblem{A rray{Float64,1},Tuple{Float64,Float64},false,Array{Float64,1},DiffEqBase.OD EFunction{false,typeof(Main.##WeaveSandBox#253.f),LinearAlgebra.UniformScal ing{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,N othing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{ },NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem},OrdinaryDiffEq.Tsi t5,OrdinaryDiffEq.InterpolationData{DiffEqBase.ODEFunction{false,typeof(Mai n.##WeaveSandBox#253.f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing, Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Not hing},Array{Array{Float64,1},1},Array{Float64,1},Array{Array{Array{Float64, 1},1},1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}},DiffEqBase.DES tats}
Plotting the first 250 trajectories produces
plot(ensemblesol, vars = (0,1), lw=1,alpha=0.1, label=nothing, idxs = 1:250)