Search code examples
schemeelispsicp

Riemann integral formula to compute high-order function


SICP introduced Riemann integral formula in Chapter 1.3.1

(define (integral f a b dx)
  (define (add-dx x) (+ x dx))
  (* (sum f (+ a (/ dx 2.0)) add-dx b)
     dx))

Apply it to a particular case

#+name: case-1.3.1-integral.scm
#+BEGIN_SRC scheme :session sicp
(define pi 3.141592653589793)

(define (integral2 f a b dx)
  (define (add-dx x) (+ x dx))
  (* (sum (f b)
          (+ a (/ dx 2.0))
          (lambda (x) (+ x dx))
          b)
     dx))

(define (f b)
  (lambda (x) (/ 1 (sqrt
                    (- (sin x)
                       (sin b))))))


(* (integral2 f 0 (/ pi 6) 0.00001)
   (sqrt (/ 40
            (* 3 9.8))))
#+END_SRC

#+RESULTS: case-1.3.1-integral.scm
: 0.0-1.777598336021436i

Got a perfect answer 1.777598336021436

Then translate it to elisp

Start from small:

#+name: case-1.3.1-integral.el
#+begin_src emacs-lisp :session sicp :lexical t
(defun integral (f a b dx)
  (* (sum f
     (+ a (/ dx 2.0))
     (lambda (x) (+ x dx))
     b)
     dx))

(defun sum(term a next b)
  (if (> a b)
      0
    (+ (funcall term a)
       (sum term (funcall next a) next b))))

(integral #'cube 0 1 0.01)
#+end_src

#+RESULTS: case-1.3.1-integral.el
: 0.24998750000000042

It works and thus use it to solve the previous problem

#+begin_src emacs-lisp :session sicp :lexical t 
(defvar pi 3.141592653589793)

(defun integral (f a b dx)
  (* (sum f
     (+ a (/ dx 2.0))
     (lambda (x) (+ x dx))
     b)
     dx))

(defun f (b)
  (lambda (x) (/ 1 (sqrt
                    (- (sin x)
                       (sin b))))))

(defun integral2 (f a b dx)
  (* (sum (funcall f b)
          (+ a (/ dx 2.0))
          (lambda (x) (+ x dx))
          b)
      dx))

(integral2 #'f 0 (/ pi 6) 0.01)

#+end_src

But it return a meaningless result

ELISP> (integral2 #'f 0 (/ pi 6) 0.01)
-0.0e+NaN

What's the problem?


Solution

  • The answer you obtained when using Scheme is a complex number, the result of calling sqrt (are you sure the Scheme code was correct in the first place? you should double-check it):

    0.0-1.777598336021436i
    

    Unfortunately, Elisp doesn't support complex numbers, that's why we get a NaN in there. But that's not the real problem; you should investigate why are you getting complex results in the Scheme code, an integral should not return complex values!