problem09

9. Canon ball. A cannon ball m is launched at angle θ and speed v0. It is acted on by gravity g and a viscous drag with magnitude \(\|c \vec{v}\|\).

(a) Find position vs time analytically.  
(b) Find a numerical solution using θ = π/4, v0 = 1 m/s, g = 1 m/s 2 , m = 1 kg, c = 1 kg/ s.  
(c) Compare the numeric and analytic solutions. At t = 2 how big is the error? How does the error depend on specified tolerances or step sizes?  
(d) Use larger and larger values of v0 and for each trajectory choose a time interval so the canon at least gets back to the ground. Plot the trajectories (using equal scale for the x and y axis. Plot all curves on one plot. As v → ∞ what is the eventual shape? [Hint: the answer is simple and interesting.]  
(e) For any given v0 there is a best launch angle θ ∗ for maximizing the range. As v0 → ∞ to what angle does θ ∗ tend? Justify your answer as best you can with careful numerics, analytical work, or both.  

../../media/problem09/problem09.png

Find position vs time analytically.

The corresponding code for this in file ./Ballistics/src/Ballistics.jl:

function analytical_sol(t)
    m, g, c = p.mass, p.gravity, p.viscosity

    u = zeros((4,))
    u[1] = vx0 * (-c/m) * (exp((-c/m)*(t-t0)) - 1)
    u[2] = ((vy0+(m*g/c))*(-c/m)*(exp((-c/m)*(t-t0)) - 1)) + ((-m*g/c)*(t-t0))
    u[3] = (vx0) * (exp((-c/m)*(t-t0)))
    u[4] = ((vy0) + (m*g/c)) * exp((-c/m)*(t-t0)) - (m*g/c)

    return u
end

This turned out to be wrong, it was an integral mistake (integration mistake). The corrected version in the same file:

function analytical_sol(t)
    m, g, c = p.mass, p.gravity, p.viscosity

    x = x0 + vx0 * ((-m / c) * exp((c / m) * (t0))) * (exp((-c / m) * t) - exp((-c / m) * t0))
    y = y0 + ((vy0 + (m * g / c)) * (-m / c) * (exp((c / m) * t0)) * (exp((-c / m) * t) - exp((-c / m) * t0))) + ((-m * g / c) * (t - t0))
    vx = (vx0) * (exp((-c / m) * (t - t0)))
    vy = ((vy0) + (m * g / c)) * exp((-c / m) * (t - t0)) - (m * g / c)

    u = zeros((4,))
    u[1] = x
    u[2] = y
    u[3] = vx
    u[4] = vy

    return u
end

TODO: Could I have created this analytical solution symbolically?

Find a numerical solution using given parameters

The corresponding code for this in file ./Ballistics/src/Ballistics.jl:

using DifferentialEquations
#...
# problem setup
x0, y0 = 0.0, 0.0
r0 = [x0; y0]
speed0 = 1
launchangle = pi / 4
vx0, vy0 = speed0 * [cos(launchangle); sin(launchangle)]
v0 = [vx0; vy0]
u0 = [r0; v0]

t0 = 0.0
tend = 1.0
tspan = (t0, tend)
p = Parameters.Param(m=1, g=1, c=1)
prob = ODEProblem(Physics.ballistic!, u0, tspan, p)
#...
sol_numeric = solve(prob)
#...
Visualization.plot_trajectory(sol_numeric.u)

This numerical solution gives the following trajectory:

../../media/problem09/numerical_trajectory.png

Compare the numeric and analytics solutions? Plot time vs error. Plot step size and tolerances vs error.

For various initial conditions, the numeric and analytical solutions on the same graph, I was shifting the analytic in x direction by 0.1 unit to be able view both, otherwise they were overlapping with the default tolerances, but then I ended up using a non-dynamic step size explicit solver:

In file ./Ballistics/src/Benchmarks.jl:

module Benchmarks

using DifferentialEquations
using LinearAlgebra

# Abs error for various step sizes
function abs_error_vs_step_sizes(prob, analytical_solve)
    # step_sizes = 10.0 .^ range(-1, -15, step=-1)
    step_sizes = 10.0 .^ range(-1, -5, step=-1)
    abs_errors_midpoint = Vector{Union{Missing, Float64}}(missing, length(step_sizes))
    abs_errors_rk4 = Vector{Union{Missing, Float64}}(missing, length(step_sizes))

    for (i,Δh) in enumerate(step_sizes)
        sol_numeric_midpoint = solve(prob, Midpoint(), dt=Δh, adaptive=false, save_everystep=false)[end]
        sol_numeric_rk4 = solve(prob, RK4(), dt=Δh, adaptive=false, save_everystep=false)[end]
        tend = prob.tspan[2]
        sol_analytic = analytical_solve(tend)

        separation_midpoint = sol_numeric_midpoint - sol_analytic
        separation_rk4 = sol_numeric_rk4 - sol_analytic
        error_midpoint = norm(separation_midpoint)
        error_rk4 = norm(separation_rk4)
    
        abs_errors_midpoint[i] = error_midpoint
        abs_errors_rk4[i] = error_rk4
    end

    log_errors_midpoint = log10.(abs_errors_midpoint)
    log_errors_rk4 = log10.(abs_errors_rk4)
    log_errors_dict = Dict("midpoint"=>log_errors_midpoint, "rk4"=> log_errors_rk4)
    log_steps = log10.(step_sizes)

    return log_steps, log_errors_dict
end

# Abs error for all times
function error_vs_time(prob, analytical_solve)
    N = 5
    step_sizes = 10.0 .^ range(-1, -N, step=-1)
    time_histories = Vector{Union{Missing, Vector{Float64}}}(missing, length(step_sizes))
    error_histories = Vector{Union{Missing, Vector{Float64}}}(missing, length(step_sizes))
    for (i,Δh) in enumerate(step_sizes)
        sol = solve(prob, Midpoint(), dt=Δh, adaptive=false, save_everystep=true)

        timestamps = sol.t
        sol_numeric_midpoint = sol.u
        sol_analytic = analytical_solve.(timestamps)

        separations = sol_numeric_midpoint .- sol_analytic
        manhatten_errors = norm.(separations, Ref(1))
        error_histories[i] = log10.(manhatten_errors)
        time_histories[i] = timestamps
    end
    return time_histories, error_histories
end


end # module Benchmarks

../../media/problem09/numeric_and_analytical_trajectories.png

Now, I will plot the steps sizes vs error (euclidean norm of state vector at tend):

../../media/problem09/steps_vs_error.png

For a better method, like RK4, on the same graph:

../../media/problem09/steps_vs_error_two_methods_compare.png

Now, for time vs error compilation:

../../media/problem09/time_vs_running_error.png

Now, I will change the relative tolerance and absolute tolerance. In file ./Ballistics/src/Benchmarks.jl:

function tolerances_vs_error(prob, analytical_solve)
    N = 15
    error_matrix = Matrix{Float64}(undef, N, N)
    num_granularity = N+7
    abstols = 10.0 .^ range(7, -N, length=num_granularity)
    reltols = 10.0 .^ range(7, -N, length=num_granularity)

    num_samples = length(abstols)*length(reltols)

    xs = Vector{Float64}(undef, num_samples)
    ys = Vector{Float64}(undef, num_samples)
    zs = Vector{Float64}(undef, num_samples)

    for (i, abstol) in enumerate(abstols)
        for (j, reltol) in enumerate(reltols)
            sol_numeric = solve(prob, Midpoint(), adaptive=true, save_everystep=false, abstol=abstol, reltol=reltol)[end]
            sol_analytic = analytical_solve(prob.tspan[2])

            separation = sol_numeric - sol_analytic
            error = norm(separation)
            xs[i+num_granularity*(j-1)] = abstol
            ys[i+num_granularity*(j-1)] = reltol
            zs[i+num_granularity*(j-1)] = error
        end
    end
    xs,ys,zs
end

Just by mistake I ran this for large (>1.0) tolerances. Below is the plot:

../../media/problem09/large_tolerances.png

Similarly for tiny (<<1.0) tolerances the plot:

../../media/problem09/small_tolerances.png

Both scales on the same graph:

../../media/problem09/tolerances_vs_error_combined.png

Simply put, it seems reducing both abstol and reltol generally reduces error from the analytic solution very quickly. What is the meaning of both is not yet entirely clear.

Use larger and larger values of v0, plot all of their trajectory till ball hits ground. What happens to the eventual shape as v -> ∞ ?

The code for this in file ./Ballistics/src/Ballistics.jl:

# ...

function solve_till_empty(prob)
    conditionatbeggining(u, t, integrator) = t > 0
    removetstop!(integrator) = add_tstop!(integrator, ∞)
    removetstopwhenbeginning = DiscreteCallback(conditionatbeggining,removetstop!)

    condition(u, t, integrator) = u[2] < 0
    affect!(integrator) = terminate!(integrator)
    stopwhenonground = DiscreteCallback(condition, affect!)

    cb = CallbackSet(removetstopwhenbeginning, stopwhenonground)
    sol = solve(prob, RK4(), callback = cb, tstops = [0.1,])

    sol
end

function problem_setup(launchangle, speed0, c)
    x0, y0 = 0.0, 0.0
    r0 = [x0; y0]
    vx0, vy0 = speed0 * [cos(launchangle); sin(launchangle)]
    v0 = [vx0; vy0]
    u0 = [r0; v0]

    t0 = 0.0
    tend = 1.0
    tspan = (t0, tend)
    p = Parameters.Param(m = 1, g = 1, c = c)
    prob = ODEProblem(Physics.ballistic!, u0, tspan, p)
    return prob
end

speeds = range(0, 100)
problems = problem_setup.(Ref(pi/4), speeds, Ref(1))
sols = solve_till_empty.(problems)
Visualization.plot_all_trajectories(sols)

# ...

The plot corresponding to this:

Increasing intial velocities

Plot speed verses best launch angle, numerically or analyticaly or both

I expect 45 degrees to still be consistantly good for launch even with drag.

Doing this numerically, in file ./Ballistics/src/Benchmarks.jl:

# ...

function best_angle_for_speed(speed0, c)
    launchangles = pi/4 .* ((range(-1.0, 1.0, 300)) .^ 3.0 .+ 1.0)

    max_range = 0
    argmax_theta = 0
    for launchangle in launchangles
        # problem setup
        prob = problem_setup(launchangle, speed0, c)

        sol = solve_till_empty(prob)
        cur_range = sol.u[end][1]
        if cur_range > max_range
            max_range = cur_range
            argmax_theta = launchangle
        end
    end
    argmax_theta
end

speeds = 10.0 .^ range(0.0, 3.1, 100)
bestangles = []
for speed0 in speeds
    bestangle = best_angle_for_speed(speed0, 1)
    push!(bestangles, bestangle)
end
Visualization.plot_best_angle_for_speeds(speeds, bestangles)

# ...

The weird choice for launchangles is for concentrating numeric choices for launchangles near \(\pi \over 4\). This produces:

../../media/problem09/speed_vs_bestangle_c_1.png

../../media/problem09/speed_vs_bestangle_c_10.png

It seems that my intuition was way off. It seems the best angle quickly becomes closer to zero as the initial velocity increases.

Let us plot the relation of drag with the best angle for a high velocity (say v = 1200), in file ./Ballistics/src/Ballistics.jl:

speed0 = 1200.0
bestangles = []
drags = 2.0 .^ range(1.0, 5.0)
for c in drags
    bestangle = best_angle_for_speed(speed0, c)
    push!(bestangles, bestangle)
end
Visualization.plot_best_angle_vs_drag(drags, bestangles)

TODO: There isn’t much change in the eventual best angle for large velocities as drag changes. It does not seem right, but the best angle in such case seems to be some constant close to zero. I will explore this more later. The plot generated by above code is:

../../media/problem09/best_angle_vs_drag.png

This concludes my attempt of problem09.