39) Review of ODEs and practice exercises

39) Review of ODEs and practice exercises#

We will implement a ODE simulation of a simple gravity pendulum. At various points in this exercise, you will have the ability to modify the code to explore different variations on the project.

The math#

Let us consider a weight at the end of a massless arm attached to a frictionless pivot. We will only consider motion in two dimensions for this example.

The differential equation which governs the motion of this pendulum is given by

\[ \frac{d^2 \theta}{d t^2} + \frac{g}{l} \sin \left( \theta \right) = 0 \]

where \(g\) is the gravitational constant, \(l\) is the length of the arm, and \(\theta\) is the angle from vertical for the pendulum.

Note: Wikipedia has a nice pair of derivations for this equation. It is a good example of deriving the equation for our simulation from first principles (in this case, physical laws around force, torque, or energy).

Extension#

Note that we can change the constants \(g\) and \(l\) here to modify our simulation. We could also change our differential equation to represent more complex behavior.

A complication#

The time steppers for ODEs we saw in class only support first derivatives with respect to time, not second derivatives. There is a common trick we can employ though.

Let \(v = \frac{d \theta}{d t}\) and note that \(\frac{d v}{d t} = \frac{d^2 \theta}{d t^2}\). Then

\[\begin{split} \begin{split} \frac{d \theta}{d t} &= v\\ \frac{d v}{d t} &= - \frac{g}{l} \sin \left( \theta \right) \end{split} \end{split}\]

We now have a system of first order equations and can use our techniques from class.

# Some simulation constants
# Note that we can change these
g = 9.81 # m/s^2
l = 10   # m

# Function to represent our physics
function dFdt(x)
    θ = x[1]
    v = x[2]
    return [v; (-g/l)*sin(θ)]
end

dFdt (generic function with 1 method)

Initial conditions#

We needed initial conditions for our timesteppers in class. Here, we will need initial values for \(\theta\) and \(v\), so an initial displacement and initial velocity.

# Initial conditions
θ_0 = π/4 # 180°
v_0 = -1  # rad/s
x_0 = [θ_0; v_0]

2-element Vector{Float64}:
  0.7853981633974483
 -1.0

Visualize#

It would be really nice if we can visualize out results. We could plot \(\theta\) as a function of \(t\), but that wouldn’t make a nice graphic of a swinging pendulum. To make a graphic of a swinging pendulum, we need the coordinates \(\left( x, y \right)\) of our pendulum.

We should pick a coordinate system first. Then we would indicate where \(\theta\) is relative to our coordinate system and derive the coordinates of the weight from this diagram.

image.png

Given this diagram, we can now say that

\[ x = l \sin \left( \theta \right), y = - l \cos \left( \theta \right) \]

With this, we can plot the initial location of our pendulum.

using Plots
default(linewidth=4, legendfontsize=12)

x =  l * sin(θ_0)
y = -l * cos(θ_0)
plot([(0, 0), (x, y)],
    xlims=[-l, l], ylims=[-l, l],
    marker=:circle,
    legend=:none, aspect_ratio=:equal)

Time Stepper#

Now we need to pick one of our for advancing our solution forward in time. This is the area where we can experiment the most. For now, let’s use a simple forward Euler.

\[ x_{n + 1} = x_n + h \, dF \left( t_n, x_n \right) \]

where \(t_n\) is the current time, \(x_n\) is the current solution, \(h\) is our step size in time, \(dF\) is our differential equation, and \(x_{n + 1}\) is our new solution.

Let’s code up forward Euler and take a single step.

# Forward euler
function euler_forward(x, h, dF)
    return x + h * dF(x)
end

# And take a single step
h  = 1e-1
@show x_1 = euler_forward(x_0, h, dFdt)

# And plot
x =  l * sin(θ_0)
y = -l * cos(θ_0)
plot([(0, 0), (x, y)],
    xlims=[-l, l], ylims=[-l, l],
    marker=:circle,
    legend=:none, aspect_ratio=:equal)

θ_1 = x_1[1]
x =  l * sin(θ_1)
y = -l * cos(θ_1)
plot!([(0, 0), (x, y)], marker=:circle)

x_1 = euler_forward(x_0, h, dFdt) = [0.6853981633974483, -1.0693671752344003]

We can also do this with backward Euler.

using LinearAlgebra

# Backward euler
function euler_backward(x, h, dF)
    # Fixed point iteration
    x_n1 = x + h * dF(x)
    x_n2 = x + h * dF(x_n1)
    # Note - this is a simple fixed point iteration - could we do this better?
    tol = 1e-10
    while (norm(x_n1 - x_n2) > tol)
        x_n1 = x_n2
        x_n2 = x + h * dF(x_n1)
    end
    return x_n1
end

# And take a single step
h  = 1e-1
@show x_1 = euler_backward(x_0, h, dFdt)

# And plot
x =  l * sin(θ_0)
y = -l * cos(θ_0)
plot([(0, 0), (x, y)],
    xlims=[-l, l], ylims=[-l, l],
    marker=:circle,
    legend=:none, aspect_ratio=:equal)

θ_1 = x_1[1]
x =  l * sin(θ_1)
y = -l * cos(θ_1)
plot!([(0, 0), (x, y)], marker=:circle)

x_1 = euler_backward(x_0, h, dFdt) = [0.6792355369698221, -1.061626264511185]

Exploration#

Implement the \(\theta\) method.

Building it out#

Ok, now it seems like we need two things right now.

  1. Something to take our time stepper (euler_forward) and our function (dFdt) and create a series of outputs.

  2. Something to automate and animate the plotting process.

Let’s start with 1). We’ll make a basic function run_simulation to produce our values \(\theta\) and \(v\) via our ODE and time stepper.

function run_simulation(x_0, h, dF, time_stepper, t_f)
    # Setup
    num_steps = Int(ceil(t_f / h))
    x = zeros((num_steps + 1, maximum(size(x_0))))
    x[1, :] = x_0

    # Simulate
    for i in 1:num_steps
        x[i + 1, :] = time_stepper(x[i, :], h, dF)
    end
    x
end

# And run it
h   = 1e-3
t_f = 30 # s
x  = run_simulation(x_0, h, dFdt, euler_forward, t_f)

# Plotting the θ and v as a sanity check
plot([x[:, 1], x[:, 2]], xlabel="\$n\$", ylabel="\$\\theta\$", label=["\$\\theta\$" "\$v\$"])

Take a look at our plots of \(\theta\) and \(v = d \theta / d t\).

Is anything concerning to you?

What happens if we use backward Euler instead?

What about a \(\theta\) method instead?

Animation#

The plotting isn’t too bad once we have the solution array.

# We wrap the plotting in an annotated loop
animation = @animate for i in 1:maximum(size(x))
    θ_n = x[i, 1]
    x_n =  l * sin(θ_n)
    y_n = -l * cos(θ_n)
    plot([(0, 0), (x_n, y_n)],
        xlims=[-l, l], ylims=[-l, l],
        marker=:circle,
        legend=:none, aspect_ratio=:equal)
end every 100;

gif(animation, "pedulum.gif", fps=30)

[ Info: Saved animation to /home/vbarra/Dropbox/SDSU/Teaching/Comp526/fall25/lectures/pedulum.gif
Contents