对火星轨道变化问题的最后解释
书名：死在火星上 作者：天瑞说符
对火星轨道变化问题的最后解释
作者君在作品相关中其实已经解释过这个问题。
不过仍然有人质疑——“你说得太含糊了”，“火星轨道的变化比你想象要大得多！”
那好吧，既然作者君的简单解释不够有力，那咱们就看看严肃的东西，反正这本书写到现在，嚷嚷着本书BUG一大堆，用初高中物理在书中挑刺的人也不少。
以下是文章内容：
Long-term
integrations
and
stability
of
planetary
orbits
in
our
Solar
system
Abstract
We
present
the
results
of
very
long-term
numerical
integrations
of
planetary
orbital
motions
over
109
-yr
time-spans
including
all
nine
planets.
A
quick
inspection
of
our
numerical
data
shows
that
the
planetary
motion，
at
least
in
our
simple
dynamical
model，
seems
to
be
quite
stable
even
over
this
very
long
time-span.
A
closer
look
at
the
lowest-frequency
oscillations
using
a
low-pass
filter
shows
us
the
potentially
diffusive
character
of
terrestrial
planetary
motion，
especially
that
of
Mercury.
The
behaviour
of
the
eccentricity
of
Mercury
in
our
integrations
is
qualitatively
similar
to
the
results
from
Jacques
Laskar's
secular
perturbation
theory
(e.g.
emax∼
0.35
over
∼±
4
Gyr).
However，
there
are
no
apparent
secular
increases
of
eccentricity
or
inclination
in
any
orbital
elements
of
the
planets，
which
may
be
revealed
by
still
longer-term
numerical
integrations.
We
have
also
performed
a
couple
of
trial
integrations
including
motions
of
the
outer
five
planets
over
the
duration
of
±
5
×
1010
yr.
The
result
indicates
that
the
three
major
resonances
in
the
Neptune–Pluto
system
have
been
maintained
over
the
1011-yr
time-span.
1
Introduction
1.1Definition
of
the
problem
The
question
of
the
stability
of
our
Solar
system
has
been
debated
over
several
hundred
years，
since
the
era
of
Newton.
The
problem
has
attracted
many
famous
mathematicians
over
the
years
and
has
played
a
central
role
in
the
development
of
non-linear
dynamics
and
chaos
theory.
However，
we
do
not
yet
have
a
definite
answer
to
the
question
of
whether
our
Solar
system
is
stable
or
not.
This
is
partly
a
result
of
the
fact
that
the
definition
of
the
term
‘stability’
is
vague
when
it
is
used
in
relation
to
the
problem
of
planetary
motion
in
the
Solar
system.
Actually
it
is
not
easy
to
give
a
clear，
rigorous
and
physically
meaningful
definition
of
the
stability
of
our
Solar
system.
Among
many
definitions
of
stability，
here
we
adopt
the
Hill
definition
(Gladman
1993):
actually
this
is
not
a
definition
of
stability，
but
of
instability.
We
define
a
system
as
becoming
unstable
when
a
close
encounter
occurs
somewhere
in
the
system，
starting
from
a
certain
initial
configuration
(Chambers，
Wetherill
&
Boss
1996;
Ito
&
Tanikawa
1999).
A
system
is
defined
as
experiencing
a
close
encounter
when
two
bodies
approach
one
another
within
an
area
of
the
larger
Hill
radius.
Otherwise
the
system
is
defined
as
being
stable.
Henceforward
we
state
that
our
planetary
system
is
dynamically
stable
if
no
close
encounter
happens
during
the
age
of
our
Solar
system，
about
±5
Gyr.
Incidentally，
this
definition
may
be
replaced
by
one
in
which
an
occurrence
of
any
orbital
crossing
between
either
of
a
pair
of
planets
takes
place.
This
is
because
we
know
from
experience
that
an
orbital
crossing
is
very
likely
to
lead
to
a
close
encounter
in
planetary
and
protoplanetary
systems
(Yoshinaga，
Kokubo
&
Makino
1999).
Of
course
this
statement
cannot
be
simply
applied
to
systems
with
stable
orbital
resonances
such
as
the
Neptune–Pluto
system.
1.2Previous
studies
and
aims
of
this
research
In
addition
to
the
vagueness
of
the
concept
of
stability，
the
planets
in
our
Solar
system
show
a
character
typical
of
dynamical
chaos
(Sussman
&
Wisdom
1988，
1992).
The
cause
of
this
chaotic
behaviour
is
now
partly
understood
as
being
a
result
of
resonance
overlapping
(Murray
&
Holman
1999;
Lecar，
Franklin
&
Holman
2001).
However，
it
would
require
integrating
over
an
ensemble
of
planetary
systems
including
all
nine
planets
for
a
period
covering
several
10
Gyr
to
thoroughly
understand
the
long-term
evolution
of
planetary
orbits，
since
chaotic
dynamical
systems
are
characterized
by
their
strong
dependence
on
initial
conditions.
From
that
point
of
view，
many
of
the
previous
long-term
numerical
integrations
included
only
the
outer
five
planets
(Sussman
&
Wisdom
1988;
Kinoshita
&
Nakai
1996).
This
is
because
the
orbital
periods
of
the
outer
planets
are
so
much
longer
than
those
of
the
inner
four
planets
that
it
is
much
easier
to
follow
the
system
for
a
given
integration
period.
At
present，
the
longest
numerical
integrations
published
in
journals
are
those
of
Duncan
&
Lissauer
(1998).
Although
their
main
target
was
the
effect
of
post-main-sequence
solar
mass
loss
on
the
stability
of
planetary
orbits，
they
performed
many
integrations
covering
up
to
∼1011
yr
of
the
orbital
motions
of
the
four
jovian
planets.
The
initial
orbital
elements
and
masses
of
planets
are
the
same
as
those
of
our
Solar
system
in
Duncan
&
Lissauer's
paper，
but
they
decrease
the
mass
of
the
Sun
gradually
in
their
numerical
experiments.
This
is
because
they
consider
the
effect
of
post-main-sequence
solar
mass
loss
in
the
paper.
Consequently，
they
found
that
the
crossing
time-scale
of
planetary
orbits，
which
can
be
a
typical
indicator
of
the
instability
time-scale，
is
quite
sensitive
to
the
rate
of
mass
decrease
of
the
Sun.
When
the
mass
of
the
Sun
is
close
to
its
present
value，
the
jovian
planets
remain
stable
over
1010
yr，
or
perhaps
longer.
Duncan
&
Lissauer
also
performed
four
similar
experiments
on
the
orbital
motion
of
seven
planets
(Venus
to
Neptune)，
which
cover
a
span
of
∼109
yr.
Their
experiments
on
the
seven
planets
are
not
yet
comprehensive，
but
it
seems
that
the
terrestrial
planets
also
remain
stable
during
the
integration
period，
maintaining
almost
regular
oscillations.
On
the
other
hand，
in
his
accurate
semi-analytical
secular
perturbation
theory
(Laskar
1988)，
Laskar
finds
that
large
and
irregular
variations
can
appear
in
the
eccentricities
and
inclinations
of
the
terrestrial
planets，
especially
of
Mercury
and
Mars
on
a
time-scale
of
several
109
yr
(Laskar
1996).
The
results
of
Laskar's
secular
perturbation
theory
should
be
confirmed
and
investigated
by
fully
numerical
integrations.
In
this
paper
we
present
preliminary
results
of
six
long-term
numerical
integrations
on
all
nine
planetary
orbits，
covering
a
span
of
several
109
yr，
and
of
two
other
integrations
covering
a
span
of
±
5
×
1010
yr.
The
total
elapsed
time
for
all
integrations
is
more
than
5
yr，
using
several
dedicated
PCs
and
workstations.
One
of
the
fundamental
conclusions
of
our
long-term
integrations
is
that
Solar
system
planetary
motion
seems
to
be
stable
in
terms
of
the
Hill
stability
mentioned
above，
at
least
over
a
time-span
of
±
4
Gyr.
Actually，
in
our
numerical
integrations
the
system
was
far
more
stable
than
what
is
defined
by
the
Hill
stability
criterion:
not
only
did
no
close
encounter
happen
during
the
integration
period，
but
also
all
the
planetary
orbital
elements
have
been
confined
in
a
narrow
region
both
in
time
and
frequency
domain，
though
planetary
motions
are
stochastic.
Since
the
purpose
of
this
paper
is
to
exhibit
and
overview
the
results
of
our
long-term
numerical
integrations，
we
show
typical
example
figures
as
evidence
of
the
very
long-term
stability
of
Solar
system
planetary
motion.
For
readers
who
have
more
specific
and
deeper
interests
in
our
numerical
results，
we
have
prepared
a
webpage
(access
)，
where
we
show
raw
orbital
elements，
their
low-pass
filtered
results，
variation
of
Delaunay
elements
and
angular
momentum
deficit，
and
results
of
our
simple
time–frequency
analysis
on
all
of
our
integrations.
In
Section
2
we
briefly
explain
our
dynamical
model，
numerical
method
and
initial
conditions
used
in
our
integrations.
Section
3
is
devoted
to
a
description
of
the
quick
results
of
the
numerical
integrations.
Very
long-term
stability
of
Solar
system
planetary
motion
is
apparent
both
in
planetary
positions
and
orbital
elements.
A
rough
estimation
of
numerical
errors
is
also
given.
Section
4
goes
on
to
a
discussion
of
the
longest-term
variation
of
planetary
orbits
using
a
low-pass
filter
and
includes
a
discussion
of
angular
momentum
deficit.
In
Section
5，
we
present
a
set
of
numerical
integrations
for
the
outer
five
planets
that
spans
±
5
×
1010
yr.
In
Section
6
we
also
discuss
the
long-term
stability
of
the
planetary
motion
and
its
possible
cause.
2
Description
of
the
numerical
integrations
（本部分涉及比较复杂的积分计算，作者君就不贴上来了，贴上来了起点也不一定能成功显示。）
2.3
Numerical
method
We
utilize
a
second-order
Wisdom–Holman
symplectic
map
as
our
main
integration
method
(Wisdom
&
Holman
1991;
Kinoshita，
Yoshida
&
Nakai
1991)
with
a
special
start-up
procedure
to
reduce
the
truncation
error
of
angle
variables，‘warm
start’(Saha
&
Tremaine
1992，
1994).
The
stepsize
for
the
numerical
integrations
is
8
d
throughout
all
integrations
of
the
nine
planets
(N±1，2，3)，
which
is
about
1/11
of
the
orbital
period
of
the
innermost
planet
(Mercury).
As
for
the
determination
of
stepsize，
we
partly
follow
the
previous
numerical
integration
of
all
nine
planets
in
Sussman
&
Wisdom
(1988，
7.2
d)
and
Saha
&
Tremaine
(1994，
225/32
d).
We
rounded
the
decimal
part
of
the
their
stepsizes
to
8
to
make
the
stepsize
a
multiple
of
2
in
order
to
reduce
the
accumulation
of
round-off
error
in
the
computation
processes.
In
relation
to
this，
Wisdom
&
Holman
(1991)
performed
numerical
integrations
of
the
outer
five
planetary
orbits
using
the
symplectic
map
with
a
stepsize
of
400
d，
1/10.83
of
the
orbital
period
of
Jupiter.
Their
result
seems
to
be
accurate
enough，
which
partly
justifies
our
method
of
determining
the
stepsize.
However，
since
the
eccentricity
of
Jupiter
(∼0.05)
is
much
smaller
than
that
of
Mercury
(∼0.2)，
we
need
some
care
when
we
compare
these
integrations
simply
in
terms
of
stepsizes.
In
the
integration
of
the
outer
five
planets
(F±)，
we
fixed
the
stepsize
at
400
d.
We
adopt
Gauss'
f
and
g
functions
in
the
symplectic
map
together
with
the
third-order
Halley
method
(Danby
1992)
as
a
solver
for
Kepler
equations.
The
number
of
maximum
iterations
we
set
in
Halley's
method
is
15，
but
they
never
reached
the
maximum
in
any
of
our
integrations.
The
interval
of
the
data
output
is
200
000
d
(∼547
yr)
for
the
calculations
of
all
nine
planets
(N±1，2，3)，
and
about
8000
000
d
(∼21
903
yr)
for
the
integration
of
the
outer
five
planets
(F±).
Although
no
output
filtering
was
done
when
the
numerical
integrations
were
in
process，
we
applied
a
low-pass
filter
to
the
raw
orbital
data
after
we
had
completed
all
the
calculations.
See
Section
4.1
for
more
detail.
2.4
Error
estimation
2.4.1
Relative
errors
in
total
energy
and
angular
momentum
According
to
one
of
the
basic
properties
of
symplectic
integrators，
which
conserve
the
physically
conservative
quantities
well
(total
orbital
energy
and
angular
momentum)，
our
long-term
numerical
integrations
seem
to
have
been
performed
with
very
small
errors.
The
averaged
relative
errors
of
total
energy
(∼10−9)
and
of
total
angular
momentum
(∼10−11)
have
remained
nearly
constant
throughout
the
integration
period
(Fig.
1).
The
special
startup
procedure，
warm
start，
would
have
reduced
the
averaged
relative
error
in
total
energy
by
about
one
order
of
magnitude
or
more.
Relative
numerical
error
of
the
total
angular
momentum
δA/A0
and
the
total
energy
δE/E0
in
our
numerical
integrationsN±
1，2，3，
where
δE
and
δA
are
the
absolute
change
of
the
total
energy
and
total
angular
momentum，
respectively，
andE0andA0are
their
initial
values.
The
horizontal
unit
is
Gyr.
Note
that
different
operating
systems，
different
mathematical
libraries，
and
different
hardware
architectures
result
in
different
numerical
errors，
through
the
variations
in
round-off
error
handling
and
numerical
algorithms.
In
the
upper
panel
of
Fig.
1，
we
can
recognize
this
situation
in
the
secular
numerical
error
in
the
total
angular
momentum，
which
should
be
rigorously
preserved
up
to
machine-ε
precision.
2.4.2
Error
in
planetary
longitudes
Since
the
symplectic
maps
preserve
total
energy
and
total
angular
momentum
of
N-body
dynamical
systems
inherently
well，
the
degree
of
their
preservation
may
not
be
a
good
measure
of
the
accuracy
of
numerical
integrations，
especially
as
a
measure
of
the
positional
error
of
planets，
i.e.
the
error
in
planetary
longitudes.
To
estimate
the
numerical
error
in
the
planetary
longitudes，
we
performed
the
following
procedures.
We
compared
the
result
of
our
main
long-term
integrations
with
some
test
integrations，
which
span
much
shorter
periods
but
with
much
higher
accuracy
than
the
main
integrations.
For
this
purpose，
we
performed
a
much
more
accurate
integration
with
a
stepsize
of
0.125
d
(1/64
of
the
main
integrations)
spanning
3
×
105
yr，
starting
with
the
same
initial
conditions
as
in
the
N−1
integration.
We
consider
that
this
test
integration
provides
us
with
a
‘pseudo-true’
solution
of
planetary
orbital
evolution.
Next，
we
compare
the
test
integration
with
the
main
integration，
N−1.
For
the
period
of
3
×
105
yr，
we
see
a
difference
in
mean
anomalies
of
the
Earth
between
the
two
integrations
of
∼0.52°(in
the
case
of
the
N−1
integration).
This
difference
can
be
extrapolated
to
the
value
∼8700°，
about
25
rotations
of
Earth
after
5
Gyr，
since
the
error
of
longitudes
increases
linearly
with
time
in
the
symplectic
map.
Similarly，
the
longitude
error
of
Pluto
can
be
estimated
as
∼12°.
This
value
for
Pluto
is
much
better
than
the
result
in
Kinoshita
&
Nakai
(1996)
where
the
difference
is
estimated
as
∼60°.
3
Numerical
results
–
I.
Glance
at
the
raw
data
In
this
section
we
briefly
review
the
long-term
stability
of
planetary
orbital
motion
through
some
snapshots
of
raw
numerical
data.
The
orbital
motion
of
planets
indicates
long-term
stability
in
all
of
our
numerical
integrations:
no
orbital
crossings
nor
close
encounters
between
any
pair
of
planets
took
place.
3.1
General
description
of
the
stability
of
planetary
orbits
First，
we
briefly
look
at
the
general
character
of
the
long-term
stability
of
planetary
orbits.
Our
interest
here
focuses
particularly
on
the
inner
four
terrestrial
planets
for
which
the
orbital
time-scales
are
much
shorter
than
those
of
the
outer
five
planets.
As
we
can
see
clearly
from
the
planar
orbital
configurations
shown
in
Figs
2
and
3，
orbital
positions
of
the
terrestrial
planets
differ
little
between
the
initial
and
final
part
of
each
numerical
integration，
which
spans
several
Gyr.
The
solid
lines
denoting
the
present
orbits
of
the
planets
lie
almost
within
the
swarm
of
dots
even
in
the
final
part
of
integrations
(b)
and
(d).
This
indicates
that
throughout
the
entire
integration
period
the
almost
regular
variations
of
planetary
orbital
motion
remain
nearly
the
same
as
they
are
at
present.
Vertical
view
of
the
four
inner
planetary
orbits
(from
the
z
-axis
direction)
at
the
initial
and
final
parts
of
the
integrationsN±1.
The
axes
units
are
au.
The
xy
-plane
is
set
to
the
invariant
plane
of
Solar
system
total
angular
momentum.(a)
The
initial
part
ofN+1
(
t
=
0
to
0.0547
×
10
9
yr).(b)
The
final
part
ofN+1
(
t
=
4.9339
×
10
8
to
4.9886
×
10
9
yr).(c)
The
initial
part
of
N−1
(t=
0
to
−0.0547
×
109
yr).(d)
The
final
part
ofN−1
(
t
=−3.9180
×
10
9
to
−3.9727
×
10
9
yr).
In
each
panel，
a
total
of
23
684
points
are
plotted
with
an
interval
of
about
2190
yr
over
5.47
×
107
yr
.
Solid
lines
in
each
panel
denote
the
present
orbits
of
the
four
terrestrial
planets
(taken
from
DE245).
The
variation
of
eccentricities
and
orbital
inclinations
for
the
inner
four
planets
in
the
initial
and
final
part
of
the
integration
N+1
is
shown
in
Fig.
4.
As
expected，
the
character
of
the
variation
of
planetary
orbital
elements
does
not
differ
significantly
between
the
initial
and
final
part
of
each
integration，
at
least
for
Venus，
Earth
and
Mars.
The
elements
of
Mercury，
especially
its
eccentricity，
seem
to
change
to
a
significant
extent.
This
is
partly
because
the
orbital
time-scale
of
the
planet
is
the
shortest
of
all
the
planets，
which
leads
to
a
more
rapid
orbital
evolution
than
other
planets;
the
innermost
planet
may
be
nearest
to
instability.
This
result
appears
to
be
in
some
agreement
with
Laskar's
(1994，
1996)
expectations
that
large
and
irregular
variations
appear
in
the
eccentricities
and
inclinations
of
Mercury
on
a
time-scale
of
several
109
yr.
However，
the
effect
of
the
possible
instability
of
the
orbit
of
Mercury
may
not
fatally
affect
the
global
stability
of
the
whole
planetary
system
owing
to
the
small
mass
of
Mercury.
We
will
mention
briefly
the
long-term
orbital
evolution
of
Mercury
later
in
Section
4
using
low-pass
filtered
orbital
elements.
The
orbital
motion
of
the
outer
five
planets
seems
rigorously
stable
and
quite
regular
over
this
time-span
(see
also
Section
5).
3.2
Time–frequency
maps
Although
the
planetary
motion
exhibits
very
long-term
stability
defined
as
the
non-existence
of
close
encounter
events，
the
chaotic
nature
of
planetary
dynamics
can
change
the
oscillatory
period
and
amplitude
of
planetary
orbital
motion
gradually
over
such
long
time-spans.
Even
such
slight
fluctuations
of
orbital
variation
in
the
frequency
domain，
particularly
in
the
case
of
Earth，
can
potentially
have
a
significant
effect
on
its
surface
climate
system
through
solar
insolation
variation
(cf.
Berger
1988).
To
give
an
overview
of
the
long-term
change
in
periodicity
in
planetary
orbital
motion，
we
performed
many
fast
Fourier
transformations
(FFTs)
along
the
time
axis，
and
superposed
the
resulting
periodgrams
to
draw
two-dimensional
time–frequency
maps.
The
specific
approach
to
drawing
these
time–frequency
maps
in
this
paper
is
very
simple
–
much
simpler
than
the
wavelet
analysis
or
Laskar's
(1990，
1993)
frequency
analysis.
Divide
the
low-pass
filtered
orbital
data
into
many
fragments
of
the
same
length.
The
length
of
each
data
segment
should
be
a
multiple
of
2
in
order
to
apply
the
FFT.
Each
fragment
of
the
data
has
a
large
overlapping
part:
for
example，
when
the
ith
data
begins
from
t=ti
and
ends
at
t=ti+T，
the
next
data
segment
ranges
from
ti+δT≤ti+δT+T，
where
δT?T.
We
continue
this
division
until
we
reach
a
certain
number
N
by
which
tn+T
reaches
the
total
integration
length.
We
apply
an
FFT
to
each
of
the
data
fragments，
and
obtain
n
frequency
diagrams.
In
each
frequency
diagram
obtained
above，
the
strength
of
periodicity
can
be
replaced
by
a
grey-scale
(or
colour)
chart.
We
perform
the
replacement，
and
connect
all
the
grey-scale
(or
colour)
charts
into
one
graph
for
each
integration.
The
horizontal
axis
of
these
new
graphs
should
be
the
time，
i.e.
the
starting
times
of
each
fragment
of
data
(ti，
where
i=
1，…，
n).
The
vertical
axis
represents
the
period
(or
frequency)
of
the
oscillation
of
orbital
elements.
We
have
adopted
an
FFT
because
of
its
overwhelming
speed，
since
the
amount
of
numerical
data
to
be
decomposed
into
frequency
components
is
terribly
huge
(several
tens
of
Gbytes).
A
typical
example
of
the
time–frequency
map
created
by
the
above
procedures
is
shown
in
a
grey-scale
diagram
as
Fig.
5，
which
shows
the
variation
of
periodicity
in
the
eccentricity
and
inclination
of
Earth
in
N+2
integration.
In
Fig.
5，
the
dark
area
shows
that
at
the
time
indicated
by
the
value
on
the
abscissa，
the
periodicity
indicated
by
the
ordinate
is
stronger
than
in
the
lighter
area
around
it.
We
can
recognize
from
this
map
that
the
periodicity
of
the
eccentricity
and
inclination
of
Earth
only
changes
slightly
over
the
entire
period
covered
by
the
N+2
integration.
This
nearly
regular
trend
is
qualitatively
the
same
in
other
integrations
and
for
other
planets，
although
typical
frequencies
differ
planet
by
planet
and
element
by
element.
4.2
Long-term
exchange
of
orbital
energy
and
angular
momentum
We
calculate
very
long-periodic
variation
and
exchange
of
planetary
orbital
energy
and
angular
momentum
using
filtered
Delaunay
elements
L，
G，
H.
G
and
H
are
equivalent
to
the
planetary
orbital
angular
momentum
and
its
vertical
component
per
unit
mass.
L
is
related
to
the
planetary
orbital
energy
E
per
unit
mass
as
E=−μ2/2L2.
If
the
system
is
completely
linear，
the
orbital
energy
and
the
angular
momentum
in
each
frequency
bin
must
be
constant.
Non-linearity
in
the
planetary
system
can
cause
an
exchange
of
energy
and
angular
momentum
in
the
frequency
domain.
The
amplitude
of
the
lowest-frequency
oscillation
should
increase
if
the
system
is
unstable
and
breaks
down
gradually.
However，
such
a
symptom
of
instability
is
not
prominent
in
our
long-term
integrations.
In
Fig.
7，
the
total
orbital
energy
and
angular
momentum
of
the
four
inner
planets
and
all
nine
planets
are
shown
for
integration
N+2.
The
upper
three
panels
show
the
long-periodic
variation
of
total
energy
(denoted
asE-
E0)，
total
angular
momentum
(
G-
G0)，
and
the
vertical
component
(
H-
H0)
of
the
inner
four
planets
calculated
from
the
low-pass
filtered
Delaunay
elements.E0，
G0，
H0
denote
the
initial
values
of
each
quantity.
The
absolute
difference
from
the
initial
values
is
plotted
in
the
panels.
The
lower
three
panels
in
each
figure
showE-E0，G-G0
andH-H0
of
the
total
of
nine
planets.
The
fluctuation
shown
in
the
lower
panels
is
virtually
entirely
a
result
of
the
massive
jovian
planets.
Comparing
the
variations
of
energy
and
angular
momentum
of
the
inner
four
planets
and
all
nine
planets，
it
is
apparent
that
the
amplitudes
of
those
of
the
inner
planets
are
much
smaller
than
those
of
all
nine
planets:
the
amplitudes
of
the
outer
five
planets
are
much
larger
than
those
of
the
inner
planets.
This
does
not
mean
that
the
inner
terrestrial
planetary
subsystem
is
more
stable
than
the
outer
one:
this
is
simply
a
result
of
the
relative
smallness
of
the
masses
of
the
four
terrestrial
planets
compared
with
those
of
the
outer
jovian
planets.
Another
thing
we
notice
is
that
the
inner
planetary
subsystem
may
become
unstable
more
rapidly
than
the
outer
one
because
of
its
shorter
orbital
time-scales.
This
can
be
seen
in
the
panels
denoted
asinner
4
in
Fig.
7
where
the
longer-periodic
and
irregular
oscillations
are
more
apparent
than
in
the
panels
denoted
astotal
9.
Actually，
the
fluctuations
in
theinner
4
panels
are
to
a
large
extent
as
a
result
of
the
orbital
variation
of
the
Mercury.
However，
we
cannot
neglect
the
contribution
from
other
terrestrial
planets，
as
we
will
see
in
subsequent
sections.
4.4
Long-term
coupling
of
several
neighbouring
planet
pairs
Let
us
see
some
individual
variations
of
planetary
orbital
energy
and
angular
momentum
expressed
by
the
low-pass
filtered
Delaunay
elements.
Figs
10
and
11
show
long-term
evolution
of
the
orbital
energy
of
each
planet
and
the
angular
momentum
in
N+1
and
N−2
integrations.
We
notice
that
some
planets
form
apparent
pairs
in
terms
of
orbital
energy
and
angular
momentum
exchange.
In
particular，
Venus
and
Earth
make
a
typical
pair.
In
the
figures，
they
show
negative
correlations
in
exchange
of
energy
and
positive
correlations
in
exchange
of
angular
momentum.
The
negative
correlation
in
exchange
of
orbital
energy
means
that
the
two
planets
form
a
closed
dynamical
system
in
terms
of
the
orbital
energy.
The
positive
correlation
in
exchange
of
angular
momentum
means
that
the
two
planets
are
simultaneously
under
certain
long-term
perturbations.
Candidates
for
perturbers
are
Jupiter
and
Saturn.
Also
in
Fig.
11，
we
can
see
that
Mars
shows
a
positive
correlation
in
the
angular
momentum
variation
to
the
Venus–Earth
system.
Mercury
exhibits
certain
negative
correlations
in
the
angular
momentum
versus
the
Venus–Earth
system，
which
seems
to
be
a
reaction
caused
by
the
conservation
of
angular
momentum
in
the
terrestrial
planetary
subsystem.
It
is
not
clear
at
the
moment
why
the
Venus–Earth
pair
exhibits
a
negative
correlation
in
energy
exchange
and
a
positive
correlation
in
angular
momentum
exchange.
We
may
possibly
explain
this
through
observing
the
general
fact
that
there
are
no
secular
terms
in
planetary
semimajor
axes
up
to
second-order
perturbation
theories
(cf.
Brouwer
&
Clemence
1961;
Boccaletti
&
Pucacco
1998).
This
means
that
the
planetary
orbital
energy
(which
is
directly
related
to
the
semimajor
axis
a)
might
be
much
less
affected
by
perturbing
planets
than
is
the
angular
momentum
exchange
(which
relates
to
e).
Hence，
the
eccentricities
of
Venus
and
Earth
can
be
disturbed
easily
by
Jupiter
and
Saturn，
which
results
in
a
positive
correlation
in
the
angular
momentum
exchange.
On
the
other
hand，
the
semimajor
axes
of
Venus
and
Earth
are
less
likely
to
be
disturbed
by
the
jovian
planets.
Thus
the
energy
exchange
may
be
limited
only
within
the
Venus–Earth
pair，
which
results
in
a
negative
correlation
in
the
exchange
of
orbital
energy
in
the
pair.
As
for
the
outer
jovian
planetary
subsystem，
Jupiter–Saturn
and
Uranus–Neptune
seem
to
make
dynamical
pairs.
However，
the
strength
of
their
coupling
is
not
as
strong
compared
with
that
of
the
Venus–Earth
pair.
5
±
5
×
1010-yr
integrations
of
outer
planetary
orbits
Since
the
jovian
planetary
masses
are
much
larger
than
the
terrestrial
planetary
masses，
we
treat
the
jovian
planetary
system
as
an
independent
planetary
system
in
terms
of
the
study
of
its
dynamical
stability.
Hence，
we
added
a
couple
of
trial
integrations
that
span
±
5
×
1010
yr，
including
only
the
outer
five
planets
(the
four
jovian
planets
plus
Pluto).
The
results
exhibit
the
rigorous
stability
of
the
outer
planetary
system
over
this
long
time-span.
Orbital
configurations
(Fig.
12)，
and
variation
of
eccentricities
and
inclinations
(Fig.
13)
show
this
very
long-term
stability
of
the
outer
five
planets
in
both
the
time
and
the
frequency
domains.
Although
we
do
not
show
maps
here，
the
typical
frequency
of
the
orbital
oscillation
of
Pluto
and
the
other
outer
planets
is
almost
constant
during
these
very
long-term
integration
periods，
which
is
demonstrated
in
the
time–frequency
maps
on
our
webpage.
In
these
two
integrations，
the
relative
numerical
error
in
the
total
energy
was
∼10−6
and
that
of
the
total
angular
momentum
was
∼10−10.
5.1
Resonances
in
the
Neptune–Pluto
system
Kinoshita
&
Nakai
(1996)
integrated
the
outer
five
planetary
orbits
over
±
5.5
×
109
yr
.
They
found
that
four
major
resonances
between
Neptune
and
Pluto
are
maintained
during
the
whole
integration
period，
and
that
the
resonances
may
be
the
main
causes
of
the
stability
of
the
orbit
of
Pluto.
The
major
four
resonances
found
in
previous
research
are
as
follows.
In
the
following
description，λ
denotes
the
mean
longitude，Ω
is
the
longitude
of
the
ascending
node
and
ϖ
is
the
longitude
of
perihelion.
Subscripts
P
and
N
denote
Pluto
and
Neptune.
Mean
motion
resonance
between
Neptune
and
Pluto
(3:2).
The
critical
argument
θ1=
3
λP−
2
λN−ϖP
librates
around
180°
with
an
amplitude
of
about
80°
and
a
libration
period
of
about
2
×
104
yr.
The
argument
of
perihelion
of
Pluto
ωP=θ2=ϖP−ΩP
librates
around
90°
with
a
period
of
about
3.8
×
106
yr.
The
dominant
periodic
variations
of
the
eccentricity
and
inclination
of
Pluto
are
synchronized
with
the
libration
of
its
argument
of
perihelion.
This
is
anticipated
in
the
secular
perturbation
theory
constructed
by
Kozai
(1962).
The
longitude
of
the
node
of
Pluto
referred
to
the
longitude
of
the
node
of
Neptune，θ3=ΩP−ΩN，
circulates
and
the
period
of
this
circulation
is
equal
to
the
period
of
θ2
libration.
When
θ3
becomes
zero，
i.e.
the
longitudes
of
ascending
nodes
of
Neptune
and
Pluto
overlap，
the
inclination
of
Pluto
becomes
maximum，
the
eccentricity
becomes
minimum
and
the
argument
of
perihelion
becomes
90°.
When
θ3
becomes
180°，
the
inclination
of
Pluto
becomes
minimum，
the
eccentricity
becomes
maximum
and
the
argument
of
perihelion
becomes
90°
again.
Williams
&
Benson
(1971)
anticipated
this
type
of
resonance，
later
confirmed
by
Milani，
Nobili
&
Carpino
(1989).
An
argument
θ4=ϖP−ϖN+
3
(ΩP−ΩN)
librates
around
180°
with
a
long
period，∼
5.7
×
108
yr.
In
our
numerical
integrations，
the
resonances
(i)–(iii)
are
well
maintained，
and
variation
of
the
critical
arguments
θ1，θ2，θ3
remain
similar
during
the
whole
integration
period
(Figs
14–16
).
However，
the
fourth
resonance
(iv)
appears
to
be
different:
the
critical
argument
θ4
alternates
libration
and
circulation
over
a
1010-yr
time-scale
(Fig.
17).
This
is
an
interesting
fact
that
Kinoshita
&
Nakai's
(1995，
1996)
shorter
integrations
were
not
able
to
disclose.
6
Discussion
What
kind
of
dynamical
mechanism
maintains
this
long-term
stability
of
the
planetary
system?
We
can
immediately
think
of
two
major
features
that
may
be
responsible
for
the
long-term
stability.
First，
there
seem
to
be
no
significant
lower-order
resonances
(mean
motion
and
secular)
between
any
pair
among
the
nine
planets.
Jupiter
and
Saturn
are
close
to
a
5:2
mean
motion
resonance
(the
famous
‘great
inequality’)，
but
not
just
in
the
resonance
zone.
Higher-order
resonances
may
cause
the
chaotic
nature
of
the
planetary
dynamical
motion，
but
they
are
not
so
strong
as
to
destroy
the
stable
planetary
motion
within
the
lifetime
of
the
real
Solar
system.
The
second
feature，
which
we
think
is
more
important
for
the
long-term
stability
of
our
planetary
system，
is
the
difference
in
dynamical
distance
between
terrestrial
and
jovian
planetary
subsystems
(Ito
&
Tanikawa
1999，
2001).
When
we
measure
planetary
separations
by
the
mutual
Hill
radii
(R_)，
separations
among
terrestrial
planets
are
greater
than
26RH，
whereas
those
among
jovian
planets
are
less
than
14RH.
This
difference
is
directly
related
to
the
difference
between
dynamical
features
of
terrestrial
and
jovian
planets.
Terrestrial
planets
have
smaller
masses，
shorter
orbital
periods
and
wider
dynamical
separation.
They
are
strongly
perturbed
by
jovian
planets
that
have
larger
masses，
longer
orbital
periods
and
narrower
dynamical
separation.
Jovian
planets
are
not
perturbed
by
any
other
massive
bodies.
The
present
terrestrial
planetary
system
is
still
being
disturbed
by
the
massive
jovian
planets.
However，
the
wide
separation
and
mutual
interaction
among
the
terrestrial
planets
renders
the
disturbance
ineffective;
the
degree
of
disturbance
by
jovian
planets
is
O(eJ)(order
of
magnitude
of
the
eccentricity
of
Jupiter)，
since
the
disturbance
caused
by
jovian
planets
is
a
forced
oscillation
having
an
amplitude
of
O(eJ).
Heightening
of
eccentricity，
for
example
O(eJ)∼0.05，
is
far
from
sufficient
to
provoke
instability
in
the
terrestrial
planets
having
such
a
wide
separation
as
26RH.
Thus
we
assume
that
the
present
wide
dynamical
separation
among
terrestrial
planets
(>
26RH)
is
probably
one
of
the
most
significant
conditions
for
maintaining
the
stability
of
the
planetary
system
over
a
109-yr
time-span.
Our
detailed
analysis
of
the
relationship
between
dynamical
distance
between
planets
and
the
instability
time-scale
of
Solar
system
planetary
motion
is
now
on-going.
Although
our
numerical
integrations
span
the
lifetime
of
the
Solar
system，
the
number
of
integrations
is
far
from
sufficient
to
fill
the
initial
phase
space.
It
is
necessary
to
perform
more
and
more
numerical
integrations
to
confirm
and
examine
in
detail
the
long-term
stability
of
our
planetary
dynamics.
——以上文段引自
Ito，
T.&
Tanikawa，
K.
Long-term
integrations
and
stability
of
planetary
orbits
in
our
Solar
System.
Mon.
Not.
R.
Astron.
Soc.
336，
483–500
(2002)
这只是作者君参考的一篇文章，关于太阳系的稳定性。
还有其他论文，不过也都是英文的，相关课题的中文文献很少，那些论文下载一篇要九美元（《Nature》真是暴利），作者君写这篇文章的时候已经回家，不在检测中心，所以没有数据库的使用权，下不起，就不贴上来了。