; by Alex Klinkhamer  (:handle grencez)
; Update: 2008.07.12
;
;  Legal info:
; You may exploit this code,
;    and it may exploit you.
; I am not responsible for either.


; The stuff in this file is my overindulgence in a differential equations
; homework assignment which compares the different methods for approximating
; integrals given the derivative's formula.
;
; These don't actually give the definite integrals, only the integrated
; function value at a certain x... so the actual integral is just a subtraction...
; Integral of y'(x,y) from x0 to x1 is y(x1) - y(x0)
; or in lisp form... (- (runge-kutta #'y-prime x1 0.05d0 x0 y0) y0)
; where we use the runge-kutta approximation method with 0.05 stepsize
; and y-prime is y'(x,y)

(in-package :cl-user)

; Set to nil if you don't want output.
; But if you're using this for efficiency,
; just remove all the format stuff.
;
; The examples specifically bind it.
(defvar *show-iterations* 5)

(defmacro defapprox (name (y-prime x-find h x0 y0) . body)
  `(defun ,name (,y-prime ,x-find ,h ,x0 ,y0)
     (declare (type double-float ,x-find)
              (type double-float ,h)
              (type double-float ,x0)
              (type double-float ,y0))
     (when *show-iterations*
       (format t "~%~%~A~%" (quote ,name)))
     (when (zerop ,h) (return-from ,name ,y0))
     (when *show-iterations*
       ,(car body))
     ,@(cdr body)))

;V The simplest, least effective approximation.V
(defapprox
  euler-approx (y-prime x-find h x0 y0)
  (format t "    ~9,a ~9,a ~9,a ~9,a~%" 'x 'y 'k 'y_n+1)
  (loop :with iters = 1
        :for y = y0 :then (+ y (* h k))
        :for x = x0 :then (+ x h)
        :for k = (funcall y-prime x y)
        :repeat (round (/ (- x-find x0) h))
        :when (and *show-iterations* (zerop (decf iters)))
        :do (format t "~9,5f ~9,5f ~9,5f ~9,5f~%" x y0 k y)
        :and :do (setq iters *show-iterations*)
        :finally
        (when *show-iterations*
          (format t "~9,5f ~9,5f ~9,5f ~9,5f~%" x y0 k y))
        (return y)))

;V Improved Euler Method.V
(defapprox
  imp-euler-approx (y-prime x-find h x0 y0)
  (format t "    ~9,a ~9,a ~9,a ~9,a ~9,a~%"
          'x 'y 'k1 'k2 'y_n+1)
  (loop :with iters = 1
        :for k1 = (coerce 0 'double-float) 
        :then (funcall y-prime x y)
        :for x = x0 :then (+ x h)
        :for k2 = (coerce 0 'double-float)
        :then (funcall y-prime x (+ y (* h k1)))
        :for y = y0 :then (+ y (/ (* h (+ k1 k2)) 2))
        :repeat (round (/ (- x-find x0) h))
        :when (and *show-iterations* (zerop (decf iters)))
        :do (format t "~9,5f ~9,5f ~9,5f~%~
                    ~9,5f ~9,5f " k1 k2 y x y)
        :and :do (setq iters *show-iterations*)
        :finally
        (when *show-iterations*
          (format t "~9,5f ~9,5f ~9,5f~%~
                  ~9,5f ~9,5f " k1 k2 y x y))
        (return y)))

;V Best method.V
(defapprox
  runge-kutta (y-prime x-find h x0 y0)
  (format t "    ~9,a ~9,a ~9,a ~9,a ~9,a ~9,a ~9,a~%"
          'x 'y 'k1 'k2 'k3 'k4 'y_n+1)
  (loop :with half-h = (/ h 2)
        :with iters = 1
        :for x = x0 :then (+ x h)
        :for k4 = 0d0
        :then (funcall y-prime x (+ y (* h k3)))
        :for y = y0 
        :then (+ y (* half-h (/ (+ k1 (* 2 k2) (* 2 k3) k4) 3)))
        :for k1 = (funcall y-prime x y)
        :for k2 = (funcall y-prime (+ x half-h)
                           (+ y (* half-h k1)))
        :for k3 = (funcall y-prime (+ x half-h)
                           (+ y (* half-h k2)))
        :repeat (round (/ (- x-find x0) h))
        :when (and *show-iterations* (zerop (decf iters)))
        :do (format t "~9,5f ~9,5f~%~
                    ~9,5f ~9,5f ~9,5f ~9,5f ~9,5f "
                    k4 y
                    x y k1 k2 k3)
        :and :do (setq iters *show-iterations*)
        :finally
        (when *show-iterations*
          (format t "~9,5f ~9,5f~%~
                  ~9,5f ~9,5f ~9,5f ~9,5f ~9,5f "
                  k4 y
                  x y k1 k2 k3))
          (return y)))

(defun showall
  (y-prime x-find h x0 y0
           &optional (exact-eqn nil)
           &aux (exact (when exact-eqn
                         (funcall exact-eqn x-find))))
  (mapc
    (lambda (approx-fn)
      (let ((approx (funcall approx-fn y-prime x-find h x0 y0)))
        (format t "~%y(~3f) = ~f" x-find approx)
        (when exact
          (format t "~%Error = ~f" (/ (- approx exact) exact)))))
    (list #'euler-approx #'imp-euler-approx #'runge-kutta))
  (fresh-line))

;--- BEGIN EXERCISES ---

; #5
; y' = y - x - 1
(defun n5-deriv (x y)
  (- y x 1))
(defun n5 (x)
  (- (+ 2 x) (exp x)))
; y(0) = 1
; y(0.5) = ?
(defun show-all-n5 ()
  (let ((*show-iterations* 1))
    (showall #'n5-deriv 0.5d0 0.25d0 0d0 1d0 #'n5)))

; #9
; y' = (y^2 + 1) / 4
(defun n9-deriv (x y)
  (declare (ignore x))
  (/ (1+ (expt y 2)) 4))
(defun n9 (x)
  (tan (/ (+ x pi) 4)))
; y(0) = 1
; y(0.5) = ?
(defun show-all-n9 ()
  (let ((*show-iterations* 1))
    (showall #'n9-deriv 0.5d0 0.25d0 0d0 1d0 #'n9)))

; #13
; yy' = 2x^3
; so...
; y' = 2x^3 / y
(defun n13-deriv (x y)
  (/ (* 2 (expt x 3)) y))
(defun n13 (x)
  (sqrt (+ (expt x 4) 8)))
; y(1) = 3
; y(2) = ?
(defun show-all-n13 ()
  (let ((*show-iterations* 20))
    (showall #'n13-deriv 2d0 0.01d0 1d0 3d0 #'n13)))

; #15
; xy' = 3x - 2y
; so...
; y' = 3 -  2y / x
(defun n15-deriv (x y)
  (- 3 (/ (* 2 y) x)))
(defun n15 (x)
  (+ x (/ 4 (expt x 2))))
; y(2) = 3
; y(3) = ?
(defun show-all-n15 ()
  (let ((*show-iterations* 40))
    (showall #'n15-deriv 3d0 0.005d0 2d0 3d0 #'n15)))

; y' = y
(defun e-fn-deriv (x y)
  (declare (ignore x))
  y)
; y = e^x
(defun e-fn (x) (exp x))
; y(0) = 1
; y(1) = e = ?
(defun show-all-e-fn ()
  (let ((*show-iterations* 40))
    (showall #'e-fn-deriv 1d0 0.005d0 0d0 1d0 #'e-fn)))

; y' = 1 / x
(defun ln-fn-deriv (x y)
  (declare (ignore y))
  (/ x))
; y = ln(x)
(defun ln-fn (x) (log x))
; y(1) = 0
; y(2) = ln(2) = ?
(defun show-all-ln-fn ()
  (let ((*show-iterations* 4))
    (showall #'ln-fn-deriv 2d0 0.05d0 1d0 0d0 #'ln-fn)))

; y' = 4 / (1 + x^2)
(defun pi-fn-deriv (x y)
  (declare (ignore y))
  (/ 4 (1+ (expt x 2))))
; y = 4 arctan(x)
(defun pi-fn (x) (* 4 (atan x)))
; y(0) = 0
; y(1) = pi = ?
(defun show-all-pi-fn ()
  (let ((*show-iterations* 40))
    (showall #'pi-fn-deriv 1d0 0.005d0 0d0 0d0 #'pi-fn)))

;--- END EXERCISES ---

#-(or)(show-all-n5)
#-(or)(show-all-n9)
#-(or)(show-all-n13)
#-(or)(show-all-n15)
#-(or)(show-all-e-fn)
#-(or)(show-all-ln-fn)
#-(or)(show-all-pi-fn)