Their  experiments  on  the  seven  s  are  not  yet  comprehensive�

�

�  but  it  seems  that  the  terrestrial  s  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  s,  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  ary  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  ary  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  ary  orbital  elements  have  been  confined  in  a  narrow  region  both  in  time  and  frequency  domain�

�

�  though  ary  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  ary  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  ary  motion  is  apparent  both  in  ary  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  ary  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  s  that  spans  ±  5  ×  1010  yr.  In  Section  6  we  also  discuss  the  long-term  stability  of  the  ary  motion  and  its  possible  cause.

2 Description  of  the  numerical  integrations

(This part involves relatively complex points calculations. The author will not post them. If you post them, the starting point may not be displayed successfully.)

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  all  integrations  of  the nine  s  (N±1�

�

�2�

�

�3)�

�

�  which is  about  1/11  of  the  orbital  period of  the  innermost    (Mercury).

As for  the  determination of  stepsize�

�

�  we  partly  follow  the  previous  numerical  integration of  all  nine s in  Sussman  &  Wisdom  (1988�

�

�  7.

2 d)  and   Saha &  Tremaine (1994�

�

�  225/32  d).

We  rounded  the  decimal  part of  the  the  stepsize  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  ary  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  s  (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  s  (N±1,2,3),  and  about  8000  000  d  (～21  903  yr)  for  the  integration  of  the  outer  five  s  (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  ary  longitudes
To be continued...