BashForth で Adams-Bashforth 法

いろいろと手抜きですが、BashForth を使って、 を時間刻み δt=0.01、2次の Adams-Bashforth 法で解いてみました。
図は解析解と比較したグラフです。

プログラム adams.fs

\ adams.fs
\ fixed point real
variable fpprec
10000 fpprec !
: i>fr
  fpprec @ *
;
: fr.
  dup 0< if
    45 emit
    negate
  then
  fpprec @ /mod
  . 46 emit
  5 0 do
    10 * fpprec @ /mod 48 + emit
  loop
  drop
;
: fr*
  fpprec @ */
;
: fr/
  fpprec @ swap */
;
\ dx/dt=-x
: func ( x -- dx/dt )
  negate
;
variable dt
: euler-method ( x[n] -- fx[n] x[n+1] )
  dup >r func
  dup dt @ fr* r> +
;
: adams-bashforth-2nd ( fx[n-1] x[n] -- fx[n] x[n+1] )
  dup >r func
  \ fx[n-1] fx[n]  | x[n]
  tuck 3 i>fr fr* swap -1 i>fr fr* + 2 i>fr fr/
  \ fx[n] (fx[n]*3+fx[n-1]*-1)/2  | x[n]
  dt @ fr* r> +
  \ fx[n] x[n+1]
;
: solve ( x[0] num -- )
  >r
  dup fr. cr
  euler-method
  dup fr. cr
  r> 0 do
    adams-bashforth-2nd
    dup fr. cr
  loop
;
\ set dt
1 i>fr 100 i>fr fr/ dt !
\ x[0]=1.0
1 i>fr
100 solve
bye

結果

% ./bashforth-0.58a adams.fs | sed 's/ //'
1.00000
0.99000
0.98020
0.97050
0.96090
0.95140
0.94200
0.93270
0.92350
0.91440
0.90540
0.89640
0.88750
0.87870
0.87000
0.86140
0.85290
0.84450
0.83610
0.82780
0.81960
0.81150
0.80350
0.79560
0.78770
0.77990
0.77220
0.76460
0.75700
0.74950
0.74210
0.73480
0.72750
0.72030
0.71320
0.70620
0.69920
0.69230
0.68550
0.67870
0.67200
0.66540
0.65880
0.65230
0.64590
0.63950
0.63320
0.62690
0.62070
0.61460
0.60850
0.60250
0.59660
0.59070
0.58490
0.57910
0.57340
0.56770
0.56210
0.55660
0.55110
0.54570
0.54030
0.53500
0.52970
0.52450
0.51930
0.51420
0.50910
0.50410
0.49910
0.49420
0.48930
0.48450
0.47970
0.47500
0.47030
0.46570
0.46110
0.45660
0.45210
0.44770
0.44330
0.43890
0.43460
0.43030
0.42610
0.42190
0.41780
0.41370
0.40960
0.40560
0.40160
0.39770
0.39380
0.38990
0.38610
0.38230
0.37850
0.37480
0.37110
0.36750