# Jet transport for the simple pendulum

In this notebook we will use `TaylorSeries.jl` and `TaylorIntegration.jl` in order to integrate the the simple pendulum, defined by the Hamiltonian:

$$
H(x, p) = \frac{1}{2}p^2-\cos x 
$$

The corresponding equations of motion are:

$$
    \begin{align}
        \dot x&=p \\
        \dot p&=-\sin x
    \end{align}
$$

We will integrate this problem for the initial condition $x(t_0)=x_0$, $p(t_0)=p_0$, as well as a neighbourhood $U_0$ around this initial condition, using the jet transport technique. For simplicity, we will take $p_0=0$, and $t_0=0$. Furthermore, we will choose $x_0$ such that the pendulum librates; that is, we will choose the numerical value of the energy, $E=H(x_0,p_0)=-\cos x_0$, such that the pendulum's motion in phase space is "below" (inside) the region bounded by the separatrix. Then, for the initial conditions given above, the librational period $T$ of the pendulum is

$$
T=\frac{4}{\sqrt{2}}\int_0^{x_0}\frac{dx}{\sqrt{\cos x_0-\cos x}}
$$

Or else, in terms of the complete elliptic integral of the first kind, $K$:

$$
T=4K(\sin(x_0/2))
$$

The jet transport technique consists in considering the neighbourhood $U_0$ as being parametrized by the sum $q_0+\xi$, where $q_0=(x_0,p_0)$ represents the coordinates of the initial condition in phase space, and $\xi=(\xi_1,\xi_2)$ represents an small variation with respect to this initial condition. We re-interpret each component of the sum $q_0+\xi$ as a multivariate polynomial in the variables $\xi_1$ and $\xi_2$, and then we integrate these multivariate polynomials in time using Taylor's method. 

The packages we need to load are:


```julia
using TaylorIntegration, Elliptic
```

The Hamiltonian for the simple pendulum is:


```julia
H(x) = 0.5x[2].^2-cos(x[1])
```

The LHS of the equations of motion is:


```julia
function pendulum!(t, x, dx)
    dx[1] = x[2] 
    dx[2] = -sin(x[1])
end
```

We must setup the the `TaylorN` variables necessary to perform the jet transport.; `varorder` represents the order of expansion in the variations $\xi$.


```julia
const varorder = 8
p = set_variables("ξ", numvars=2, order=varorder)
```

The nominal initial condition is:


```julia
q0 = [1.3, 0.0]
```

The initial value of the energy is:


```julia
H0 = H(q0)
```

The parametrization of the neighbourhood $U_0$ is represented by the following array of `TaylorN`'s:


```julia
q0TN = q0 + p
```

To get a feel of how the jet transport technique works, we will evaluate the
Hamiltonian at $q_0+\xi$ in order to obtain the 8-th order Taylor expansion of
the Hamiltonian wrt the variations $\xi$, around the initial condition $q_0$:


```julia
H(q0TN)
```

Below, we setup some parameters for the Taylor integration. In particular, we will use a method of `taylorinteg` which returns the solution only at `t0`, `t0+integstep`, `t0+2integstep`,...,`tmax`, where `t0` and `tmax` are the initial and final times of integration, whereas `integstep` is the time interval between successive evaluations of the solution vector, chosen by the user. In this case, we are choosing `integstep` to be $\frac{1}{8}$th of the period of the pendulum.


```julia
const order = 28 #the order of the Taylor expansion wrt time
const abstol = 1e-20 #the absolute tolerance of the integration
const T = 4*Elliptic.K(sin(q0[1]/2)^2) #the librational period
const t0 = 0.0 #the initial time
const tmax = 6T #the final time
const integstep = T/8 #the time interval between successive evaluations of the solution vector
```

Then, we perform a Taylor integration over the initial condition `x0TN`:


```julia
#this is a "warmup lap"
@time xv = taylorinteg(pendulum!, q0TN, t0:integstep:tmax, order, abstol, maxsteps=1);
```


```julia
tv = t0:integstep:tmax;
@time xv = taylorinteg(pendulum!, q0TN, t0:integstep:tmax, order, abstol);
```

Note that we the integration above "works" for any initial neighborhood $U_0$ around the nominal initial condition $x_0$, of sufficiently small radius.

We will consider the particular case where $U_0$ is a disk of radius $r$, centered at $q_0$; that is $U_0=\{ q_0+\xi:\xi=(r\cos\phi,r\sin\phi); \phi\in[0,2\pi) \}$ for a given radius $r>0$. We will denote by $U_t$ the propagation of the initial neighbourhood $U_0$ evaluated at time $t$. Also, we denote by $q(t)$ the coordinates of the nominal solution at time $t$: $q(t)=(x(t),p(t))$. Likewise, we will denote the propagation at time $t$ of a given initial variation $\xi_0$ by $\xi(t)$. Then, we can compute the propagation of the boundary $\partial U_t$ of the neighbourhood $U_t$.


```julia
polar2cart(r, ϕ) = [r*cos(ϕ), r*sin(ϕ)]
r = 0.05 #the radius of the neighbourhood
ϕ = 0.0:0.1:(2π+0.1) #the values of the angle
ξv = polar2cart.(r, ϕ);
```

Evaluate the jet at $\partial U_{x(t)}$ (the boundary of $U_{x(t)}$) at each value of the solution vector `xv`. We organize these values such that we can plot them later:


```julia
@time xjet_plot = map(λ->λ.(ξv), xv[:,1])
@time pjet_plot = map(λ->λ.(ξv), xv[:,2]);
```

Evaluate the jet at the nominal solution, which corresponds to $\xi=(0,0)$, at each value of the solution vector `xv`:


```julia
x_nom = xv[:,1]()
p_nom = xv[:,2]();
```

Now, we plot the nominal solution (black dots), as well as the evolution of the neighbourhood $U_0$ (in colors), each $\frac{1}{8}$th of a period. The initial condition corresponds to the black dot situated at $q_0=(1.3,0)$


```julia
using Plots, LaTeXStrings
gr()
```


```julia
plot(
    xjet_plot,
    pjet_plot,
    xaxis=(L"$x$",),
    yaxis=(L"$p$",),
    # title="Simple pendulum phase space",
    leg=false,
    aspect_ratio=1.0
)
scatter!(
    x_nom,
    p_nom, 
    color=:black,
    m=(1,2.8,stroke(0))
)
#savefig("pendulum.pdf")
```


```julia

```

ae987ddbafbe3232db1b9c65cf54a17aeb948428
ipynb
Jupyter Notebook
examples/Jet-transport-simple-pendulum.ipynb
SebastianM-C/TaylorIntegration.jl
f3575ee1caba43e21312062d960613ec2ccba325
MIT
2016-09-22T22:32:12.000Z
2022-03-23T13:35:18.000Z
examples/Jet-transport-simple-pendulum.ipynb
SebastianM-C/TaylorIntegration.jl
f3575ee1caba43e21312062d960613ec2ccba325
MIT
2016-09-21T05:43:08.000Z
2022-03-15T02:55:17.000Z
examples/Jet-transport-simple-pendulum.ipynb
SebastianM-C/TaylorIntegration.jl
f3575ee1caba43e21312062d960613ec2ccba325
MIT
2016-09-24T04:37:11.000Z
2022-03-25T13:48:07.000Z
Qwen/Qwen-72B
1. YES
2. YES