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)