;; Module for testing Runge-Kutta integration of ODEs
;; Based on code by  Uwe Hollerbach <uh@alumni.caltech.edu>

(use runge-kutta)

(define (scaler x y)  (fp* x y))
(define (summer x y)  (fp+ x y))

;; Solve the test problem dy/dt = -t^3 y^3 
;;
;;   which has the exact solution y = 1/sqrt(C + t^4/2)
;;   where C = 1/y_0^2 - t_0^4/2

(define con -0.4)
(define (deriv t y) (fp* con y))
(define t0  0.0)
(define y0  1.75)
(define (exact t) (fp* y0 (exp (fp* con (fp- t t0)))))


(define (printst t y) (printf "~A ~A~%" y (fp- y (exact t))))

(define (soln1 stepr t st)
  (let ((nstate (stepr t st)))
    (let ((tn (car nstate))
	  (stn (cadr nstate)))
      (if (fp>= t 5.0) 
	  (printst t st)
	  (begin
	    (soln1 stepr tn stn))))))

(define (soln2 stepr t st)
  (let ((nstate (stepr t st)))
    (let ((tn (car nstate))
	  (stn (cadr nstate)))
      (if (fp>= t 5.0) 
	  (printst t st)
	  (begin
	    (soln2 stepr tn stn))))))

(define (case1 stepr)
  (lambda (n)
    (let ((h (if (negative? n) (expt 2.0 (abs n)) (expt 0.5 n))))
      (printf "~A " h)
      (soln1 (stepr h) t0 y0))))

(define (case2 stepr)
  (lambda (n)
    (let ((h (if (negative? n) (expt 2.0 (abs n)) (expt 0.5 n))))
      (printf "~A " h)
      (soln2 (stepr h) t0 y0))))

(define (solver1 stepr)
  (print "# step yf err")
  (for-each (case1 (stepr scaler summer deriv))
	    (list-tabulate 15 (lambda (x) (- x 2))))
  (print "All done!"))

(define (solver2 stepr)
  (print "# step yf err")
  (for-each (case2 (stepr scaler summer deriv))
	    (list-tabulate 15 (lambda (x) (- x 2))))
  (print "All done!"))


(print "#### Non-Adaptive Solvers")
(print "Forward Euler")
(solver1 (make-rkfe))
(print "Kutta's third order method")
(solver1 (make-rk3))
(print "Classic fourth order method")
(solver1 (make-rk4a))
(print "Kutta's other classic fourth order method")
(solver1 (make-rk4b))
(print "#### Adaptive Solvers")
;;(print "Heun-Euler, order 2/1")
;;(solver2 (make-he))
(print "Bogacki-Shampine, order 3/2")
(solver2 (make-bs))
(print "Runge-Kutta-Fehlberg, order 4/5")
(solver2 (make-rkf45))