Computing the derivatives

We have already argued that using derivatives gives optimization algorithms a boost. There are three methods to compute derivatives (and gradients, Jacobians, Hessians):

symbolic methods,
numerical finite-difference (FD) methods,
algorithmic (also automatic) differentiation (AD) methods.

Symbolic methods

These are essentially the methods that we have all learnt to apply using a pen and paper. A bunch of rules. The input for the procedure is a function and the output from the procedure is another function. For example, for f(x) = x^2, the derivative is f'(x) = 2x.

Although straightforward and automatable, symbolic methods are not always the best choice. When does this happen?

The function to be differentiated is already rather complicated, and the derivative will typically be even more complicated. Its evaluation then may be computationally expensive. We will see that in the example.
The function to be differentiated is not available in the closed form (as a formula), but only as a software implementation, however open-source.

Example 1 (Symbolic computation of the gradient of a function of the simulated trajectory) Consider the state-space model of a pendulum \underbrace{ \begin{bmatrix} \dot \theta\\ \dot \omega \end{bmatrix}}_{\dot{\bm x}} = \underbrace{ \begin{bmatrix} \omega\\ -\frac{g}{l}\sin\theta \end{bmatrix}}_{\mathbf f(\bm x)}, where l=1\,\mathrm{m} is the length of the pendulum, g=9.81\,\mathrm{m}/\mathrm{s}^2 is the acceleration due to gravity, \theta and \omega are the angle and angular velocity of the pendulum, respectively. We are going to simulate the trajectory of the pendulum that is initially at some nonzero angle, say, \theta(0) = \pi/4 = \theta_0, and zero velocity, that is, \omega(0) = 0 = \omega_0. And we are going to consider the 2-norm (actually its square for convenience) of the state vector at the end of the simulation interval as the cost function to be minimized, for which we need to evaluate the gradient at the initial state.

First, we write an ODE solver that obtains an approximation of the final point \bm x(t_\mathrm{f}) of the state trajectory on the interval [0,t_\mathrm{f}], and a function that computes the cost as a function of the initial state J(\bm x_0).

Alternative (and still imperfect) notation

The state at the final time is a function of the state at the initial time, hence we could also write it as \bm x(t_\mathrm{f};\bm x_0), in which case the cost cold be written as J(\bm x(t_\mathrm{f};\bm x_0)). The dependence of the cost on the initial state is still visible, but the notation is a bit more clumsy (and abusive anyway).

function solve_for_final_state_fd(f, x₀, tspan, h)
    t0, tf = tspan
    t = t0
    x = x₀
    while t < tf
        x = x + h * f(x)
        t = t + h
    end
    return x
end

function J(x₀) 
    x_final = solve_for_final_state_fd(f, x₀, tspan, h)
    return x_final[1]^2 + x_final[2]^2
end

And now we use the solver to compute the trajectory and the cost

g = 9.81
l = 1.0
f(x) = [x[2], -g/l*sin(x[1])]  
θ₀ = π/4
ω₀ = 0.0   
tspan = (0.0, 1.0)
h = 0.1

J([θ₀, ω₀])

1.6184460563562304

We now use the Symbolics.jl package to compute the gradient of the cost function at the initial state. We first define symbolic state variables and and then obtain the symbolic expression for the cost function just by evaluating the function we already have at these symbolic state variables.

using Symbolics
@variables θ₀ ω₀
J_sym = J([θ₀, ω₀])

(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀)^2 + (-0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(0.1(-0.9810000000000001sin(0.1(-0.9810000000000001sin(θ₀) + ω₀) + θ₀ + 0.1ω₀) - 0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) + ω₀) + 0.1(-0.9810000000000001sin(θ₀) - 0.9810000000000001sin(θ₀ + 0.1ω₀) + ω₀) + θ₀ + 0.1ω₀) + ω₀)^2

If the shortcomings of symbolic methods have not yet started surfacing, scroll to the right in the output field. Rather long, n’est-ce pas? And we have just started, because we now must differentiate this long expression symbolically (and then we convert it from a symbolic expression back to a standard function in Julia):

∇J_sym = Symbolics.gradient(J_sym,[θ₀, ω₀])
∇J_sym_expr = build_function(∇J_sym, [θ₀, ω₀])
∇J_sym_fun = eval(∇J_sym_expr[1])

Finally, let’s do some benchmarking of evaluation of the gradient at the given (initial) state:

using BenchmarkTools
@btime ∇J_sym_fun([π/4,0.0])

  228.940 μs (4 allocations: 160 bytes)

2-element Vector{Float64}:
 3.9096539447898264
 0.020083627236371047

As a teaser for what what is to come, we also benchmark the solution based on AD:

using ForwardDiff

@btime ForwardDiff.gradient(J, [π/4,0.0])

  3.674 μs (136 allocations: 5.42 KiB)

2-element Vector{Float64}:
 3.909653944789824
 0.020083627236369115

Note that while for the former quite some work must have been done before the timing analysis was started (namely, the formula for the derivative had to be found), for the latter we started only with the function definition. And yet the latter approach wins hands down. But before we start exploring the AD methods, we give a brief overview of numerical methods based on finite-difference (FD) approximations.

Numerical finite-difference (FD) methods

These methods approximate the derivative by computing differences between the function values at different points, hence the name finite-difference (FD) methods. The simplest FD methods follow from the definition of the derivative after omiting the limit:

\frac{\mathrm d f(x)}{\mathrm d x} \approx \frac{f(x+h)-f(x)}{h}\qquad\qquad \text{forward difference} or \frac{\mathrm d f(x)}{\mathrm d x} \approx \frac{f(x)-f(x-h)}{h}\qquad\qquad \text{backward difference} or \frac{\mathrm d f(x)}{\mathrm d x} \approx \frac{f(x+\frac{h}{2})-f(x-\frac{h}{2})}{h}\qquad\qquad \text{central difference}

For functions of vector variables, the same idea applies, but now we have to compute the difference for each component of the vector.

Dependence of the error on the step size

The finite-difference methods only approximate the derivatives. The smaller the h in the above formulas, the smaller the approximation error. Really? Not quite. Let’s explore it through an example.

Example 2 (Dependence of the error on the step size) Consider a scalar function of a vector argument f(\bm x) = \sum_{i=1}^n x_i^2.

f(x) = sum(xᵢ^2 for xᵢ in x)

f (generic function with 1 method)

Now, in order to compute the gradient \nabla f, we need to compute all the individual partial derivatives, the individual components of the vector. Let’s now restrict ourselves just to one component, say, the first one, that is, let’s compute \frac{\partial f(\mathbf x)}{\partial x_1}.

In this simple example, a formula can be written down upon inspection: \frac{\partial f(\mathbf x)}{\partial x_1} = 2x_1:

x₀ = rand(10)
∇f_1_exact = 2*x₀[1]

1.0356801966104037

but let’s pretend that this answer is not available to us (we will only use it for evaluation of approximation errors of of chosen FD methods).

We now give a function for computing the first entry in the gradient (vector) by using the forward FD method. Note that in defining the function we exploit the multiple dispatch functionality of Julia, thanks to which the function will handled the floating point model of the input appropriately. That is, the input vector could be IEEE-754 double-precision floating-point format, or IEEE-754 single-precision floating-point format (or perhaps even something else).

function ∇f_1(f,x₀::Vector{T},h::T) where T<:Real
    (f(x₀[1]+h)-f(x₀[1]))/h
end

∇f_1 (generic function with 1 method)

We can now compute the first entry of the gradient for the particular vector given in the IEEE double format (default for Julia)

∇f_1(f,x₀,h)

1.1356801966104029

We can also compute the same quantity for the same vector represented in the IEEE single format:

∇f_1(f,Vector{Float32}(x₀),Float32(h))

1.1356807f0

Obviously, both answers differ from the accurate one computed above.

Now, let’s examine the error as a function of the size of the interval h:

h_range = exp10.(range(-13, stop=-1, length=1000));
abs_err_64 = [abs((∇f_1_exact - ∇f_1(f,x₀,h))) for h in h_range];
abs_err_32 = [abs((∇f_1_exact - ∇f_1(f,Vector{Float32}(x₀),Float32(h)))) for h in h_range];

using Plots
plot(h_range, abs_err_64,xaxis=:log, yaxis=:log, xlabel="h", ylabel = "|e|", label = "IEEE double")
plot!(h_range, abs_err_32,xaxis=:log, yaxis=:log, xlabel="h", ylabel = "|e|", label = "IEEE single")

Figure 1: Dependence of the absolute error of the FD approximation of the derivative on the step size h for the 64- and 32-bit floating-point formats

If we read the graph from right to left – as h is getting smaller –, we observe that initially the error decreases for both 64- and 32-bit floating-point format, and it decreases at the same rate. The reason for this decrease is that the trunction error (essentially what we commit here by doing FD approximation is that we truncate the Taylor series) dominates here, and this error goes down with h.

The major intended takeaway from this example is that this reduction of the error only takes place down to some h below which the error no longer decreases; in fact, it actally increases as h gets smaller. The reason is that for h this small, the rounding errors dominate. Apparently, they start exhibitting themselves for smaller values with the 64-bit format than with 32-bit format. The rounding errors are getting dominant here as we are subtracting two numbers that are getting more and more similar as h decreases. This is known as the phenomenon of catastrophic cancellation.

It is known from rigorous numerical analysis that the error in the case of the simple backward or forward finite difference approximation to a scalar derivative scales with \sqrt{\epsilon}, where \epsilon is the machine epsilon. Here we can expect even worse outcome as the dimension n of the vector grows. Note that \epsilon for double precision IEEE is

2^(-53)

1.1102230246251565e-16

which is available in Julia through eps() function with the type as the input argument (if no argument is given, it assumes Float64):

sqrt(eps())

1.4901161193847656e-8

Similarly, the 32-bit version is

sqrt(eps(Float32))

0.00034526698f0

Apparently, these are roughly the “cutt-off” values of h.

Automatic (also Algorithmic) differentiation (AD) methods

Here we only sketch the key ideas. Details can be found at [1, Ch. 5], or [2, Sec. 6.6], or [3, Sec. 2].

Forward AD

Example 3 We explain the essence of the Forward AD by means of an example. Consider the scalar function f(\bm x) = x_1x_2 + \sin(x_1) of a vector argument \bm x\in \mathbb R^2. We want to compute the derivatives with respect to the two inputs \nabla f(\bm x) = \begin{bmatrix}\frac{\partial f}{\partial x_1} \\ \frac{\partial f}{\partial x_2} \end{bmatrix} at a particular value of x_1 and x_2, say, x_1=2 and x_2=3.

We will find use of the computational graph of the function as is in Fig. 2.

Figure 2: Computational tree for the example

The nodes correspond to the unary or binary operations and we label the edges corresponding to the flow of data with auxiliary v_i variables. For given values of x_1 and x_2, we can evaluate the function by evaluating the values of the auxiliary variables in the order of their dependencies, that is, first v_1 and v_2, then v_3 and v_4, and finally v_5. For the the particular values of x_1 and x_2 given above, we have v_1=2, v_2=3, v_3=6, v_4=\sin(2), and v_5=6+\sin(2) as in the figure below.

Figure 3: Computational tree for the example, the values of the auxiliary variables are shown for the particular values of x_1 and x_2.

We now go after the derivatives. In forward AD we computed the derivatives \frac{\partial f}{\partial x_1} and \frac{\partial f}{\partial x_2} at the given values of x_1 and x_2 one by one.

We invoke the chain rule for differentiation of composite functions to compute the derivative \frac{\partial y}{\partial x_1} first. Namely, \begin{aligned} \frac{\partial y}{\partial x_1} &= \frac{\partial y}{\partial v_5}\frac{\partial v_5}{\partial x_1}\\ &= \frac{\partial y}{\partial v_5}\left(\frac{\partial v_5}{\partial v_3}\frac{\partial v_3}{\partial x_1} + \frac{\partial v_5}{\partial v_4}\frac{\partial v_4}{\partial x_1}\right)\\ &= \frac{\partial y}{\partial v_5}\left(\frac{\partial v_5}{\partial v_3}\left(\frac{\partial v_3}{\partial v_1}{\color{blue}\frac{\partial v_1}{\partial x_1}} + \frac{\partial v_3}{\partial v_2}{\color{blue}\frac{\partial v_2}{\partial x_1}}\right) + \frac{\partial v_5}{\partial v_4}\frac{\partial v_4}{\partial v_1}{\color{blue}\frac{\partial v_1}{\partial x_1}}\right) \end{aligned}

We can now recognize that \frac{\partial v_1}{\partial x_1}=1 and that \frac{\partial v_2}{\partial x_1}=0. We can fill these into the above formula.

Note that we can write the other derivatives \frac{\partial v_i}{\partial v_j} upon inspection of the computational graph, since each node corresponds to a simple unary or binary elementary operation.

But we structure this computation as follows. We denote \frac{\partial v_i}{\partial x} as \dot v_i. We can therefore write \dot v_1 = 1 and \dot v_2 = 0 as the first step (this is also called the seed). And referring to the computational graph, we keep determining \dot v_i for all the higher i=3, 4, 5, see the graph and the table below.

Figure 4: Computational tree for the example. Associated with each edge is not only the value v_i but also the derivative \dot v_i.

v_i	\dot v_i
v_1 = x_1 = 2	{\color{blue}\dot v_1 = 1}
v_2 = x_2 = 3	{\color{blue}\dot v_2 = 0}
v_3 = v_1\cdot v_2 = 6	\dot v_3 = \dot v_1 v_2 + v_1 \dot v_2 = 3
v_4 = \sin(v_1) = \sin(2)	\dot v_4 = \cos(v_1)\cdot \dot v_1 = \cos(2)
v_5 = v_3 + v_4 = 6+\sin(2)	\dot v_5 = \dot v_3 + \dot v_4 = 3+\cos(2)

Voila! We have obtained the value of the derivative \frac{\partial y}{\partial x_1} = \dot v_5 = 3+\cos(2) at the given values of x_1 and x_2.

If we now want to compute the other derivative \frac{\partial y}{\partial x_2}, we need to repeat the whole procedure, but now with the seed \dot v_1 = 0 and \dot v_2 = 1. We will then obtain \frac{\partial y}{\partial x_2} = \dot v_5 = 2.

Implementation of Forward AD by dual numbers

Similar to a complex number, a dual number has two components, one corresponding to the value, the other corresponding to the derivative:

x = v + d\epsilon, where the special property of \epsilon is \epsilon^2=0.

(Compare it with the property of the imaginary unit: i^2=-1.)

Multiplication of two dual numbers y = x_1\cdot x_2 is the defined naturally as \begin{aligned} x &= (v_1 + d_1\epsilon)\cdot (v_2 + d_2\epsilon)\\ &= v_1v_2 + (v_1d_2+d_1v_2) \epsilon. \end{aligned}

Similarly for other functions. We illustrate this using the following example.

Example 4 Consider a function y(x) = \cos(x^2). Its derivative is trivially \frac{\mathrm d y}{\mathrm d x} = -2x\sin(x^2).

Let’s now evaluate this result at a particular value of x, say

x = 2

First, we are going to develop our own data class (actually struct) in Julia for dual numbers. The code below is inspired by a code from [3] (they even provide Jupyter Notebooks).

Then we show the usage of a functionality already provided in ForwardDiff.jl package implementing the forward mode AD.

struct Dual
    v        # the VALUE part
    d        # the DERIVATIVE part
end

Now we need to overload the involved basic operations such as addition and multiplication of two dual numbers, multiplication by a scalar, squaring and finding the value of cosine function.

Base.:+(a::Dual,b::Dual) = Dual(a.v+b.v,a.d+b.d)
Base.:*(a::Dual,b::Dual) = Dual(a.v*b.v,a.d*b.v+b.d*a.v)
Base.:*(a::Number,b::Dual) = Dual(a*b.v,a*b.d)
Base.:^(a::Dual,b::Int) = Dual(a.v^b,b*a.v^(b-1))
Base.:cos(a::Dual) = Dual(cos(a.v),-sin(a.v)*a.d)

Let’s now check the functionality of the individual functions

X = Dual(x,1)

Dual(2, 1)

Y = Dual(3,0)

Dual(3, 0)

3*X

Dual(6, 3)

Y*X

Dual(6, 3)

X^2

Dual(4, 4)

cos(X)

Dual(-0.4161468365471424, -0.9092974268256817)

Finally, let’s use the new functionality to compute the derivative of the assigned function \cos(x^2)

cos(X^2)

Dual(-0.6536436208636119, 3.027209981231713)

In practice, you will hardly feel a need to implement your own library for algorithmic differentiation. Instead, you may want to use one of those avaialable ones, such as ForwardDiff.jl.

using ForwardDiff
X = ForwardDiff.Dual(x,1)
Y = cos(X^2)
Y.value

-0.6536436208636119

Y.partials

1-element ForwardDiff.Partials{1, Float64}:
 3.027209981231713

Compare with the “exact”

-2*x*sin(x^2)

3.027209981231713

An alternative interpretation to the above Forward AD procedure is that instead of the passes through the computational graph represented by scalar operations, we can think of the procedure as a composed vector function of vector arguments

y = f(\bm x) = e(\mathbf g(\mathbf h(\bm x))), and \begin{aligned} \bm u &= \mathbf h(\bm x) = \begin{bmatrix}x_1\\ x_2\\ \sin(x_1)\end{bmatrix},\\ \bm v &= \mathbf g(\bm u) = \begin{bmatrix}u_1\cdot u_2\\ u_3\end{bmatrix},\\ y &= e(\bm v) = v_1 + v_2. \end{aligned} The forward AD procedure then corresponds to the following application of the chain rule

\frac{\partial e(\mathbf g (\mathbf h(\boldsymbol x)))}{\partial \boldsymbol x} =\frac{\partial }{\partial \boldsymbol x}(e\circ\mathbf g \circ \mathbf h)(\boldsymbol x) = \underbrace{\frac{\partial e(\boldsymbol v)}{\partial \boldsymbol v}}_{1\times 2} \underbrace{\frac{\partial \mathbf g(\boldsymbol u)}{\partial \boldsymbol u}}_{2\times 3} \underbrace{\frac{\partial \mathbf h(\boldsymbol x)}{\partial \boldsymbol x}}_{3\times 2}.

The three terms are Jacobians of the involved functions, generally rectangular matrices. Apparently, we can proceed in two ways in order to get the product of these three. We can first compute the product of the second and the third Jacoabians

\underbrace{\frac{\partial e(\boldsymbol v)}{\partial \boldsymbol v}}_{1\times 2} \left(\underbrace{\frac{\partial \mathbf g(\boldsymbol u)}{\partial \boldsymbol u}}_{2\times 3} \underbrace{\frac{\partial \mathbf h(\boldsymbol x)}{\partial \boldsymbol x}}_{3\times 2}\right) or we can first compute the product of the first and the second Jacobians \left(\underbrace{\frac{\partial e(\boldsymbol v)}{\partial \boldsymbol v}}_{1\times 2} \underbrace{\frac{\partial \mathbf g(\boldsymbol u)}{\partial \boldsymbol u}}_{2\times 3}\right) \underbrace{\frac{\partial \mathbf h(\boldsymbol x)}{\partial \boldsymbol x}}_{3\times 2}.

The two are mathematically equivalent, but the computational cost of the two is different. We will now leave it up to you to figure out why the latter is computationally cheaper in our particular case (hint: count the number of additions and multiplications).

We only state that the forward AD corresponds to the former. The reverse AD that we are going to introduce now corresponds to the latter.

Reverse AD

The reverse AD is a bit less intuitive, but it is more efficient than the forward AD when then number of inputs is larger than the number of outputs. The extreme case is a scalar function of a vector argument, which was our case in the example above.

Similarly as in the forward AD, we introduce an auxilliary variable correspondig to the derivative, but this time the notation is different \bar v_i = \frac{\partial y}{\partial v_i}.

We first fill in all the values of v_i starting with the inputs, proceeding to the right. Once we reach the output, we set \bar v_5 = \frac{\partial y}{\partial v_5} = 1 as the seed. We then proceed to the left, filling in all the \bar v_i until we reach the inputs. The value of \bar v_1 and \bar v_2 at the end of this procedure will give us the desired derivatives \frac{\partial y}{\partial x_1} and \frac{\partial y}{\partial x_2}, respectively.

Figure 5: Computational tree for the example. Associated with each edge is not only the value v_i but also the derivative \bar v_i.

The table documents the computation. The left column is evaluated from the top to the bottom, while the right column is evaluated from the bottom to the top.

v_i	\bar v_i
v_1 = x_1 = 2	\bar v_1 \coloneqq \frac{\partial y}{\partial v_1} = \underbrace{\frac{\partial y}{\partial v_3}}_{\bar v_3} \frac{\partial v_3}{\partial v_1} + \underbrace{\frac{\partial y}{\partial v_4}}_{\bar v_4} \frac{\partial v_4}{\partial v_1} = 1\cdot 3 + 1\cdot \cos(2) = 3 + \cos(2)
v_2 = x_2 = 3	\bar v_2 \coloneqq \frac{\partial y}{\partial v_2} = \underbrace{\frac{\partial y}{\partial v_3}}_{\bar v_3} \frac{\partial v_3}{\partial v_2} = 1 \cdot v_1 = 2
v_3 = v_1\cdot v_2 = 6	\bar v_3 \coloneqq \frac{\partial y}{\partial v_3} = \underbrace{\frac{\partial y}{\partial v_5}}_{\bar v_5} \frac{\partial v_5}{\partial v_3} = 1 \cdot 1 = 1
v_4 = \sin(v_1) = \sin(2)	\bar v_4 \coloneqq \frac{\partial y}{\partial v_4} = \underbrace{\frac{\partial y}{\partial v_5}}_{\bar v_5} \frac{\partial v_5}{\partial v_4} = 1 \cdot 1 = 1
v_5 = v_3 + v_4 = 6+\sin(2)	\bar v_5 \coloneqq \frac{\partial y}{\partial v_5} = 1

Voila! We have obtained the value of the derivative \frac{\partial y}{\partial x_1} = \bar v_1 = 3+\cos(2) and \frac{\partial y}{\partial x_2} = \bar v_2 = 2 at the given values of x_1 and x_2, both as the result of the single (backward) pass through the computational graph (preceeded by the single forward pass to compute the values).

We need to emphasize that our explanation of the AD method is really just very introductory. There are many details that we have not covered. These can be found in the references.

Implementation of AD methods

Julia

ForwardDiff.jl and Enzyme.jl for forward AD. Zygote.jl and Enzyme.jl (and possibly ReverseDiff.jl that is being superseded by Enzyme.jl) for reverse AD. There are a few more but these should certainly suffice at the beginning.

It may be worth using the DifferentiationInterface.jl, which provides a common interface to different AD packages, and we can easily switch between them.

If may also be useful to check also ChainRules.jl, which support writing your own rules for differentiation of custom functions.

Some more info can be found in the discussion thread on Julia Discourse. Although it is from 2024, it is still relevant.

Matlab

Matlab has only started supporting AD relatively recently, see Automatic Differentiation Background and it is not intended to be used directly by users but rather as a backend for the optimization solvers.

There is, however, a very popular third party package for AD (not only) in Matlab, which is CasADi. It is written in C++ but has interfaces for Matlab and Python too. Very advanced, very popular in the computational/optimal control community.

References

[1]

S. Gros and M. Diehl, “Numerical Optimal Control (Draft).” Systems Control; Optimization Laboratory IMTEK, Faculty of Engineering, University of Freiburg, Apr. 2022. Available: https://www.syscop.de/files/2020ss/NOC/book-NOCSE.pdf

[2]

J. R. R. A. Martins and A. Ning, Engineering Design Optimization. Cambridge ; New York, NY: Cambridge University Press, 2022. Available: https://mdobook.github.io/

[3]

M. J. Kochenderfer and T. A. Wheeler, Algorithms for Optimization. The MIT Press, 2019. Accessed: Dec. 29, 2020. [Online]. Available: https://algorithmsbook.com/optimization/