diff options
Diffstat (limited to 'hw1')
-rw-r--r-- | hw1/1-3.jl | 40 | ||||
-rw-r--r-- | hw1/1-6.jl | 68 | ||||
-rw-r--r-- | hw1/hw1-writeup.pdf | bin | 0 -> 760944 bytes |
3 files changed, 108 insertions, 0 deletions
diff --git a/hw1/1-3.jl b/hw1/1-3.jl new file mode 100644 index 0000000..ad3f1b0 --- /dev/null +++ b/hw1/1-3.jl @@ -0,0 +1,40 @@ +using Plots # for plotting trajectory +using DifferentialEquations # for solving ODEs + + +# INITIAL CONDITIONALS AND PARAMETERS +a = 10.0 # 1st drag coefficient +b = 1.0 # 2nd drag coefficient (on v) +v0 = 0.0 # initial velocity in m/s +t_final = 10.0 # time of trajectory in s + + +# EULER'S METHOD +dt = 0.01 # time step +steps = Int64(t_final/dt) # number of time steps + +v = zeros(steps+1) # initial array of velocities +t = zeros(steps+1) # initial array of time intervals + +function dynamics!(v::Vector{Float64}, t::Vector{Float64}) + for i in 1:steps + # equation: dv = dt(a - b*v) + dv = dt*(a - b*v[i]) + + v[i+1] = v[i] + dv + t[i+1] = t[i] + dt + end +end + +# do the calculation, store into arrays +v[1] = v0 +t[1] = 0.0 +dynamics!(v, t) + +# print the parameters & terminal velocity +println("Parameters (a, b, v0, t_final): ", a, ", ", b, ", ", v0, ", ", t_final) +println("Terminal velocity: ", v[end]) + +# plot the velocities over time +plot(t, v, xlabel="time (s)", ylabel="velocity (m/s)", title="Velocity vs. time (affected by drag, w/ a,b=10,3)", label="", lw=2, color=:blue) + diff --git a/hw1/1-6.jl b/hw1/1-6.jl new file mode 100644 index 0000000..4704add --- /dev/null +++ b/hw1/1-6.jl @@ -0,0 +1,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)
\ No newline at end of file diff --git a/hw1/hw1-writeup.pdf b/hw1/hw1-writeup.pdf Binary files differnew file mode 100644 index 0000000..cc8e38f --- /dev/null +++ b/hw1/hw1-writeup.pdf |