usingDifferentialEquationsfunctionparameterized_pendulum!(du,u,p,t) x,y = u α,β = p du[1] = dx = y du[2] = dy =-(α+β*cos(2*π*t))*sin(x)endnum_final=20u0 = [0.0,1.0] # Initial Conditionstspan = (0.0,num_final)p = [0.5,0.1]f=ODEFunction(parameterized_pendulum!, syms = [:x, :y])prob =ODEProblem(f,u0,tspan,p,saveat=[0,1,num_final])sol =solve(prob)
Ensemble Simulation
Now that we have tested that the system computes one solution, we can assemble the simulation experiment
When we want to simulate a bunch of solutions, called and ensemble simulation, we need to define a function which will assign new problems or initial conditions (one could also change the parameters) to the same ODEs.
A simple scenario is to modify the initial condition by a small random perturbation, to see the effect of dependence with respect to the initial conditions. To this end we take a random initial condition with a normal distribution around our the point \((\pi,0)\), which is a hyperbolic saddle pont, and build a simulation:
We can test that the function works applying it to the problem
julia> prob
ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 100.0)
u0: 2-element Vector{Float64}:
0.0
1.0
julia> prob_func(prob,1,true)
ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 100.0)
u0: 2-element Vector{Float64}:
3.1469246668695297
0.007749192580242798
and we see that the initial condition has changed from \([1,0]\) above to the random perturbation of \((\pi,0)\) (numbers will differ at every run).
One could also specify a predefined set of initial conditions to cover regularly an interval or a condition to stop based on previous computations. See the documentation at [https://docs.sciml.ai/DiffEqDocs/stable/features/ensemble/] for more details.
Now it is time to define which values will be recorded from every simulation. This can be the full solution, in which case we do not need to specify anything or a special function. For instance, in the following function we will save the first and the last values:
Althoug there is a plot recipe for plotting the solutions, we decided to save the simulation first as a DataFrame/CSV file and then loaded for plotting with your program of choice (later I will post how to do it in Julia/PlotlyJS). I think it is good practice to separate computation from postprocessing and many computing facilities only have terminal access.
We create an empty dataframe with the right column names and then fill it with the trajectories recorded
Instead of saving only the last iterate in x1 and y1, we could produce long dataframe which includes the number of the iteration, the initial condition \((x_0,y_0)\), the parameters and the computed \((x,y)\) for every iteration with is computed using the num_iterates value. In this case we need to modify the output_sol so that we return an array where each element is an array again with the saved solution:
functionoutput_func(sol,i) dim_output=length(sol.t) outsol=([sol.t, #save tmap(mod2pi,sol[1,:]), #save x sol[2,:], #save yfill(sol.prob.u0[1],dim_output), #save ic xfill(sol.prob.u0[2],dim_output), #save ic yfill(sol.prob.p[1],dim_output), #save param afill(sol.prob.p[2],dim_output), #save param bfill(i,dim_output)], #Save simulation numberfalse) #end outsolendprob =ODEProblem(f,u0,tspan,p, saveat=range(0,num_final,num_final+1))sol =solve(prob)
We check that the solution is computed at the right place and now we set the experiment
