I am learning how to use DifferentialEquations.jl, and wrote a simple first piece to code to simulate an undamped SHO. When I solve for the motion and plot the x and v curves using DifferentialEquations.jl, I get a beautiful smooth set of sinusoids (as expected), even though the RK4 method is using adaptive time stepping. So, it must be doing an interpolation of the returned data in order to make such a smooth plot.
The trouble I run into is that I tried to use the solutions for x (sol[1, : ]) and v (sol[2, : ] ) to then plot the kinetic, potential and total energies. These curves look horrible (because of the irregular time sampling). So, I went to Interpolations.jl and tried to resample the data to a regular grid. I get a better result, but not nearly as good as DifferentialEquations.jl.
I’m wondering if there are any experts out there that can either (a) tell me a better interpolation method to use, or (b) tell me which one DifferentialEquations.jl uses. The code is below (and it’s easiest to look at in a Jupyterlab notebook so that you can see both plots.
using Plots, Interpolations, DifferentialEquations
gr()
m = 1.0
k = 1.0
function sho!(ddu, du, u, p, t)
ddu[1] = -(k/m)*u[1]
end
du0 = [0.0]
u0 = [3.0]
tspan = (0.0,30.0)
prob = SecondOrderODEProblem(sho!,du0,u0,tspan)
sol = solve(prob, RK4())
# Plot solution (this makes a smooth plot)
plot(sol)
# now interpolate the solution data and compute the energies:
dt = 0.001
t_grid = 0:dt:30
x_interp = interpolate((sol.t,), sol[1,:], Gridded(Linear()))
x_fixed = [x_interp(t) for t in t_grid]
v_interp = interpolate((sol.t,), sol[2,:], Gridded(Linear()))
v_fixed = [v_interp(t) for t in t_grid]
# Compute the potential and kinetic energies at each time step
E_pot = 0.5 * k .* x_fixed.^2
E_kin = 0.5 * m .* v_fixed.^2
E_tot = E_pot + E_kin
# Plot the energies as a function of time (doesn't look so great!)
plot(t_grid, E_tot, xlabel="Time", ylabel="Total Energy", label="Total Energy")
plot!(t_grid, E_kin, xlabel="Time", ylabel="Total Energy", label="Kinetic Energy")
plot!(t_grid, E_pot, xlabel="Time", ylabel="Total Energy", label="Potential Energy")
Thanks in advance for any help understanding where I’m going wrong here.
-paul