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
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
|
using Plots # for plotting trajectory
using DifferentialEquations # for solving ODEs
# INITIAL CONDITIONALS AND PARAMETERS
a = 10.0 # birth of new members
b = 3.0 # death of members
n0 = 10.0 # initial number of members
t_final = 1.0 # final time of simulation
p = 0.0 # parameters (not used here)
dt = 0.01 # time step for euler's
# EULER'S METHOD -> APPROXIMATE ANSWER
steps = Int64(t_final/dt) # number of time steps
n = zeros(steps+1) # initial array of members
v = zeros(steps+1) # initial array of rate of change of members
t = zeros(steps+1) # initial array of time intervals
function dynamics!(n::Vector{Float64}, v::Vector{Float64}, t::Vector{Float64})
for i in 1:steps
# equation: dn = dt(aN - bN^2)
dn = dt*(a*n[i] - b*n[i]*n[i])
vn = dt*(a-2.0*b*n[i])
n[i+1] = n[i] + dn
v[i+1] = v[i] + vn
t[i+1] = t[i] + dt
end
end
# calcuate with current dt, store into arrays
n[1] = n0
v[1] = 0.0
t[1] = 0.0
dynamics!(n, v, t)
# USING ODE SOLVER -> EXACT ANSWER
function tendency!(dnv::Vector{Float64}, nv::Vector{Float64}, p, t::Float64) # ! notation tells us that arguments will be modified
n = nv[1] # 2D phase space; use vcat(x, v) to combine 2 vectors
v = nv[2] # dn/dt = v
dnv[1] = a*n - b*n*n
dnv[2] = a - 2.0*b*n
end
i0 = [n0, v0] # set initial conditions
tspan = (0.0, t_final) # span of time to simulate
prob = ODEProblem(tendency!, i0, tspan, p) # specify ODE
sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8) # solve using Tsit5 algorithm to specified accuracy
n_exact = sol[1, :] # extract the population values over time
# PLOTTING AND COMPARISON
println("Parameters (a, b, n0, t_final): ", a, ", ", b, ", ", n0, ", ", t_final)
println("Final population (at ", t_final, ") via Euler's Method:\t", n[end])
println("Final population (at ", t_final, ") via ODE Solver:\t", n_exact[end])
plot_title = "Population v. time (w/ n0,a,b= " * string(n0) * ", " * string(a) * ", " * string(b) * ")"
plot(t, n, label="Euler's Method (dt = .01)", title=plot_title, lw=2, xlabel="time", ylabel="population")
plot!(sol.t, n_exact, label="Exact Solution (ode solver)", lw=2) # plot!() to add to existing plot
# NOTE: uncomment the two lines below if you also want to plot the next derativate
# plot!(t, v, label="d^2n/dt^2 (Euler's Method)", lw=2)
# plot!(sol.t, sol[2, :], label="d^2n/dt^2(ode solver)", lw=2)
|
