DECLARE SUB Macheps ()
DECLARE SUB Euler ()
DECLARE SUB Runge ()
DECLARE SUB RK4 ()
' r‚solutions d'‚quations diff‚rentielle
' Macheps ' pour tester la pr‚cision des calcls
' Euler
' Runge
 RK4


SUB Euler
' M‚thode de Euler
PRINT
PRINT "  M‚thode d'Euler"
GM# = -6.67E-11 * 5.97E+24
' conditions initiales
x0# = 1E+08
v0# = 0
t0# = 0
tLimit# = 55000  ' on s'int‚resse … la position en tLimit

dtBase# = 2 ^ 11' incr‚ment en temps
DO
dtBase# = dtBase# / 2
dt# = dtBase#
t# = t0#
x# = x0#
v# = v0#
DO
' pour s'arrˆter exactement … tLimit
IF (t# + dt# > tLimit#) THEN
  dt# = tLimit# - t#
  t# = tLimit#
ELSE
t# = t# + dt#
END IF

xt# = x#  ' valeur temporaire
x# = x# + v# * dt#
v# = v# + (GM# / xt# ^ 2) * dt#
LOOP UNTIL t# >= tLimit#
PRINT USING "  dt= #####.#####  err= #.###^^^^  x= #######.#########  v= ######.#########"; dtBase#; ABS(x# / 1000 - 9049.721370064#); x# / 1000; v#
LOOP UNTIL dtBase# < .03

END SUB

SUB Macheps
' calcul la pr‚cision des calculs en virgules flottantes
eps# = 1
DO
  eps# = eps# / 2
  a# = 1 + eps#
LOOP UNTIL a# = 1
PRINT "  eps="; eps#
END SUB

SUB RK4
' M‚thode de Runge-Kutta d'ordre 4
' Equation diff‚rentielle : (x, v)' = (v, GM/x^2) = F(x, v)
' F(x, v) = (v, GM/x^2)

GM# = -6.67E-11 * 5.97E+24
' conditions initiales
x0# = 1E+08
v0# = 0
t0# = 0
tLimit# = 55000   ' on s'int‚resse … la position en tLimit

PRINT
PRINT "  M‚thode de Runge-Kutta d'ordre 4"

dtBase# = 2 ^ 15  ' incr‚ment en temps
DO
dtBase# = dtBase# / 2
dt# = dtBase#
t# = t0#
x# = x0#
v# = v0#
DO
' pour arriver exactement sur tLimit
IF (t# + dt# > tLimit#) THEN
   dt# = tLimit# - t#
   t# = tLimit#
ELSE
   t# = t# + dt#
END IF

' F(x, v) = (v, GM/x^2)
' ‚tages de Runge-Kutta
xk1# = dt# * v#
vk1# = dt# * (GM# / x# ^ 2)
xg1# = x# + xk1# / 2
vg1# = v# + vk1# / 2
xk2# = dt# * vg1#
vk2# = dt# * (GM# / xg1# ^ 2)
xg2# = x# + xk2# / 2
vg2# = v# + vk2# / 2
xk3# = dt# * vg2#
vk3# = dt# * (GM# / xg2# ^ 2)
xg3# = x# + xk3#
vg3# = v# + vk3#
xk4# = dt# * vg3#
vk4# = dt# * (GM# / xg3# ^ 2)
x# = x# + (xk1# + 2 * xk2# + 2 * xk3# + xk4#) / 6
v# = v# + (vk1# + 2 * vk2# + 2 * vk3# + vk4#) / 6
LOOP UNTIL t# >= tLimit#
PRINT USING "  dt= #####.#####  err= #.###^^^^  x= #######.#########  v= ######.#########"; dtBase#; ABS(x# / 1000 - 9049.721370064#); x# / 1000; v#
LOOP UNTIL dtBase# < .1

END SUB

SUB Runge
' M‚thode de Runge d'ordre 2
GM# = -6.67E-11 * 5.97E+24
' conditions initiales
x0# = 1E+08
v0# = 0
t0# = 0
tLimit# = 55000  ' on s'int‚resse … la position en tLimit

PRINT
PRINT "  M‚thode de Runge d'ordre 2"

dtBase# = 2 ^ 11' incr‚ment en temps
DO
dtBase# = dtBase# / 2
dt# = dtBase#
t# = t0#
x# = x0#
v# = v0#
DO
' assure de calculer jusqu'a exactement  tLimit.
IF (t# + dt# > tLimit#) THEN
  dt# = tLimit# - t#
  t# = tLimit#
ELSE
  t# = t# + dt#
END IF

' premier demi pas de Runge
x1# = x# + v# * dt# / 2
v1# = v# + (GM# / x# ^ 2) * dt# / 2
' pas complet, utilisant le demi-pas pr‚c‚dent.
x# = x# + v1# * dt#
v# = v# + (GM# / x1# ^ 2) * dt#
LOOP UNTIL t# >= tLimit#
PRINT USING "  dt= #####.#####  err= #.###^^^^  x= #######.#########  v= ######.#########"; dtBase#; ABS(x# / 1000 - 9049.721370064#); x# / 1000; v#
LOOP UNTIL dtBase# < .03

END SUB