k = 100

100

m = 1

1

ω = sqrt(k/m)

10.0

# Simple Harmonic Oscillator Problem
import OrdinaryDiffEq as ODE, Plots

#Initial Conditions

x0 = 0

0

v0 = 2

2

tspan = (0, 4*pi/ω)

(0, 1.2566370614359172)

A = x0

0

B = v0/ω

#Define the problem

0.2

function harmonicoscillator(ddu, du, u, ω, t)
    ddu .= -ω^2 * u
end

#Pass to solvers

harmonicoscillator (generic function with 1 method)

prob = ODE.SecondOrderODEProblem(harmonicoscillator, 
                                 [v0], 
                                 [x0], 
                                 tspan, 
                                 ω)

�[38;2;86;182;194mODEProblem�[0m with uType �[38;2;86;182;194mRecursiveArrayTools.ArrayPartition{Int64, Tuple{Vector{Int64}, Vector{Int64}}}�[0m and tType �[38;2;86;182;194mFloat64�[0m. In-place: �[38;2;86;182;194mtrue�[0m
timespan: (0.0, 1.2566370614359172)
u0: ([2], [0])

sol = ODE.solve(prob, ODE.DPRKN6())

	timestamp	value1	value2
1	0.0	2.0	0.0
2	0.00049975	1.99998	0.000999496
3	0.00336487	1.99887	0.00672848
4	0.0103095	1.98938	0.0205824
5	0.0216327	1.95338	0.0429287
6	0.0367291	1.86661	0.0718177
7	0.0555494	1.69928	0.105473
8	0.0786199	1.41308	0.141535
9	0.10616	0.974955	0.174627
10	0.13867	0.366124	0.19662
...

t = sol.t  # extracting time from solution

31-element Vector{Float64}:
 0.0
 0.0004997501249375313
 0.00336487244651807
 0.010309471834588583
 0.02163266731231664
 0.03672910309780761
 0.0555494478947054
 ⋮
 0.9648505252275339
 1.0166043914224208
 1.0730274208898432
 1.1322722494744484
 1.2013597923291945
 1.2566370614359172

sin.(t)
#Plot

31-element Vector{Float64}:
 0.0
 0.0004997501041354171
 0.0033648660968017747
 0.010309289211496783
 0.021630980103633787
 0.03672084556276043
 0.055520883765826645
 ⋮
 0.8219637937686642
 0.850325966986029
 0.8786500202851768
 0.905379331276634
 0.9325309554239456
 0.9510565162951535

Plots.plot(sol, idxs = [2, 1], 
           linewidth = 1,
           marker = 'o',
           title = "Simple Harmonic Oscillator",
           xaxis = "Time", 
           yaxis = "Elongation", 
           label = ["x" "v"])

Plots.plot!(t, A * cos.(ω * t) + B * sin.(ω * t), 
            lw = 3, 
            ls = :dash, 
            label = "Analytical Solution x")

Plots.plot!(t, ω * (-A * sin.(ω * t) + B * cos.(ω * t)), 
            lw = 3, 
            ls = :dash, 
            label = "Analytical Solution v")