Simple example:
function pendulum!(du,u,p,t)
g = 9.81
L = 2
alpha = u[1]
omega = u[2]
du[1] = omega
du[2] = -g/L*sin(alpha)
end
#
u0 = [pi/4,0.1]
tspan = (0.,5)
prob = ODEProblem(pendulum!,u0,tspan)
sol = solve(prob)
Here, sol
is a fancy data structure that contains the time instances for the solution in field t
, sol.t
. It also contains the solution vector of variable 1
in sol[1,:]
, so I can plot the first variable as follows:
plot(sol.t,sol[1,:])
However, sol
also contains an interpolation function, and by using function syntax, sol(0.5)
contains the state vector at time = 0.5. A better way to do the plotting is the following:
plot(sol,vars=(0,1))
which (i) plots variable 1 as a function of variable 0 (= time), and (ii) uses interpolation so the solution looks smoother than if you use plot(sol.t,sol[1,:])
.
You can extract the solutions as a matrix with variable 1 in row 1, variable 2 in row 2, etc.:
sol[1:end,:]
Here, column n
corresponds to the time instance sol.t[n]
.
You can achieve the same by:
reduce(hcat,sol.(sol.t))
Here, sol.(sol.t))
computes the solution in the points given by time vector sol.t
, and reduce(hcat,.)
wraps the result into the same matrix as sol[1:end,:]
.
Finally, the solver decides at which points you want to compute the solution, i.e., the time instances in sol.t
. If you, for some reason, want to compute the solution at other time instances, say at regular time intervals tsamp = tspan[1]:0.1:tspan[2]
for my pendulum, you can create the solution matrix in those points by:
reduce(hcat,sol.(tsamp))