{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# The Lyapunov spectrum of the Lorenz system"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here, we present a calculation of the Lyapunov spectrum of the Lorenz system, using `TaylorIntegration.jl`. Our calculation involves evaluating the 1st order variational equations $\\dot \\xi = J \\cdot \\xi$ for this system, where $J = \\operatorname{D}f$ is the Jacobian. By default, the numerical value of the Jacobian is computed using automatic differentiation techniques implemented in `TaylorSeries.jl`. Automatic differentation helps us to avoid writing down manually the Jacobian; instead, `TaylorIntegration.jl` does that for us! For complex problems, or when the number of dependent variables is large, this may help save typos, which sometimes are hard to trace. Otherwise, if performance is critical, then the user may provide the Jacobian function, as will be shown below."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The Lorenz system is the ODE defined as:\n",
"\n",
"$$\n",
"\\begin{align}\n",
" \\begin{split}\n",
" \\dot x_1 &= \\sigma(x_2-x_1) \\\\\n",
" \\dot x_2 &= x_1(\\rho-x_3)-x_2 \\\\\n",
" \\dot x_3 &= x_1x_2-\\beta x_3 \n",
" \\end{split}\n",
"\\end{align}\n",
"$$\n",
"\n",
"where $\\sigma$, $\\rho$ and $\\beta$ are constant parameters."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## The setup"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"First, we write a Julia function which evaluates (in-place) the Lorenz system equations:"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"lorenz! (generic function with 1 method)"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"#Lorenz system ODE:\n",
"function lorenz!(t, x, dx)\n",
" dx[1] = σ*(x[2]-x[1])\n",
" dx[2] = x[1]*(ρ-x[3])-x[2]\n",
" dx[3] = x[1]*x[2]-β*x[3]\n",
" nothing\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next, we define the parameters:"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"45.92"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"#Lorenz system parameters\n",
"#we use the `const` prefix in order to help the compiler speed things up\n",
"const σ = 16.0\n",
"const β = 4.0\n",
"const ρ = 45.92"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The initial conditions, the initial and final time are:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"100.0"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"const x0 = [19.0, 20.0, 50.0] #the initial condition\n",
"const t0 = 0.0 #the initial time\n",
"const tmax = 100.0 #final time of integration"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We know that the sum of the Lyapunov spectrum has to be equal to the trace of the Jacobian of the equations of motion. We will calculate this trace using `TaylorSeries.jl`, and after the numerical integration, we will come back to this value to check that this is indeed the case:"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-21.0"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Note that TaylorSeries.jl exported variables and methods are @reexport-ed by TaylorIntegration.jl\n",
"#Calculate trace of Lorenz system Jacobian via TaylorSeries.jacobian:\n",
"import LinearAlgebra: tr\n",
"using TaylorIntegration\n",
"xi = set_variables(\"δ\", order=1, numvars=length(x0))\n",
"x0TN = [ x0[1]+xi[1], x0[2]+xi[2], x0[3]+xi[3] ]\n",
"dx0TN = similar(x0TN)\n",
"lorenz!(t0, x0TN, dx0TN)\n",
"lorenztr = tr(TaylorSeries.jacobian(dx0TN)) #trace of Lorenz system Jacobian matrix"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# The integration"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now, we are ready to perform the integration, using `TaylorIntegration.jl`. Usually, we would use the `taylorinteg` method in order to integrate the equations of motion; but, since we are interested in evaluating the Lyapunov spectrum, we will use the `lyap_taylorinteg` method, which calculates the Lyapunov spectrum via the variational equations and Oseledets' theorem. The expansion order will be $28$; the local absolute tolerance will be $10^{-20}$. `lyap_taylorinteg` will return three arrays: one with the evaluation times, one with the value of the dependent variables at the time of evaluation, and another one with the vaues of the Lyapunov spectrum at each time of the evaluation."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 3.038651 seconds (15.28 M allocations: 1.371 GiB, 13.00% gc time)\n"
]
}
],
"source": [
"@time tv, xv, λv = lyap_taylorinteg(lorenz!, x0, t0, tmax, 28, 1e-20; maxsteps=2000000);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As explained above, instead of using automatic differentiation to compute the numerical value of the Jacobian, the user may provide a function which computes the Jacobian in-place, as follows:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"lorenz_jac! (generic function with 1 method)"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"#Lorenz system Jacobian (in-place):\n",
"function lorenz_jac!(jac, t, x)\n",
" jac[1,1] = -σ+zero(x[1]); jac[1,2] = σ+zero(x[1]); jac[1,3] = zero(x[1])\n",
" jac[2,1] = ρ-x[3]; jac[2,2] = -1.0+zero(x[1]); jac[2,3] = -x[1]\n",
" jac[3,1] = x[2]; jac[3,2] = x[1]; jac[3,3] = -β+zero(x[1])\n",
" nothing\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Then, we may perform the same integration as before using `lorenz_jac!`:"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 0.856854 seconds (4.35 M allocations: 575.767 MiB, 20.75% gc time)\n"
]
}
],
"source": [
"@time tv_, xv_, λv_ = lyap_taylorinteg(lorenz!, x0, t0, tmax, 28, 1e-20, lorenz_jac!; maxsteps=2000000);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In general, we may check the consistency of our integration checking that both solutions are either equal, or at least approximately equal with differences comparable to roundoff errors. In the particular case of the integrations performed above, we obtain the same solution, *exactly*:"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(true, true, true)"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"tv == tv_, xv == xv_, λv == λv_"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The number of steps taken is:"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"4093"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"length(tv)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Is the final time actually the requested value?"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"true"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"tv[end] == tmax"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"What is the minimum, maximum and mean time-step?"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(0.009471209659452029, 0.0608034582365633, 0.024437927663734114)"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import Statistics: mean\n",
"minimum(diff(tv)), maximum(diff(tv)), mean(diff(tv))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"What is the standard deviation of the time-step size distribution?"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.007428238447533755"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import Statistics: std\n",
"std(diff(tv))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"What is the final value for each of the exponents?"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3-element Array{Float64,1}:\n",
" 1.4856676067102397 \n",
" -0.0012149610701265275\n",
" -22.484452645640157 "
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"λv[end,:]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"What is the value of the sum of the value we obtained for the spectrum?"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-21.000000000000043"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sum(λv[end,:])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now, is the sum of the exponents exactly equal to the trace of the Jacobian?"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"false"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sum(λv[end,:]) == lorenztr"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"So the sum is not *exactly* equal, but is it approximately equal? What is the relative error in the computed vs the expected value?"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2.0301221021717147e-15"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"rel_error = (sum(λv[end,:])-lorenztr)/lorenztr"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"How does this error compare to the machine epsilon?"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"9.142857142857142"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"rel_error/eps()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"So the relative error in our computation for the sum of the Lyapunov spectrum is comparable to the machine epsilon."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Visualization"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally, we plot our results:"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"using Plots\n",
"#gr()\n",
"#using LaTeXStrings\n",
"#plotlyjs()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Lyapunov exponents vs time plot"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {
"scrolled": false
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(tv, λv[:,1], label=\"L_1\")\n",
"plot!(tv, λv[:,2], label=\"L_2\")\n",
"plot!(tv, λv[:,3], label=\"L_3\")\n",
"xlabel!(\"time\")\n",
"ylabel!(\"L_i, i=1,2,3\")\n",
"title!(\"Lyapunov exponents vs time\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Lyapunov exponents vs time plot (semi-log)"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {
"scrolled": false
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(tv[2:end], λv[:,1][2:end], xscale=:log10, label=\"L_1\")\n",
"plot!(tv[2:end], λv[:,2][2:end], label=\"L_2\")\n",
"plot!(tv[2:end], λv[:,3][2:end], label=\"L_3\")\n",
"xlabel!(\"time\")\n",
"ylabel!(\"L_i, i=1,2,3\")\n",
"title!(\"Lyapunov exponents vs time (semi-log)\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Convergence of Lyapunov exponents vs time plot (log-log)"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {
"scrolled": false
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(tv[2:end], abs.(diff(λv[:,1])), yscale=:log10, xscale=:log10, label=\"dL_1\")\n",
"plot!(tv[2:end], abs.(diff(λv[:,2])), label=\"dL_2\")\n",
"plot!(tv[2:end], abs.(diff(λv[:,3])), label=\"dL_3\")\n",
"xlabel!(\"time\")\n",
"ylabel!(\"dL_i, i=1,2,3\")\n",
"title!(\"Convergence of Lyapunov exponents vs time (log-log)\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Absolute difference wrt trace of Jacobian vs time plot (log-log)"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {
"scrolled": false
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(tv[2:end], log10.(abs.(λv[2:end,1]+λv[2:end,2]+λv[2:end,3].-lorenztr)), xscale=:log10, leg=false)\n",
"#ylims!(-15,-12)\n",
"xlabel!(\"time (log10)\")\n",
"ylabel!(\"Absolute difference wrt Tr(J) (log10)\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## The Lorenz attractor"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(xv[:,1], xv[:,2], xv[:,3], leg=false)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 1.0.1",
"language": "julia",
"name": "julia-1.0"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.0.1"
}
},
"nbformat": 4,
"nbformat_minor": 1
}