31:T5f2,

The orbit of a planet that would orbit alone around the Sun would be an immutable ellipse having the Sun located at one of its foci. This orbit would be completely determined by the following parameters: (i) the semi-major axis ${tex`a`} of the ellipse, (ii) the eccentricity ${tex`e`} of the ellipse, (iii) the longitude of the periapsis ${tex`\varpi`} (i.e., longitude of the perihelion in the case of the solar system), (iv) the inclination ${tex`I`} of its plane, and (v) the longitude of the ascending node ${tex`\Omega`}.

However, the presence of several planets around the Sun makes the precise determination of their motions much less simple. According to Newton's law of universal graviation, like the Sun, the planets act on each other proportionally to their masses and in inverse proportion to the squares of their respective distances.

When considered over very long periods of time – over centuries or millions of years – interplanetary gravitational interactions have a profound influence on planetary orbits. In this notebook, a theory of celestial mechanics known as the orbital perturbation theory is used to examine the secular interactions between Jupiter and Saturn. This theory is also commonly referred to as Laplace-Lagrange secular perturbation theory.

32:T262f,illustration = { const boxHeight = 420; const svg = d3.create('svg').attr('viewBox', [0, 0, width, boxHeight]); // Define the circle variables const X_offset = 80, Y_offset = 170; const X_origin = 430 + X_offset, Y_origin = 220 + Y_offset; const X_Sun = 0 + X_offset, Y_Sun = 0 + Y_offset; const X_Jupiter = 450 + X_offset, Y_Jupiter = 50 + Y_offset; const X_Saturn = 600 + X_offset, Y_Saturn = -120 + Y_offset; const radius_Sun = 40, radius_Jupiter = 25, radius_Saturn = 20; const psi_start = 100*Math.PI/180, psi_end= 75*Math.PI/180 ; const arc_radius = 105; const color_arrow_origin = '#B0B0B0', color_arrow = 'black', color_arc = 'purple'; // Define the lines variables const pos_Origin = [X_origin, Y_origin] const pos_Sun = [X_Sun, Y_Sun]; const pos_Jupiter = [X_Jupiter, Y_Jupiter]; const pos_Saturn = [X_Saturn, Y_Saturn]; const lineEnd_Sun = [X_Sun + radius_Sun -1, Y_Sun + radius_Sun -20]; const lineEnd_Jupiter1 = [X_Jupiter - radius_Jupiter -5, Y_Jupiter -3.3]; const lineEnd_Jupiter2 = [X_Jupiter -3.5, Y_Jupiter + radius_Jupiter + 3]; const lineEnd_Saturn1 = [X_Saturn - radius_Saturn -5, Y_Saturn +5]; const lineEnd_Saturn2 = [X_Saturn -16, Y_Saturn + radius_Saturn -2]; const lineEnd_Saturn3 = [X_Saturn -11, Y_Saturn + radius_Saturn + 2]; // Arrowhead marker variables const markerBoxWidth = 50; const markerBoxHeight = 12; const refX = markerBoxWidth / 2; const refY = markerBoxHeight / 2; const markerWidth = markerBoxWidth / 2; const markerHeight = markerBoxHeight / 2; const arrowPoints = [[0, 0], [refX, refY], [0, markerBoxHeight]]; // First circle: Sun svg .append('circle') .attr('cx', X_Sun) .attr('cy', Y_Sun) .attr('r', radius_Sun) .style('fill', color_Sun); // Planet: Jupiter svg .append('circle') .attr('cx', X_Jupiter) .attr('cy', Y_Jupiter) .attr('r', radius_Jupiter) .style('fill', color_Jupiter); // Planet: Saturn svg .append('circle') .attr('cx', X_Saturn) .attr('cy', Y_Saturn) .attr('r', radius_Saturn) .style('fill', color_Saturn); // Add the arrowhead marker definition to the svg element svg .append('defs') .append('marker') .attr('id', 'arrow_origin') .attr('viewBox', [0, 0, markerBoxWidth, markerBoxHeight]) .attr('refX', refX) .attr('refY', refY) .attr('markerWidth', markerBoxWidth) .attr('markerHeight', markerBoxHeight) .attr('orient', 'auto-start-reverse') .append('path') .attr('d', d3.line()(arrowPoints)) .attr('fill', color_arrow_origin) .attr('stroke', color_arrow_origin); svg .append('defs') .append('marker') .attr('id', 'arrow') .attr('viewBox', [0, 0, markerBoxWidth, markerBoxHeight]) .attr('refX', refX) .attr('refY', refY) .attr('markerWidth', markerBoxWidth) .attr('markerHeight', markerBoxHeight) .attr('orient', 'auto-start-reverse') .append('path') .attr('d', d3.line()(arrowPoints)) .attr('fill', color_arrow) .attr('stroke', color_arrow); // Add the lines svg .append('path') .attr('d', d3.line()([pos_Origin, pos_Sun])) .attr('stroke', color_arrow_origin) .attr('opacity', 0.6) .attr('fill', 'none'); svg .append('path') .attr('d', d3.line()([pos_Origin, pos_Jupiter])) .attr('stroke', color_arrow_origin) .attr('opacity', 0.6) .attr('fill', 'none'); svg .append('path') .attr('d', d3.line()([pos_Origin, pos_Saturn])) .attr('stroke', color_arrow_origin) .attr('opacity', 0.6) .attr('fill', 'none'); svg .append('path') .attr('d', d3.line()([pos_Sun, pos_Jupiter])) .attr('stroke', color_arrow) .attr('opacity', 0.6) .attr('fill', 'none'); svg .append('path') .attr('d', d3.line()([pos_Sun, pos_Saturn])) .attr('stroke', color_arrow) .attr('opacity', 0.6) .attr('fill', 'none'); svg .append('path') .attr('d', d3.line()([pos_Jupiter, pos_Saturn])) .attr('stroke', color_arrow) .attr('opacity', 0.6) .attr('fill', 'none'); // Add the lines with arrowhead at the end svg .append('path') .attr('d', d3.line()([pos_Origin, lineEnd_Sun])) .attr('stroke', color_arrow_origin) .attr('marker-end', 'url(#arrow_origin)'); svg .append('path') .attr('d', d3.line()([pos_Origin, lineEnd_Jupiter2])) .attr('stroke', color_arrow_origin) .attr('marker-end', 'url(#arrow_origin)'); svg .append('path') .attr('d', d3.line()([pos_Origin, lineEnd_Saturn3])) .attr('stroke', color_arrow_origin) .attr('marker-end', 'url(#arrow_origin)'); svg .append('path') .attr('d', d3.line()([pos_Sun, lineEnd_Jupiter1])) .attr('stroke', color_arrow) .attr('marker-end', 'url(#arrow)') .attr('fill', 'none'); svg .append('path') .attr('d', d3.line()([pos_Sun, lineEnd_Saturn1])) .attr('stroke', color_arrow) .attr('marker-end', 'url(#arrow)') .attr('fill', 'none'); svg .append('path') .attr('d', d3.line()([pos_Jupiter, lineEnd_Saturn2])) .attr('stroke', color_arrow) .attr('marker-end', 'url(#arrow)') .attr('fill', 'none'); // Add text svg .append("text") .attr("x",pos_Origin[0]) .attr("y",pos_Origin[1]) .attr("dx", -10) .attr("dy", 25) .attr("font-size", '25px') .attr("font-weight", 500) .attr("font-family", 'Arial') .attr("fill", color_arrow_origin) .text('O'); svg .append("text") .attr("x",pos_Sun[0]) .attr("y",pos_Sun[1]) .attr("dx", -72) .attr("dy", 5) .attr("font-size", '23px') .attr("font-weight", 300) .attr("font-family", 'Arial') .attr("fill", color_Sun) .text('m₀'); svg .append("text") .attr("x",pos_Jupiter[0]) .attr("y",pos_Jupiter[1]) .attr("dx", 30) .attr("dy", 8) .attr("font-size", '23px') .attr("font-weight", 300) .attr("font-family", 'Arial') .attr("fill", color_Jupiter) .text('m₅'); svg .append("text") .attr("x",pos_Saturn[0]) .attr("y",pos_Saturn[1]) .attr("dx", 25) .attr("dy", 8) .attr("font-size", '23px') .attr("font-weight", 300) .attr("font-family", 'Arial') .attr("fill", color_Saturn) .text('m₆'); svg .append("text") .attr("x",pos_Sun[0]) .attr("y",pos_Sun[1]) .attr("dx", -45) .attr("dy", -40) .attr("font-size", '25px') .attr("font-weight", 500) .attr("font-family", 'Arial') .style('fill', color_Sun) .text('⊙'); svg .append("text") .attr("x",pos_Jupiter[0]) .attr("y",pos_Jupiter[1]) .attr("dx", -35) .attr("dy", -30) .attr("font-size", '30px') .attr("font-weight", 500) .attr("font-family", 'Arial') .style('fill', color_Jupiter) .text('♃'); svg .append("text") .attr("x",pos_Saturn[0]) .attr("y",pos_Saturn[1]) .attr("dx", -22) .attr("dy", -25) .attr("font-size", '30px') .attr("font-weight", 500) .attr("font-family", 'Arial') .style('fill', color_Saturn) .text('♄'); svg .append("text") .attr("x",pos_Origin[0]/2 + pos_Sun[0]/2) .attr("y",pos_Origin[1]/2 + pos_Sun[1]/2) .attr("dx", -10) .attr("dy", 35) .attr("font-size", '30px') .attr("font-weight", 700) .attr("font-family", 'Arial') .style('fill', color_arrow_origin) .text('R₀'); svg .append("text") .attr("x",pos_Origin[0]/2 + pos_Jupiter[0]/2) .attr("y",pos_Origin[1]/2 + pos_Jupiter[1]/2) .attr("dx", -40) .attr("dy", 10) .attr("font-size", '30px') .attr("font-weight", 700) .attr("font-family", 'Arial') .style('fill', color_arrow_origin) .text('R₅'); svg .append("text") .attr("x",pos_Origin[0]/2 + pos_Saturn[0]/2) .attr("y",pos_Origin[1]/2 + pos_Saturn[1]/2) .attr("dx", -5) .attr("dy", 45) .attr("font-size", '30px') .attr("font-weight", 700) .attr("font-family", 'Arial') .style('fill', color_arrow_origin) .text('R₆'); svg .append("text") .attr("x",pos_Sun[0]/2 + pos_Jupiter[0]/2) .attr("y",pos_Sun[1]/2 + pos_Jupiter[1]/2) .attr("dy", -10) .attr("font-size", '30px') .attr("font-weight", 700) .attr("font-family", 'Arial') .text('r₅'); svg .append("text") .attr("x",pos_Sun[0]/2 + pos_Saturn[0]/2) .attr("y",pos_Sun[1]/2 + pos_Saturn[1]/2) .attr("dy", -15) .attr("font-size", '30px') .attr("font-weight", 700) .attr("font-family", 'Arial') .text('r₆'); svg .append("text") .attr("x",pos_Jupiter[0]/2 + pos_Saturn[0]/2) .attr("y",pos_Jupiter[1]/2 + pos_Saturn[1]/2) .attr("dx", -80) .attr("dy", -5) .attr("font-size", '30px') .attr("font-weight", 700) .attr("font-family", 'Arial') .text('r₆ – r₅'); // Angle psi var arc = d3.arc() .innerRadius(arc_radius) .outerRadius(arc_radius+1) .startAngle(psi_start) .endAngle(psi_end)(); svg.append("path") .attr("transform", `translate(80,170)`) //move from 0,0 coordinates .attr("fill", color_arc) .attr("d", arc); svg .append("text") .attr("x",pos_Sun[0]) .attr("y",pos_Sun[1]) .attr("dx", arc_radius + 10) .attr("dy", 0) .attr("font-size", '22px') .attr("font-weight", 300) .attr("font-style", 'italic') .attr("font-family", 'Arial') .attr("fill", color_arc) .text('ψ'); return svg.node(); }33:T643,Inputs.table(orbital_elements, { columns: [ "index", "name", "relative_mass", "semi_major_axis", "eccentricity", "angular_velocity", "longitude_perihelion", "inclination", "longitude_ascending_node" ], format: { "relative_mass": x => x.toFixed(9), "semi_major_axis": x => x.toFixed(6), "eccentricity": x => x.toFixed(7), "angular_velocity": x => x.toFixed(4), "longitude_perihelion": x => x.toFixed(6), "inclination": x => x.toFixed(5), "longitude_ascending_node": x => x.toFixed(4), }, header: { index: html`Index, i`, name: "Planet", relative_mass: html`m_i / m₀`, semi_major_axis: html`a_i [AU]`, eccentricity: html`e_i`, angular_velocity: html`n_i [°/y]`, longitude_perihelion: html`ϖ_i [°]`, inclination: html`I_i [°]`, longitude_ascending_node: html`Ω_i [°]`, }, align: { index: "left", name: "left", relative_mass: "left", semi_major_axis: "left", eccentricity: "left", angular_velocity: "left", longitude_perihelion: "left", inclination: "left", longitude_ascending_node: "left", }, width: { index: 50, name: 60, relative_mass: 95, semi_major_axis: 80, eccentricity: 85, angular_velocity: 70, longitude_perihelion: 80, inclination: 70, longitude_ascending_node: 70, }, maxWidth: 820, maxHeight: 80, height: 80 })34:[] 2f:{"title":"Three-body problem","cells":[{"id":0,"value":"
Three-body problem
\n
(Perturbation theory in celestial mechanics)
","pinned":false,"hidden":false,"mode":"html"},{"id":1761,"value":"$31","pinned":false,"hidden":false,"mode":"html"},{"id":1973,"value":"three_body = section(\"I.\", \"Three-body system\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":1320,"value":"
Our solar system consists of the Sun (our star, hereafter identified with index 0), and eight major planets: \n
\n
Mercury
\n
Venus
\n
Earth
\n
Mars
\n
Jupiter
\n
Saturn
\n
Uranus
\n
Neptune
\n
\n
\n
\nJupiter (planet No. 5) and Saturn (No. 6) are by far the heaviest planets of the solar system. To study the 3-body problem, it is therefore appropriate to focus on these two planets since they are comparatively much less perturbed by the other planets.\n
","pinned":false,"hidden":false,"mode":"html"},{"id":2049,"value":"
I.1. Evolution equations for the three-body system
\n
Let's consider the simplified 3-body system consisting of the Sun (⊙), Jupiter (♃), and Saturn (♄). We assume their respective masses to be ${tex`m_0 \\gg m_5,m_6`}. With respect to a fixed origin O, the position vectors of the three bodies of masses ${tex`m_0`}, ${tex`m_5`}, and ${tex`m_6`} are ${tex`\\textbf{R₀}`}, ${tex`\\textbf{R₅}`}, and ${tex`\\textbf{R₆}`}, respectively. We denote by ${tex`\\textbf{r₅} = \\textbf{R₅} - \\textbf{R₀}`} and ${tex`\\textbf{r₆} = \\textbf{R₆} - \\textbf{R₀}`} the position vectors of these two planets relative to the Sun.
","pinned":false,"hidden":false,"mode":"html"},{"id":1041,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":2028,"value":"
\n
Figure 1: Model of the simplified 3-body system consisting of the Sun (⊙), Jupiter (♃), and Saturn (♄).
\n
","pinned":false,"hidden":false,"mode":"html"},{"id":1362,"value":"
From Newton's laws of gravitation and motion, the equations of motion of the three masses in the reference frame are
","pinned":false,"hidden":false,"mode":"html"},{"id":1416,"value":"tex.block`\n m_{0} \\textbf{\\\"{R}}_{0}= \\mathcal{G}m_{0}m_{5}\\frac{\\textbf{r}_{5}}{{r_{5}}^3} + \\mathcal{G}m_{0}m_{6}\\frac{\\textbf{r}_{6}}{{r_{6}}^3} \\:,\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1475,"value":"tex.block`\n m_{5} \\textbf{\\\"{R}}_{5}= \\mathcal{G}m_{5}m_{6}\\frac{\\textbf{r}_{6}-\\textbf{r}_{5}} {|r_{6}-r_{5}|^3} - \\mathcal{G}m_{5}m_{0}\\frac{\\textbf{r}_{5}}{{r_{5}}^3} \\:,\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1463,"value":"tex.block`\n m_{6} \\textbf{\\\"{R}}_{6}= \\mathcal{G}m_{6}m_{5}\\frac{\\textbf{r}_{5}-\\textbf{r}_{6}} {|r_{5}-r_{6}|^3} - \\mathcal{G}m_{6}m_{0}\\frac{\\textbf{r}_{6}}{{r_{6}}^3} \\:,\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1480,"value":"
where ${tex`\\mathcal{G}`} is the universal gravitational constant.
\n
The accelarations of the planets are given by
","pinned":false,"hidden":false,"mode":"html"},{"id":1489,"value":"tex.block`\n \\textbf{\\\"{r}}_{5} = \\textbf{\\\"{R}}_{5} - \\textbf{\\\"{R}}_{0} \\:,\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1498,"value":"tex.block`\n \\textbf{\\\"{r}}_{6} = \\textbf{\\\"{R}}_{6} - \\textbf{\\\"{R}}_{0} \\:.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1504,"value":"
Combining the equations given above, it follows that
","pinned":false,"hidden":false,"mode":"html"},{"id":1509,"value":"tex.block`\n \\textbf{\\\"{r}}_{5} + \\mathcal{G}(m_{0}+m_{5})\\frac{\\textbf{r}_{5}}{{r_{5}}^3} = \\mathcal{G}m_{6}\\Big(\\frac{\\textbf{r}_{6}-\\textbf{r}_{5}} {|r_{6}-r_{5}|^3} - \\frac{\\textbf{r}_{6}}{{r_{6}}^3}\\Big) \\:,\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1529,"value":"tex.block`\n \\textbf{\\\"{r}}_{6} + \\mathcal{G}(m_{0}+m_{6})\\frac{\\textbf{r}_{6}}{{r_{6}}^3} = \\mathcal{G}m_{5}\\Big(\\frac{\\textbf{r}_{5}-\\textbf{r}_{6}} {|r_{5}-r_{6}|^3} - \\frac{\\textbf{r}_{5}}{{r_{5}}^3}\\Big) \\:.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1540,"value":"
The right-hand side of these equations can conveniently be expressed as the gradient of potentials (index i = 5, 6):
","pinned":false,"hidden":false,"mode":"html"},{"id":1542,"value":"tex.block`\n \\textbf{\\\"{r}}_{i} + \\mathcal{G}(m_{0}+m_{i})\\frac{\\textbf{r}_{i}}{{r_{i}}^3} = \\vec{\\nabla}_{i}\\mathcal{R}_{i} \\:,\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1566,"value":"
where ${tex`\\vec{\\nabla}_{i}`} is the Del operator (vector differential operator) involving the coordinates with index i; and ${tex`\\mathcal{R}_{i}`} is the disturbing\nfunction, which represents the gravitational potential that arises from the other secondary\nplanet. It can be demonstrated that
","pinned":false,"hidden":false,"mode":"html"},{"id":1582,"value":"tex.block`\n \\mathcal{R}_{5} = \\frac{\\mathcal{G}m_{6}}{|\\textbf{r}_{6}-\\textbf{r}_{5}|} - \\mathcal{G}m_{6} \\frac{\\textbf{r}_{5} \\cdot \\textbf{r}_{6}}{r_{6}^3} \\:,\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1597,"value":"tex.block`\n \\mathcal{R}_{6} = \\frac{\\mathcal{G}m_{5}}{|\\textbf{r}_{6}-\\textbf{r}_{5}|} - \\mathcal{G}m_{5} \\frac{\\textbf{r}_{5} \\cdot \\textbf{r}_{6}}{r_{5}^3} \\:.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1689,"value":"
We denote ${tex`\\psi`} the angle between the position vectors ${tex`\\textbf{r}_{5}`} and ${tex`\\textbf{r}_{6}`}. From the cosine rule, we have
","pinned":false,"hidden":false,"mode":"html"},{"id":1730,"value":"tex.block`\n \\textbf{r}_{5} \\cdot \\textbf{r}_{6} = r_{5}r_{6}\\cos(\\psi) \\:,\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1738,"value":"
and
","pinned":false,"hidden":false,"mode":"html"},{"id":1695,"value":"tex.block`\n |\\textbf{r}_{6}-\\textbf{r}_{5}|^2 = r_{6}^2 +r_{5}^2 -2r_{5}r_{6}\\cos(\\psi) \\:.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1711,"value":"
Hence, we have
","pinned":false,"hidden":false,"mode":"html"},{"id":1701,"value":"tex.block`\n \\begin{array}{rl}\n \\frac{1}{|\\textbf{r}_{6}-\\textbf{r}_{5}|} \n &=& \\Big[r_{6}^2 +r_{5}^2 -2r_{5}r_{6}\\cos(\\psi)\\Big]^{-\\frac{1}{2}} \\\\\n &=& \\frac{1}{r_{6}}\\Big[1 +\\big(\\frac{r_{5}}{r_{6}}\\big)^2 -2\\frac{r_{5}}{r_{6}}\\cos(\\psi)\\Big]^{-\\frac{1}{2}} \\:.\n \\end{array}\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1810,"value":"
Finally, introducing a dimensionless variable ${tex`\\alpha`}, we obtain
","pinned":false,"hidden":false,"mode":"html"},{"id":1796,"value":"tex.block`\n \\frac{1}{|\\textbf{r}_{6}-\\textbf{r}_{5}|} = \\frac{1}{r_{6}}\\Big[1 +\\alpha^2 -2\\alpha\\cos(\\psi)\\Big]^{-\\frac{1}{2}} \\:, \n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1965,"value":"
with ${tex`\\alpha = \\frac{r_{5}}{r_{6}} < 1`}.
","pinned":false,"hidden":false,"mode":"html"},{"id":1742,"value":"disturbing_function = md `
It would actually require many more lines of calculations to bring the disturbing function ${tex`\\mathcal{R}_{i}`} into a usable form. Such detailed developments can be found in the handbooks listed in the [References](#references) (e.g, Murray & Dermott, Chapter 6).^{[[1](#Murray_Dermott),[2](#Fitzpatrick)]} Here, we will not push the matter any further because we have already reached an essential point: Computing the perturbation functions involves using the Laplace coefficients that are briefly presented in [Appendix A](#Laplace_coefficients).
`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2158,"value":" md `
I.2. Approximation of the disturbing function
The disturbing function can be expanded into an infinite series whose terms can be classified either as secular, resonant, or short period, depending on the physical problem under consideration.^{[[1](#Murray_Dermott)]} Here, we will investigate the influence of interplanetary gravitational interactions on planetary orbits when integrated over long periods of time (i.e., secular perturbations).\n
`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2185,"value":"
For this purpose, it is is convenient to introduce alternative orbital\nelements for the eccentricity and inclination:
","pinned":false,"hidden":false,"mode":"html"},{"id":2187,"value":"tex.block`\n h_{i} = e_{i} \\sin(\\varpi_{i}) \\:, \\quad k_{i} = e_{i} \\cos(\\varpi_{i}) \\:, \\\\\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2195,"value":"tex.block`\n p_{i} = I_{i} \\sin(\\Omega_{i}) \\:, \\quad q_{i} = I_{i} \\cos(\\Omega_{i}) \\:. \\\\\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2199,"value":"
These variables relate the eccentricity ${tex`e_{i}`} to the longitude of the perihelion ${tex`\\varpi_{i}`}, and the inclination ${tex`I_{i}`} to the longitude of the ascending node ${tex`\\Omega_{i}`} for each planet ${tex`i`} (where ${tex`i = 5, 6`}, for Jupiter and Saturn, respectively).
","pinned":false,"hidden":false,"mode":"html"},{"id":2213,"value":"
It is also helpful to define the following parameters:
","pinned":false,"hidden":false,"mode":"html"},{"id":2215,"value":"tex.block`\n \\alpha_{ij} = \n\\begin{cases}\n \\frac{a_{i}}{a_j} \\quad &\\textrm{if} \\quad a_{i}a_{j} \\; (i>j) \\\\\n\\end{cases} \\:, \n\\\\\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2235,"value":"
and
","pinned":false,"hidden":false,"mode":"html"},{"id":2237,"value":"tex.block`\n \\bar{\\alpha}_{ij} = \n\\begin{cases}\n \\frac{a_{i}}{a_j} \\quad &\\textrm{if} \\quad a_{i}a_{j} \\; (i>j) \\\\\n\\end{cases} \\:, \n\\\\\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2243,"value":"
where ${tex`a_{i}`} and ${tex`a_{j}`} are the semi-major axes of planets ${tex`i`} and ${tex`j`}.
","pinned":false,"hidden":false,"mode":"html"},{"id":2247,"value":"
The general secular part of the disturbing function can be written:
","pinned":false,"hidden":false,"mode":"html"},{"id":2250,"value":"tex.block`\n \\mathcal{R}_{i} = n_{i}a_{i}^2 \\Big[\\frac{1}{2}A_{ii}(h_{i}^2 + k_{i}^2)+ A_{ij}(h_{i}h_{j} + k_{i}k_{j}) \\\\\n \\qquad \\qquad \\qquad + \\frac{1}{2}B_{ii}(p_{i}^2 + q_{i}^2)+ B_{ij}(p_{i}p_{j} + q_{i}q_{j} \\Big] \\:,\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2269,"value":"
where ${tex`n_{i}`} is the orbital angular velocity of planet ${tex`i`}. The remaining quantities are constant elements from two matrices ${tex`\\textbf{A}`} and ${tex`\\textbf{B}`}:
","pinned":false,"hidden":false,"mode":"html"},{"id":2275,"value":"tex.block`\n A_{ii} = \\sum_{j, \\: j \\neq i} \\frac{n_i}{4}\\frac{m_{j}}{m_{0}+m_{i}}\\alpha_{ij}\\bar{\\alpha}_{ij} b_{3/2}^{(1)}(\\alpha_{ij}) \\:,\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2293,"value":"tex.block`\n A_{ij} = -\\frac{n_i}{4}\\frac{m_{j}}{m_{0}+m_{i}}\\alpha_{ij}\\bar{\\alpha}_{ij} b_{3/2}^{(2)}(\\alpha_{ij}) \\:,\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2298,"value":"tex.block`\n B_{ii} = -\\sum_{j, \\: j \\neq i} \\frac{n_i}{4}\\frac{m_{j}}{m_{0}+m_{i}}\\alpha_{ij}\\bar{\\alpha}_{ij} b_{3/2}^{(1)}(\\alpha_{ij}) \\:,\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2305,"value":"tex.block`\n B_{ij} = \\frac{n_i}{4}\\frac{m_{j}}{m_{0}+m_{i}}\\alpha_{ij}\\bar{\\alpha}_{ij} b_{3/2}^{(1)}(\\alpha_{ij}) \\:.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2311,"value":"md `
In these equations, one can recognize the Laplace coefficients ${tex`b_s^{(m)}(\\alpha)`} discussed in [Appendix A](#Laplace_coefficients).\nNote that these expressions are generalised forms applicable for the N-body problem (i.e., for ${tex`i,j = 1 ... 8`}, ${tex`i \\neq j`}).
`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2335,"value":" md `
I.3. Eigenfrequencies and eigenvectors
The solutions of the equations of motions are
`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2341,"value":"tex.block`\n h_{i} = \\sum_{j} e_{ij} \\sin(g_{j}t + \\beta_{j}) \\:, \\quad k_{i} = \\sum_{j} e_{ij} \\cos(g_{j}t + \\beta_{j}) \\:,\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2361,"value":"tex.block`\n p_{i} = \\sum_{j} I_{ij} \\sin(f_{j}t + \\gamma_{j}) \\:, \\quad q_{i} = \\sum_{j} I_{ij} \\cos(f_{j}t + \\gamma_{j}) \\:,\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2346,"value":"
where ${tex`t \\;`} is the time variable, \n
${tex`g_{j} \\;`} are the eigenvalues of matrix ${tex`\\textbf{A}`},\n
${tex`e_{ij} \\;`} are the eigenvectors of matrix ${tex`\\textbf{A}`},\n
${tex`f_{j} \\;`} are the eigenvalues of matrix ${tex`\\textbf{B}`},\n
${tex`I_{ij} \\;`} are the eigenvectors of matrix ${tex`\\textbf{B}`},\n
and ${tex`\\beta_{j}`} and ${tex`\\gamma_{j}`} are the phase angles.\n
","pinned":false,"hidden":false,"mode":"html"},{"id":2380,"value":"
The determination of the secular evolution of planetary orbits is now reduced to a numerical problem consisting in finding the eigenvalues and eigenvectors of matrices ${tex`\\textbf{A}`} and ${tex`\\textbf{B}`}. The final solution can be found by solving these equations for the initial condition (${tex`t=0`}) using the observation values for each planet ${tex`i`}: ${tex`e_{i}`}, ${tex`\\varpi_{i}`}, ${tex`I_{i}`}, and ${tex`\\Omega_{i}`}.
","pinned":false,"hidden":false,"mode":"html"},{"id":825,"value":"
","pinned":false,"hidden":false,"mode":"html"},{"id":1978,"value":"secular_perturbations = section(\"II.\", \"Secular perturbations for two planets\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":2396,"value":"md`
As already suggested in the begining of this notebook, we will limit the study to the 3-body system. \nWe will therefore apply the theory given above to the case of Jupiter and Saturn orbiting the Sun. \nFor comparison purposes, the parameters used below are those given by Murray & Dermott.^{[[1](#Murray_Dermott)]}\n

In 1983 the orbital parameters of Jupiter (planet No. 5) and Saturn (No. 6) were as follows:
`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2420,"value":"$33","pinned":false,"hidden":false,"mode":"ojs"},{"id":3507,"value":"
In section II.1, we complete all the calculations needed to solve the equations of motion.\nYou can jump directly to [section II.2](#graphs) to see some of the insightful graphs that can be ploted from the Laplace-Lagrande perturbation theory.
","pinned":false,"hidden":false,"mode":"md"},{"id":2137,"value":"
II.1. Jupiter and Saturn: Detailed calculations
\n","pinned":false,"hidden":false,"mode":"html"},{"id":2767,"value":"md `
The ratio of their semi-major axes is ${tex`\\alpha = \\alpha_{56} = ${alpha_56.toFixed(6)}`}. From the definition of Laplace coefficients given in [Appendix A](#Laplace_coefficients), we obtain:
`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2763,"value":"tex.block`\n b_{3/2}^{(1)}(\\alpha_{56})= ${(b_m_s_alpha(1,1.5,alpha_56)[0]).toFixed(5)} \\:, \\qquad \n b_{3/2}^{(2)}(\\alpha_{56})= ${(b_m_s_alpha(2,1.5,alpha_56)[0]).toFixed(5)} \\:. \\\\\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2778,"value":"
The elements of matrices ${tex`\\textbf{A}`} and ${tex`\\textbf{B}`} are:
","pinned":false,"hidden":false,"mode":"html"},{"id":2782,"value":"tex.block`\n\\textbf{A} = \\begin{bmatrix}\n A_{55} & A_{56} \\\\\n A_{65} & A_{66} \n\\end{bmatrix} =\n\\begin{bmatrix}\n ${A_55.toFixed(8)} & ${A_56.toFixed(8)} \\\\\n ${A_65.toFixed(8)} & ${A_66.toFixed(8)}\n\\end{bmatrix}\n\\,, \\qquad \n \\textbf{B} = \\begin{bmatrix}\n B_{55} & B_{56} \\\\\n B_{65} & B_{66}\n\\end{bmatrix} =\n\\begin{bmatrix}\n ${B_55.toFixed(8)} & ${B_56.toFixed(8)} \\\\\n ${B_65.toFixed(8)} & ${B_66.toFixed(8)}\n\\end{bmatrix}.\n\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2877,"value":"
The eigenvalues of these matrix are:
","pinned":false,"hidden":false,"mode":"html"},{"id":2886,"value":"tex.block`\n g_{5} = ${(g_5*1E4).toFixed(5)} \\times 10^{-4} \\; \\textrm{°/year} = ${(g_5*3600).toFixed(3)} \\; \\textrm{''/year} \\:, \\\\\n g_{6} = ${(g_6*1E3).toFixed(5)} \\times 10^{-3} \\; \\textrm{°/year} = ${(g_6*3600).toFixed(3)} \\; \\textrm{''/year} \\:, \\\\\n f_{5} = ${f_5.toFixed(0)} \\textrm{°/year} = ${(f_5*3600).toFixed(0)} \\; \\textrm{''/year} \\:, \\\\\n f_{6} = ${(f_6*1E3).toFixed(5)} \\times 10^{-3} \\; \\textrm{°/year} = ${(f_6*3600).toFixed(3)} \\; \\textrm{''/year} \\:. \n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2928,"value":"
Note that the javascript library Math.js is used to compute these eigenvalues. Using that library, we also compute the four associated unscaled eigenvectors with components ${tex`\\tilde{e}_{ij}`} and ${tex`\\tilde{I}_{ij}`}:
","pinned":false,"hidden":false,"mode":"html"},{"id":2982,"value":"tex.block`\n\\begin{pmatrix}\n \\tilde{e}_{55} \\\\\n \\tilde{e}_{65}\n\\end{pmatrix} = \n\\begin{pmatrix}\n ${(eigs_A[1][0][0]).toFixed(8)} \\\\\n ${(eigs_A[1][1][0]).toFixed(8)}\n\\end{pmatrix} , \\qquad\n\n\\begin{pmatrix}\n \\tilde{e}_{56} \\\\\n \\tilde{e}_{66}\n\\end{pmatrix} = \n\\begin{pmatrix}\n ${(eigs_A[1][0][1]).toFixed(8)} \\\\\n ${(eigs_A[1][1][1]).toFixed(8)}\n\\end{pmatrix} , \\\\\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":3007,"value":"tex.block`\n\\begin{pmatrix}\n \\tilde{I}_{55} \\\\\n \\tilde{I}_{65}\n\\end{pmatrix} = \n\\begin{pmatrix}\n ${(eigs_B[1][0][0]).toFixed(8)} \\\\\n ${(eigs_B[1][0][0]).toFixed(8)}\n\\end{pmatrix} , \\qquad\n\n\\begin{pmatrix}\n \\tilde{I}_{56} \\\\\n \\tilde{I}_{66}\n\\end{pmatrix} = \n\\begin{pmatrix}\n ${(eigs_B[1][0][1]).toFixed(8)} \\\\\n ${(eigs_B[1][1][1]).toFixed(8)}\n\\end{pmatrix} . \\\\\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":3158,"value":"
The scaled eigenvectors are:
","pinned":false,"hidden":false,"mode":"html"},{"id":3160,"value":"tex.block`\n\\begin{pmatrix}\n e_{55} \\\\\n e_{65}\n\\end{pmatrix} = S_{5}\n\\begin{pmatrix}\n \\tilde{e}_{55} \\\\\n \\tilde{e}_{65}\n\\end{pmatrix} , \\qquad\n\n\\begin{pmatrix}\n e_{56} \\\\\n e_{66}\n\\end{pmatrix} = S_{6}\n\\begin{pmatrix}\n \\tilde{e}_{56} \\\\\n \\tilde{e}_{66}\n\\end{pmatrix} , \\\\\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":3165,"value":"tex.block`\n\\begin{pmatrix}\n I_{55} \\\\\n I_{65}\n\\end{pmatrix} = T_{5}\n\\begin{pmatrix}\n \\tilde{I}_{55} \\\\\n \\tilde{I}_{65}\n\\end{pmatrix} , \\qquad\n\n\\begin{pmatrix}\n I_{56} \\\\\n I_{66}\n\\end{pmatrix} = T_{6}\n\\begin{pmatrix}\n \\tilde{I}_{56} \\\\\n \\tilde{I}_{66}\n\\end{pmatrix} , \\\\\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":3037,"value":"
where ${tex`S_{5}`}, ${tex`S_{6}`}, ${tex`T_{5}`}, and ${tex`T_{6}`} are scaling factors that have to be determined.
\n
Appropriate scaling can be found from the boundary conditions. At time ${tex`t=0`}, we have
","pinned":false,"hidden":false,"mode":"html"},{"id":3041,"value":"tex.block`\n h_{5} = e_{5} \\sin(\\varpi_{5}) = ${h_5.toFixed(7)} \\,, \\quad h_{6} = e_{6} \\sin(\\varpi_{6}) = ${h_6.toFixed(7)} \\,, \\\\\n k_{5} = e_{5} \\cos(\\varpi_{5}) = ${k_5.toFixed(7)} \\,, \\quad k_{6} = e_{6} \\cos(\\varpi_{6}) = ${k_6.toFixed(8)} \\,, \\\\\n \n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":3082,"value":"
and
","pinned":false,"hidden":false,"mode":"html"},{"id":3086,"value":"tex.block`\n p_{5} = I_{5} \\sin(\\Omega_{5}) = ${p_5.toFixed(7)}~\\textrm{rad} \\,, \\quad p_{6} = I_{6} \\sin(\\Omega_{6}) = ${p_6.toFixed(7)}~\\textrm{rad} \\,,\\\\\n q_{5} = I_{5} \\cos(\\Omega_{5}) = ${q_5.toFixed(8)}~\\textrm{rad} \\,, \\quad q_{6} = I_{6} \\cos(\\Omega_{6}) = ${q_6.toFixed(7)}~\\textrm{rad}. \\\\\n \n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":3136,"value":"
The two sets of equations that have to be solved are
","pinned":false,"hidden":false,"mode":"html"},{"id":3138,"value":"tex.block`\n \\begin{cases}\n h_{5} = e_{55}\\sin(\\beta_{5}) + e_{56}\\sin(\\beta_{6}) = \\tilde{e}_{55} \\cdot S_{5}\\sin(\\beta_{5}) + \\tilde{e}_{56} \\cdot S_{6}\\sin(\\beta_{6}) \\\\\n h_{6} = e_{65}\\sin(\\beta_{5}) + e_{66}\\sin(\\beta_{6}) = \\tilde{e}_{65} \\cdot S_{5}\\sin(\\beta_{5}) + \\tilde{e}_{66} \\cdot S_{6}\\sin(\\beta_{6}) \\\\\n k_{5} = e_{55}\\cos(\\beta_{5}) + e_{56}\\cos(\\beta_{6}) = \\tilde{e}_{55} \\cdot S_{5}\\cos(\\beta_{5}) + \\tilde{e}_{56} \\cdot S_{6}\\cos(\\beta_{6})\\\\\n k_{6} = e_{65}\\cos(\\beta_{5}) + e_{66}\\cos(\\beta_{6}) = \\tilde{e}_{65} \\cdot S_{5}\\cos(\\beta_{5}) + \\tilde{e}_{66} \\cdot S_{6}\\cos(\\beta_{6}) \\\\\n\\end{cases} \\:,\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":3152,"value":"
and
","pinned":false,"hidden":false,"mode":"html"},{"id":3154,"value":"tex.block`\n \\begin{cases}\n p_{5} = I_{55}\\sin(\\gamma_{5}) + e_{56}\\sin(\\gamma_{6}) = \\tilde{I}_{55} \\cdot T_{5}\\sin(\\gamma_{5}) + \\tilde{I}_{56} \\cdot T_{6}\\sin(\\gamma_{6}) \\\\\n p_{6} = I_{65}\\sin(\\gamma_{5}) + e_{66}\\sin(\\gamma_{6}) = \\tilde{I}_{65} \\cdot T_{5}\\sin(\\gamma_{5}) + \\tilde{I}_{66} \\cdot T_{6}\\sin(\\gamma_{6}) \\\\\n q_{5} = I_{55}\\cos(\\gamma_{5}) + e_{56}\\cos(\\gamma_{6}) = \\tilde{I}_{55} \\cdot T_{5}\\cos(\\gamma_{5}) + \\tilde{I}_{56} \\cdot T_{6}\\cos(\\gamma_{6})\\\\\n q_{6} = I_{65}\\cos(\\gamma_{5}) + e_{66}\\cos(\\gamma_{6}) = \\tilde{I}_{65} \\cdot T_{5}\\cos(\\gamma_{5}) + \\tilde{I}_{66} \\cdot T_{6}\\cos(\\gamma_{6}) \\\\\n\\end{cases} \\:.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":3170,"value":"
Let's solve the first set of equations. In matrix representation, these equations are written as follows:
","pinned":false,"hidden":false,"mode":"html"},{"id":3174,"value":"tex.block`\n\\textbf{y} = \\begin{bmatrix}\n h_{5} \\\\\n h_{6} \\\\\n k_{5} \\\\\n k_{6}\n\\end{bmatrix} =\n\\begin{bmatrix}\n \\tilde{e}_{55} & \\tilde{e}_{56} & 0 & 0\\\\\n \\tilde{e}_{65} & \\tilde{e}_{66} & 0 & 0\\\\\n 0 & 0 & \\tilde{e}_{55} & \\tilde{e}_{56} \\\\\n 0 & 0 & \\tilde{e}_{65} & \\tilde{e}_{66}\n\\end{bmatrix} \n\\begin{bmatrix}\n S_{5}\\sin(\\beta_{5}) \\\\\n S_{6}\\sin(\\beta_{6}) \\\\\n S_{5}\\cos(\\beta_{5}) \\\\\n S_{6}\\cos(\\beta_{6})\n\\end{bmatrix} = \\textbf{M} \\cdot \\textbf{x}\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":3191,"value":"
The solution is found by solving
","pinned":false,"hidden":false,"mode":"html"},{"id":3193,"value":"tex.block`\n\\textbf{x} = \\big(\\textbf{M}^{\\textrm{T}}\\textbf{M}\\big)^{-1} \\textbf{M}^{\\textrm{T}} \\cdot \\textbf{y}\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":3240,"value":"
Using the library Math.js, we find:
","pinned":false,"hidden":false,"mode":"html"},{"id":3244,"value":"tex.block`\n\\textbf{x} = \\begin{bmatrix}\n S_{5}\\sin(\\beta_{5}) \\\\\n S_{6}\\sin(\\beta_{6}) \\\\\n S_{5}\\cos(\\beta_{5}) \\\\\n S_{6}\\cos(\\beta_{6})\n\\end{bmatrix} = \n\\begin{bmatrix}\n ${x[0].toFixed(8)} \\\\\n ${x[1].toFixed(8)} \\\\\n ${x[2].toFixed(8)} \\\\\n ${x[3].toFixed(8)}\n\\end{bmatrix} \n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":3253,"value":"
Thus, we deduce:
","pinned":false,"hidden":false,"mode":"html"},{"id":3257,"value":"tex.block`\n S_{5} = ${S_5.toFixed(6)} \\:, \\quad S_{6} = ${S_6.toFixed(6)} \\:,\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":3283,"value":"tex.block`\n \\beta_{5} = ${beta_5.toFixed(3)}° \\,, \\quad \\beta_{6} = ${beta_6.toFixed(4)}° \\,,\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":3271,"value":"
and
","pinned":false,"hidden":false,"mode":"html"},{"id":3273,"value":"tex.block`\n\\begin{pmatrix}\n e_{55} \\\\\n e_{65}\n\\end{pmatrix} = \n\\begin{pmatrix}\n ${e_55.toFixed(7)} \\\\\n ${e_65.toFixed(7)}\n\\end{pmatrix} , \\qquad\n\n\\begin{pmatrix}\n e_{56} \\\\\n e_{66}\n\\end{pmatrix} = \n\\begin{pmatrix}\n ${e_56.toFixed(7)} \\\\\n ${e_66.toFixed(6)}\n\\end{pmatrix} . \\\\\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":3298,"value":"
Following the same procedure to solve the second set of equations, we find:
","pinned":false,"hidden":false,"mode":"html"},{"id":3323,"value":"tex.block`\n T_{5} = ${T_5.toFixed(6)} \\:, \\quad T_{6} = ${T_6.toFixed(6)} \\:,\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":3327,"value":"tex.block`\n \\gamma_{5} = ${gamma_5.toFixed(2)}° \\,, \\quad \\gamma_{6} = ${gamma_6.toFixed(4)}° \\,,\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":3332,"value":"
and
","pinned":false,"hidden":false,"mode":"html"},{"id":3335,"value":"tex.block`\n\\begin{pmatrix}\n I_{55} \\\\\n I_{65}\n\\end{pmatrix} = \n\\begin{pmatrix}\n ${I_55.toFixed(7)}~\\textrm{rad}\\\\\n ${I_65.toFixed(7)}~\\textrm{rad}\n\\end{pmatrix} , \\qquad\n\n\\begin{pmatrix}\n I_{56} \\\\\n I_{66}\n\\end{pmatrix} = \n\\begin{pmatrix}\n ${I_56.toFixed(7)}~\\textrm{rad} \\\\\n ${I_66.toFixed(6)}~\\textrm{rad}\n\\end{pmatrix} . \\\\\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":3411,"value":"
We have now determined all the constants that intervene in the equations of motion of Jupiter and Saturn.
","pinned":false,"hidden":false,"mode":"html"},{"id":823,"value":"graphs = md`
II.2. Secular evolution of the orbits of Jupiter and Saturn
`","pinned":false,"hidden":false,"mode":"ojs"},{"id":3549,"value":"
In section II.1, we have obtained ${tex`h_{i}(t)`}, ${tex`k_{i}(t)`}, ${tex`p_{i}(t)`}, and ${tex`q_{i}(t)`} for Jupiter (${tex`i=5`}) and Saturn (${tex`i=6`}). We can therefore derive the orbital elements of the two planets at any time ${tex`t`}. For example, for the eccentricites we have
","pinned":false,"hidden":false,"mode":"html"},{"id":3560,"value":"tex.block`\n e_{i}(t) = \\sqrt{h_{i}^2(t)+k_{i}^2(t)}\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":3571,"value":"
Developing this equation, we find
","pinned":false,"hidden":false,"mode":"html"},{"id":3567,"value":"tex.block`\n e_{5}(t) = \\sqrt{e_{55}^2 + e_{56}^2 +2e_{55}e_{56}\\cos\\big((g_{6}-g_{5})t+(\\beta_{6}-\\beta_{5})\\big)}\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":3572,"value":"tex.block`\n e_{6}(t) = \\sqrt{e_{66}^2 + e_{65}^2 +2e_{66}e_{65}\\cos\\big((g_{6}-g_{5})t+(\\beta_{6}-\\beta_{5})\\big)}\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":3566,"value":"
The numerical application gives
","pinned":false,"hidden":false,"mode":"html"},{"id":3579,"value":"tex.block`\n e_{5}(t) = \\sqrt{${(math.pow(e_55,2) \n + math.pow(e_56,2)).toFixed(5)} ${math.sign(e_55*e_56)<0?'-':'+'} ${(math.abs(2*e_55*e_56)).toFixed(5)}\\cos\\big(${(g_6-g_5).toFixed(5)}\\cdot t +${(beta_6-beta_5).toFixed(2)}°\\big)} \\:,\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":3600,"value":"tex.block`\n e_{6}(t) = \\sqrt{${(math.pow(e_66,2) \n + math.pow(e_65,2)).toFixed(5)} ${math.sign(e_66*e_65)<0?'-':'+'} ${(math.abs(2*e_66*e_65)).toFixed(5)}\\cos\\big(${(g_6-g_5).toFixed(5)}\\cdot t +${(beta_6-beta_5).toFixed(2)}°\\big)} \\:,\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":3621,"value":"
where phase angles are in degrees and angular velocities are in degrees per year.
","pinned":false,"hidden":false,"mode":"html"},{"id":3626,"value":"
Similarely, we find for the inclinations (expressed in radians):
","pinned":false,"hidden":false,"mode":"html"},{"id":3630,"value":"tex.block`\n I_{5}(t) = \\sqrt{${(math.pow(I_55,2) \n + math.pow(I_56,2)).toFixed(6)} ${math.sign(I_55*I_56)<0?'-':'+'} ${(math.abs(2*I_55*I_56)).toFixed(6)}\\cos\\big(${(f_5-f_6).toFixed(5)}\\cdot t +${(gamma_5-gamma_6).toFixed(2)}°\\big)} \\:,\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":3643,"value":"tex.block`\n I_{6}(t) = \\sqrt{${(math.pow(I_66,2) \n + math.pow(I_65,2)).toFixed(6)} ${math.sign(I_66*I_65)<0?'-':'+'} ${(math.abs(2*I_66*I_65)).toFixed(6)}\\cos\\big(${(f_5-f_6).toFixed(5)}\\cdot t +${(gamma_5-gamma_6).toFixed(2)}°\\big)} \\:,\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":3638,"value":"
where, inside the cosine functions, phase angles are in degrees and angular velocities are in degrees per year.
","pinned":false,"hidden":false,"mode":"html"},{"id":4064,"value":"
The most striking feature is the strong coupling of the two planets. The eccentricities of Jupiter and Saturn vary with a period of around ${100*(360/math.abs(g_6-g_5)/100).toFixed(0)} years. The secular evolution of their inclination shows a periodicity of around ${100*(360/math.abs(f_6-f_5)/100).toFixed(0)} years.
\n
The secular variations of the orbits of Jupiter and Saturn are shown in Figure 2. Note that the solutions that we have derived from the\nthe orbital perturbation theory is only an approximation to the actual variations. Indeed, the perturbations from the planets Uranus and Neptune also exert a noticeable influence on the orbits of Jupiter and Saturn.
\n
","pinned":false,"hidden":false,"mode":"html"},{"id":3768,"value":"figure2 = Plot.plot({\n height: 300,\n width: 550,\n x: {\n grid: true,\n label: \"→ Years \"\n }, \n y: {\n grid: false,\n label: \"↑ Eccentricity\"\n },\n style: { fontSize: \"11px\" },\n marks: [\n Plot.ruleX([-100000,100000]),\n Plot.ruleX([0], {strokeDasharray: \"4,3\"}),\n Plot.ruleY([0,0.1]), \n Plot.line(plot_e5, Plot.windowY({k: 1, x: \"time\", y: \"eccentricity\", dash: '8,2', strokeWidth:2.5, stroke: color_Jupiter2})),\n Plot.line(plot_e6, Plot.windowY({k: 1, x: \"time\", y: \"eccentricity\", strokeWidth:2.5, stroke: color_Saturn2})),\n\n Plot.text(plot_e5, Plot.selectLast({\n x: \"time\",\n y: \"eccentricity\",\n text: () => \"Jupiter\",\n fontSize: 13.5,\n fontWeight: 500,\n fill: color_Jupiter2,\n dx: -34,\n dy: -34\n })),\n \n Plot.text(plot_e6, Plot.selectMaxY({\n x: \"time\",\n y: \"eccentricity\",\n text: () => \"Saturn\",\n fontSize: 13.5,\n fontWeight: 500,\n fill: color_Saturn2,\n dy: -15\n })),\n \n ]\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":3913,"value":"Plot.plot({\n height: 300,\n width: 550,\n x: {\n grid: true,\n label: \"→ Years \"\n }, \n y: {\n grid: false,\n label: \"↑ Inclination (°)\"\n },\n style: { fontSize: \"11px\" },\n marks: [\n Plot.ruleX([-100000,100000]),\n Plot.ruleX([0], {strokeDasharray: \"4,3\"}),\n Plot.ruleY([0,3]), \n Plot.line(plot_I5, Plot.windowY({k: 1, x: \"time\", y: \"inclination\", strokeWidth:2.5, stroke: color_Jupiter2})),\n Plot.line(plot_I6, Plot.windowY({k: 1, x: \"time\", y: \"inclination\", strokeWidth:2.5, stroke: color_Saturn2})),\n\n Plot.text(plot_I5, Plot.selectMaxY({\n x: \"time\",\n y: \"inclination\",\n text: () => \"Jupiter\",\n fontSize: 13.5,\n fontWeight: 500,\n fill: color_Jupiter2,\n dx: 0,\n dy: -15\n })),\n \n Plot.text(plot_I6, Plot.selectLast({\n x: \"time\",\n y: \"inclination\",\n text: () => \"Saturn\",\n fontSize: 13.5,\n fontWeight: 500,\n fill: color_Saturn2,\n dx: -114,\n dy: -24\n })),\n \n ]\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":4006,"value":"
\n
Figure 2: Eccentricities and inclinations of Jupiter and Saturn derived from\n the 3-body secular perturbation theory. The time span is 200'000 years, centered on year 1983.
\n
","pinned":false,"hidden":false,"mode":"html"},{"id":3467,"value":"
","pinned":false,"hidden":false,"mode":"html"},{"id":1874,"value":"
","pinned":false,"hidden":false,"mode":"html"},{"id":1982,"value":"references = section(\"\", \"References\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":626,"value":"Murray_Dermott = reference(\"disturbing_function\", 1, \"C. D. Murray & S. F. Dermott, \\\"Solar System Dynamics,\\\" Cambridge University Press, 2000. DOI: 10.1017/CBO9781139174817\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":636,"value":"Fitzpatrick = reference(\"disturbing_function\", 2, \"R. Fitzpatrick, \\\"An Introduction to Celestial Mechanics,\\\" Cambridge University Press, 2012. DOI: 10.1017/CBO9781139152310\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":152,"value":"Izsak1963 = reference(\"Laplace_series_expansion\", 3, \"I. G. Izsak, \\\"Laplace coefficients and their Newcomb derivatives,\\\" SAO Special Report 129, 1963. [PDF]\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":183,"value":"Mardling2013 = reference(\"Laplace_series_expansion\", 4, \"R. A. Mardling, \\\"New developments for modern celestial mechanics — I. General coplanar three-body systems. Application to exoplanets,\\\" Monthly Notices of the Royal Astronomical Society 435 (3), pp. 2187–2226, November 2013. DOI: 10.1093/mnras/stt1438\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":1898,"value":"
","pinned":false,"hidden":false,"mode":"html"},{"id":1841,"value":"Laplace_coefficients = section(\"Appendix A:\", \"Laplace coefficients\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":764,"value":"md `
As discussed in [Section I](#disturbing_function), Laplace coefficients are key functions that need to be manipulated intensively when investigating the orbit of planets in the framework of the perturbation theory of planetary motion.^{[[1](#Murray_Dermott),[2](#Fitzpatrick)]} These coefficients look as follows:
`","pinned":false,"hidden":false,"mode":"ojs"},{"id":12,"value":"tex.block`\nb_s^{(m)}(\\alpha) = \\frac{1}{\\pi} \\int_{0}^{2\\pi} \\frac{\\cos(m\\psi)}{ \\big( 1 - 2\\alpha\\cos(\\psi) + \\alpha^2 \\big)^s} \\,d\\psi \\:,\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":37,"value":"where ${tex`m`} is a natural number, \n
${tex`s`} is a half-integer,\n
and ${tex`\\alpha`} is a positive real number < 1.","pinned":false,"hidden":false,"mode":"html"},{"id":109,"value":"Laplace_series_expansion = md `
The series expansion of ${tex`b_s^{(m)}(\\alpha)`} is ^{[[3](#Izsak1963),[4](#Mardling2013)]}:
`","pinned":false,"hidden":false,"mode":"ojs"},{"id":116,"value":"tex.block`\nb_s^{(m)}(\\alpha) = 2\\alpha^m \\sum_{p=0}^{\\infin} \\; \\frac{\\Gamma(s+p)\\;\\Gamma(s+m+p)}{[\\Gamma(s)]^2 \\; p! \\; (m+p)!} \\; \\alpha^{2p} \\:.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1906,"value":"
Further details on Laplace coefficients can be found in this other Observable notebook.
","pinned":false,"hidden":false,"mode":"html"},{"id":577,"value":"
Calculation of the Laplace coefficient
\n
To calculate a Laplace coefficient, click on the sliders or enter new values in the input fields.
","pinned":false,"hidden":false,"mode":"html"},{"id":543,"value":"viewof m = Inputs.range([0, 10], {value: 1, step: 1, label: \"m = \"})","pinned":false,"hidden":false,"mode":"ojs"},{"id":315,"value":"viewof n = Inputs.range([-10, 10], {value: 1, step: 1, label: \"s = n + ½, with n = \"})","pinned":false,"hidden":false,"mode":"ojs"},{"id":302,"value":"viewof alpha = Inputs.range([0.00000001, 0.99999999], {value: 5.202545/9.554841, step: 0.00000001, label: \"α = \"})","pinned":false,"hidden":false,"mode":"ojs"},{"id":590,"value":"viewof p_max = Inputs.range([20, 200], {value: 200, step: 1, label: \"max {p} = \"})","pinned":false,"hidden":false,"mode":"ojs"},{"id":300,"value":"tex.block`\n\n\\begin{array}{rl} \n\\\\\n & m = ${m} \\:, \\quad s = ${math.sign(2*n+1)=='-1'? '-':''}\\frac{${math.abs(2*n+1)}}{2} \\:, \\quad \\alpha = ${alpha} \\:,\\\\ \\\\\n b_{${2*n+1}/2}^{(${m})}(${alpha}) &\n = \\quad ${b_m_s_alpha(m,n+0.5,alpha,p_max)[0]} \\qquad (p = 0 \\,... \\, ${b_m_s_alpha(m,n+0.5,alpha,p_max)[1]}) \\:.\n\\end{array}\n\n \n\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":573,"value":"
","pinned":false,"hidden":false,"mode":"html"},{"id":10,"value":"
","pinned":false,"hidden":false,"mode":"html"},{"id":1914,"value":"other_appendix = section(\"Appendix B:\", \"Javascript / code + variables\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":873,"value":"math = require(\"mathjs@11.0.0\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":197,"value":"function b_m_s_alpha (m,s,alpha, P_max = 100){ // P_max is an optional parameter\n let b_m_s_alpha = 0;\n let sum = 0;\n let sum_test = 0;\n let p = 0;\n for(p; p <= P_max && Number.isFinite(sum_test); p++){\n sum_test = math.gamma(s+p)/math.pow(math.gamma(s),2); \n sum_test /= math.factorial(p);\n sum_test *= math.gamma(s+m+p)*math.pow(alpha,2*p);\n sum_test /= math.factorial(m+p);\n Number.isFinite(sum_test)? sum += sum_test : null; // Check if the calculated value is a finite number\n }\n b_m_s_alpha = 2*math.pow(alpha,m)*sum;\n return [b_m_s_alpha,p-1];\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1838,"value":"section = (section, text) => {\n return html`
${section} ${text}
`;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":147,"value":"// https://observablehq.com/@jashkenas/markdown-style-footnotes\nreference = (referencingCellName, number, text) => {\n return html`
↫ ${text}
`;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":2682,"value":"function planet_data (index) {\n return orbital_elements.filter(d => d.index==index)[0];\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":2413,"value":"orbital_elements = [\n /*{\n index: ,\n name: \"\",\n relative_mass: , // Relative to the Sun\n semi_major_axis: , // in AU\n angular_velocity: , // °/year\n eccentricity: ,\n longitude_perihelion: , // °\n inclination: , // °\n longitude_ascending_node: ,\n }, */\n {\n index: 5,\n name: \"Jupiter\",\n relative_mass: 0.000954786, // Relative to the Sun\n semi_major_axis: 5.202545, // in AU\n angular_velocity: 30.3374, // °/year\n eccentricity: 0.0474622,\n longitude_perihelion: 13.983865, // °\n inclination: 1.30667, // °\n longitude_ascending_node: 100.0381,\n }, \n {\n index: 6,\n name: \"Saturn\",\n relative_mass: 0.000285837, // Relative to the Sun\n semi_major_axis: 9.554841, // in AU\n angular_velocity: 12.1890, // °/year\n eccentricity: 0.0575481,\n longitude_perihelion: 88.719425, // °\n inclination: 2.48795, // °\n longitude_ascending_node: 113.1334,\n },\n]","pinned":false,"hidden":false,"mode":"ojs"},{"id":1335,"value":"color_Sun = 'orange';","pinned":false,"hidden":false,"mode":"ojs"},{"id":1337,"value":"color_Jupiter = '#E060B0'","pinned":false,"hidden":false,"mode":"ojs"},{"id":3822,"value":"color_Jupiter2 = '#EC56B5'","pinned":false,"hidden":false,"mode":"ojs"},{"id":1340,"value":"color_Saturn = '#8090B0'","pinned":false,"hidden":false,"mode":"ojs"},{"id":3825,"value":"color_Saturn2 = '#45526D'","pinned":false,"hidden":false,"mode":"ojs"},{"id":3654,"value":"// Eccentricity e_5(t)\nplot_e5 = {\n var data = []\n var t, y\n for (let i = 0; i <= 200; i++) {\n t = -100000 + i*1000;\n y = math.pow(math.pow(e_55,2) + math.pow(e_56,2) +2*e_55*e_56*math.cos(((g_6-g_5)*t*math.PI/180 + (beta_6-beta_5)*math.PI/180)),0.5)\n data.push({time:t,eccentricity:y})\n }\n return data\n} ","pinned":false,"hidden":false,"mode":"ojs"},{"id":3793,"value":"// Eccentricity e_6(t)\nplot_e6 = {\n var data = []\n var t, y\n for (let i = 0; i <= 200; i++) {\n t = -100000 + i*1000;\n y = math.pow(math.pow(e_66,2) + math.pow(e_65,2) +2*e_66*e_65*math.cos(((g_6-g_5)*t*math.PI/180 + (beta_6-beta_5)*math.PI/180)),0.5)\n data.push({time:t,eccentricity:y})\n }\n return data\n} ","pinned":false,"hidden":false,"mode":"ojs"},{"id":3916,"value":"// Inclination I_5(t)\nplot_I5 = {\n var data = []\n var t, y\n for (let i = 0; i <= 200; i++) {\n t = -100000 + i*1000;\n y = (180/math.PI)*math.pow(math.pow(I_55,2) + math.pow(I_56,2) +2*I_55*I_56*math.cos(((f_6-f_5)*t*math.PI/180 + (gamma_6-gamma_5)*math.PI/180)),0.5)\n data.push({time:t,inclination:y})\n }\n return data\n} ","pinned":false,"hidden":false,"mode":"ojs"},{"id":3933,"value":"// Inclination I_6(t)\nplot_I6 = {\n var data = []\n var t, y\n for (let i = 0; i <= 200; i++) {\n t = -100000 + i*1000;\n y = (180/math.PI)*math.pow(math.pow(I_66,2) + math.pow(I_65,2) +2*I_66*I_65*math.cos(((f_6-f_5)*t*math.PI/180 + (gamma_6-gamma_5)*math.PI/180)),0.5)\n data.push({time:t,inclination:y})\n }\n return data\n} ","pinned":false,"hidden":false,"mode":"ojs"},{"id":2749,"value":"alpha_56 = Jupiter.semi_major_axis/Saturn.semi_major_axis","pinned":false,"hidden":false,"mode":"ojs"},{"id":3345,"value":"e_55 = S_5*eigs_A[1][0][0]","pinned":false,"hidden":false,"mode":"ojs"},{"id":3347,"value":"e_65 = S_5*eigs_A[1][1][0]","pinned":false,"hidden":false,"mode":"ojs"},{"id":3352,"value":"e_56 = S_6*eigs_A[1][0][1]","pinned":false,"hidden":false,"mode":"ojs"},{"id":3357,"value":"e_66 = S_6*eigs_A[1][1][1]","pinned":false,"hidden":false,"mode":"ojs"},{"id":3285,"value":"beta_5 = math.atan2(x[0],x[2])* 180 / Math.PI;","pinned":false,"hidden":false,"mode":"ojs"},{"id":3288,"value":"beta_6 = math.atan2(x[1],x[3])* 180 / Math.PI;","pinned":false,"hidden":false,"mode":"ojs"},{"id":3366,"value":"I_55 = T_5*eigs_B[1][0][0]","pinned":false,"hidden":false,"mode":"ojs"},{"id":3369,"value":"I_65 = T_5*eigs_B[1][1][0]","pinned":false,"hidden":false,"mode":"ojs"},{"id":3377,"value":"I_56 = T_6*eigs_B[1][0][1]","pinned":false,"hidden":false,"mode":"ojs"},{"id":3380,"value":"I_66 = T_6*eigs_B[1][1][1]","pinned":false,"hidden":false,"mode":"ojs"},{"id":3315,"value":"gamma_5 = math.atan2(x2[0],x2[2])* 180 / Math.PI;","pinned":false,"hidden":false,"mode":"ojs"},{"id":3319,"value":"gamma_6 = math.atan2(x2[1],x2[3])* 180 / Math.PI;","pinned":false,"hidden":false,"mode":"ojs"},{"id":3260,"value":"S_5 = math.pow(math.pow(x[0],2)+math.pow(x[2],2),0.5)","pinned":false,"hidden":false,"mode":"ojs"},{"id":3265,"value":"S_6 = math.pow(math.pow(x[1],2)+math.pow(x[3],2),0.5)","pinned":false,"hidden":false,"mode":"ojs"},{"id":3306,"value":"T_5 = math.pow(math.pow(x2[0],2)+math.pow(x2[2],2),0.5)","pinned":false,"hidden":false,"mode":"ojs"},{"id":3311,"value":"T_6 = math.pow(math.pow(x2[1],2)+math.pow(x2[3],2),0.5)","pinned":false,"hidden":false,"mode":"ojs"},{"id":3218,"value":"x = math.multiply(math.inv(math.multiply(math.transpose(M), M)), math.transpose(M), math.transpose([h_5,h_6,k_5,k_6]));","pinned":false,"hidden":false,"mode":"ojs"},{"id":3303,"value":"x2 = math.multiply(math.inv(math.multiply(math.transpose(N), N)), math.transpose(N), math.transpose([p_5,p_6,q_5,q_6]));","pinned":false,"hidden":false,"mode":"ojs"},{"id":3202,"value":"M = [ [ eigs_A[1][0][0], eigs_A[1][0][1], 0, 0 ],\n [ eigs_A[1][1][0], eigs_A[1][1][1], 0, 0 ],\n [ 0, 0, eigs_A[1][0][0], eigs_A[1][0][1] ],\n [ 0, 0, eigs_A[1][1][0], eigs_A[1][1][1] ]]","pinned":false,"hidden":false,"mode":"ojs"},{"id":3300,"value":"N = [ [ eigs_B[1][0][0], eigs_B[1][0][1], 0, 0 ],\n [ eigs_B[1][1][0], eigs_B[1][1][1], 0, 0 ],\n [ 0, 0, eigs_B[1][0][0], eigs_B[1][0][1] ],\n [ 0, 0, eigs_B[1][1][0], eigs_B[1][1][1] ]]","pinned":false,"hidden":false,"mode":"ojs"},{"id":3046,"value":"h_5 = Jupiter.eccentricity*math.sin(Jupiter.longitude_perihelion*math.PI/180)","pinned":false,"hidden":false,"mode":"ojs"},{"id":3062,"value":"k_5 = Jupiter.eccentricity*math.cos(Jupiter.longitude_perihelion*math.PI/180)","pinned":false,"hidden":false,"mode":"ojs"},{"id":3066,"value":"h_6 = Saturn.eccentricity*math.sin(Saturn.longitude_perihelion*math.PI/180)","pinned":false,"hidden":false,"mode":"ojs"},{"id":3069,"value":"k_6 = Saturn.eccentricity*math.cos(Saturn.longitude_perihelion*math.PI/180)","pinned":false,"hidden":false,"mode":"ojs"},{"id":3090,"value":"p_5 = Jupiter.inclination*(math.PI/180)*math.sin(Jupiter.longitude_ascending_node*math.PI/180)","pinned":false,"hidden":false,"mode":"ojs"},{"id":3096,"value":"q_5 = Jupiter.inclination*(math.PI/180)*math.cos(Jupiter.longitude_ascending_node*math.PI/180)","pinned":false,"hidden":false,"mode":"ojs"},{"id":3093,"value":"p_6 = Saturn.inclination*(math.PI/180)*math.sin(Saturn.longitude_ascending_node*math.PI/180)","pinned":false,"hidden":false,"mode":"ojs"},{"id":3100,"value":"q_6 = Saturn.inclination*(math.PI/180)*math.cos(Saturn.longitude_ascending_node*math.PI/180)","pinned":false,"hidden":false,"mode":"ojs"},{"id":3414,"value":"g_5 = eigs_A[0][0]","pinned":false,"hidden":false,"mode":"ojs"},{"id":3418,"value":"g_6 = eigs_A[0][1]","pinned":false,"hidden":false,"mode":"ojs"},{"id":3424,"value":"f_5 = eigs_B[0][0]","pinned":false,"hidden":false,"mode":"ojs"},{"id":3429,"value":"f_6 = eigs_B[0][1]","pinned":false,"hidden":false,"mode":"ojs"},{"id":2879,"value":"eigs_A = {\n const prec = 1e-11 // Reduce prec to ensure convergence of the eigenvalues\n const rounding = 5 // Rounding of the eigenvalues \n const ans = math.eigs([[A_55,A_56],[A_65,A_66]],prec) // returns {values: [E1,E2...sorted], vectors: [v1,v2.... corresponding vectors as columns]}\n let g = ans.values; // eigenvalues in °/y\n let U = ans.vectors; // ans.vectors is a matrix (vectors as columns)\n return [g,U]\n\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":2920,"value":"eigs_B = {\n const prec = 1e-11 // Reduce prec to ensure convergence of the eigenvalues\n const rounding = 5 // Rounding of the eigenvalues \n const ans = math.eigs([[B_55,B_56],[B_65,B_66]],prec) // returns {values: [E1,E2...sorted], vectors: [v1,v2.... corresponding vectors as columns]}\n let g = ans.values; // eigenvalues in °/y\n let U = ans.vectors; // ans.vectors is a matrix (vectors as columns)\n return [g,U]\n\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":2795,"value":"A_55 = \n(Jupiter.angular_velocity/4)*Saturn.relative_mass/(1+Jupiter.relative_mass)*math.pow(alpha_56,2)*b_m_s_alpha(1,1.5,alpha_56)[0];","pinned":false,"hidden":false,"mode":"ojs"},{"id":2799,"value":"A_56 = - (Jupiter.angular_velocity/4)*Saturn.relative_mass/(1+Jupiter.relative_mass)*math.pow(alpha_56,2)*b_m_s_alpha(2,1.5,alpha_56)[0];","pinned":false,"hidden":false,"mode":"ojs"},{"id":2801,"value":"A_65 = -(Saturn.angular_velocity/4)*(Jupiter.relative_mass/(1+Saturn.relative_mass))*alpha_56*b_m_s_alpha(2,1.5,alpha_56)[0];","pinned":false,"hidden":false,"mode":"ojs"},{"id":2803,"value":"A_66 = (Saturn.angular_velocity/4)*(Jupiter.relative_mass/(1+Saturn.relative_mass))*alpha_56*b_m_s_alpha(1,1.5,alpha_56)[0];","pinned":false,"hidden":false,"mode":"ojs"},{"id":2859,"value":"B_55 = -(Jupiter.angular_velocity/4)*(Saturn.relative_mass/(1+Jupiter.relative_mass))*math.pow(alpha_56,2)*b_m_s_alpha(1,1.5,alpha_56)[0];","pinned":false,"hidden":false,"mode":"ojs"},{"id":2864,"value":"B_56 = (Jupiter.angular_velocity/4)*(Saturn.relative_mass/(1+Jupiter.relative_mass))*math.pow(alpha_56,2)*b_m_s_alpha(1,1.5,alpha_56)[0];","pinned":false,"hidden":false,"mode":"ojs"},{"id":2868,"value":"B_65 = (Saturn.angular_velocity/4)*(Jupiter.relative_mass/(1+Saturn.relative_mass))*alpha_56*b_m_s_alpha(1,1.5,alpha_56)[0];","pinned":false,"hidden":false,"mode":"ojs"},{"id":2873,"value":"B_66 = -(Saturn.angular_velocity/4)*(Jupiter.relative_mass/(1+Saturn.relative_mass))*alpha_56*b_m_s_alpha(1,1.5,alpha_56)[0];","pinned":false,"hidden":false,"mode":"ojs"},{"id":2706,"value":"Jupiter = planet_data(5)","pinned":false,"hidden":false,"mode":"ojs"},{"id":2708,"value":"Saturn = planet_data(6)","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W34","stdlib":"1","titling":"auto"}