1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
|
#!/Applications/Julia-1.8.app/Contents/Resources/julia/bin/julia
using Plots # for plotting trajectory
using DifferentialEquations # for solving ODEs
g = 9.8 # acceleration of gravity in m/s^2
t_final = 1.0 # final time of trajectory
p = 0.0 # parameters (not used here)
function tendency!(dyv::Vector{Float64}, yv::Vector{Float64}, p, t) # ! notation tells us that arguments will be modified
y = yv[1] # 2D phase space; use vcat(x, v) to combine 2 vectors
v = yv[2] # dy/dt = v
a = -g # dv/dt = -g
dyv[1] = v
dyv[2] = a
println("t = ", t, " y = ", y, " v = ", v)
end
y0 = 10.0 # initial position in meters
v0 = 0.0 # initial velocity in m/s
yv0 = [y0, v0] # initial condition in phase space
tspan = (0.0, t_final) # span of time to simulate
prob = ODEProblem(tendency!, yv0, tspan, p) # specify ODE
sol = solve(prob, Tsit5(), reltol = 1e-8, abstol = 1e-8) # solve using Tsit5 algorithm to specified accuracy
println("\n\t Results")
println("final time = ", sol.t[end])
println("y = ", sol[1, end], " and v = ", sol[2, end])
println("exact v = ", v0 - g * t_final)
println("exact y = ", y0 + v0 * t_final - 0.5 * g * t_final^2.0)
plot(sol, idxs = (1)) # plot position as a function of time
|
