From 3c7d70ebd43423220b266dab218ca6d687996d08 Mon Sep 17 00:00:00 2001 From: sotech117 Date: Thu, 1 Feb 2024 12:35:05 -0500 Subject: pull examples and complete homework 1 --- examples/FallingBall.jl | 25 +++++++++++++++++ examples/FallingBall2.jl | 44 ++++++++++++++++++++++++++++++ examples/FallingBall3.jl | 34 ++++++++++++++++++++++++ hw1/1-3.jl | 40 ++++++++++++++++++++++++++++ hw1/1-6.jl | 68 +++++++++++++++++++++++++++++++++++++++++++++++ hw1/hw1-writeup.pdf | Bin 0 -> 760944 bytes 6 files changed, 211 insertions(+) create mode 100644 examples/FallingBall.jl create mode 100644 examples/FallingBall2.jl create mode 100644 examples/FallingBall3.jl create mode 100644 hw1/1-3.jl create mode 100644 hw1/1-6.jl create mode 100644 hw1/hw1-writeup.pdf diff --git a/examples/FallingBall.jl b/examples/FallingBall.jl new file mode 100644 index 0000000..35d4b31 --- /dev/null +++ b/examples/FallingBall.jl @@ -0,0 +1,25 @@ +#!/Applications/Julia-1.8.app/Contents/Resources/julia/bin/julia + +dt = 0.01 # time step in seconds +g = 9.8 # acceleration of gravity in m/s^2 + +function dynamics(y::Float64, v::Float64, t::Float64) + for i in 1:100 + y = y + v * dt + v = v - g * dt + t = t + dt + end + + return y, v, t +end + +y0 = 10.0 # initial position in meters +v0 = 0.0 # initial velocity in m/s + +y, v, t = dynamics(y0, v0, 0.0) # evolave for 100 time steps + +println("\n\t Results") +println("final time = ", t) +println("y = ", y, " and v = ", v) +println("exact v = ", v0 - g * t) +println("exact y = ", y0 + v0 * t - 0.5 * g * t^2.0) \ No newline at end of file diff --git a/examples/FallingBall2.jl b/examples/FallingBall2.jl new file mode 100644 index 0000000..535af6d --- /dev/null +++ b/examples/FallingBall2.jl @@ -0,0 +1,44 @@ +#!/Applications/Julia-1.8.app/Contents/Resources/julia/bin/julia + +using Plots # for plotting trajectory + +g = 9.8 # acceleration of gravity in m/s^2 + +dt = 0.01 # time step in seconds +t_final = 1.0 # final time of trajectory + +steps = Int64(t_final/dt) # number of time steps + +y = zeros(steps+1) # initial array of heights in meters +v = zeros(steps+1) # initial array of velocities in m/s + +function dynamics!(y, v, t::Float64) # ! notation tells us that arguments will be modified + for i in 1:steps + y[i+1] = y[i] + v[i] * dt + v[i+1] = v[i] - g * dt + #y[i+1] = y[i] + 0.5 * (v[i] + v[i+1]) * dt + t = t + dt + end + + return t +end + +y0 = 10.0 # initial position in meters +v0 = 0.0 # initial velocity in m/s + +y[1] = y0 +v[1] = v0 +t = dynamics!(y, v, 0.0) # evolve forward in time + +println("\n\t Results") +println("final time = ", t) +println("y = ", y[steps+1], " and v = ", v[steps+1]) +println("exact v = ", v0 - g * t) +println("exact y = ", y0 + v0 * t - 0.5 * g * t^2.0) + +plot(y) # plot position as a function of time + +# energy = g * y + 0.5 * v .* v # here the mass = 1 +# println("initial energy = ", energy[1]) +# println("final energy = ", energy[steps+1]) +# plot(energy) \ No newline at end of file diff --git a/examples/FallingBall3.jl b/examples/FallingBall3.jl new file mode 100644 index 0000000..123a19a --- /dev/null +++ b/examples/FallingBall3.jl @@ -0,0 +1,34 @@ +#!/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::Float64) # ! 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 +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 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 new file mode 100644 index 0000000..cc8e38f Binary files /dev/null and b/hw1/hw1-writeup.pdf differ -- cgit v1.2.3-70-g09d2