# Interpolation problem

**URL:** <https://discourse.julialang.org/t/interpolation-problem/96530>\
**Category:** Modelling & Simulations\
**Tags:** diffeq, interpolationsjl\
**Created:** [March 23, 2023, 6:54pm UTC](https://discourse.julialang.org/t/interpolation-problem/96530 "2023-03-23T18:54:34Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![PaulN](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pauln/32/47997_2.png) [@PaulN](https://discourse.julialang.org/u/PaulN)\
**Post date:** [March 23, 2023, 6:54pm UTC](https://discourse.julialang.org/t/interpolation-problem/96530/1 "2023-03-23T18:54:35Z")

</div>

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.

```plaintext
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

---

<div class="post-metadata">

**Author:** ![contradict](https://avatars.discourse-cdn.com/v4/letter/c/ac91a4/32.png) [@contradict](https://discourse.julialang.org/u/contradict)\
**Post date:** [March 23, 2023, 8:27pm UTC](https://discourse.julialang.org/t/interpolation-problem/96530/2 "2023-03-23T20:27:06Z")

</div>

The interpolation method depends on the [solver](https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/#OrdinaryDiffEq.jl-for-Non-Stiff-Equations) used, but you can get access to it using the [solution interface](https://docs.sciml.ai/DiffEqDocs/stable/basics/solution/#Interpolations-and-Calculating-Derivatives)

For your example,

```julia
x_fixed = sol.(t_grid; idxs=1)

v_fixed = sol.(t_grid; idxs=2)

```

---

<div class="post-metadata">

**Author:** ![PaulN](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pauln/32/47997_2.png) [@PaulN](https://discourse.julialang.org/u/PaulN)\
**Post date:** [March 23, 2023, 11:47pm UTC](https://discourse.julialang.org/t/interpolation-problem/96530/3 "2023-03-23T23:47:27Z")

</div>

Thank you so much! I spent way too much time going down dead ends on my own today.  
This totally fixed my problem!

Best,  
-p