[<<][meta][>>][..]
Sat Mar 9 11:41:09 EST 2013

Approach for linear exponentials

See [1].  Evaluate 2 exponentials:

exp(P_k) -> reset p_k
exp((P_k_+1-P_k)/n) -> q_k

p_k+1 = p_k q_k^n

Approx error is then:

exp(P_k+1) - p_k+1

Note that the limiting factor is the accuracy of exp(P_k).


Find out:
- what the pure numerical error is with exact exp()
- what the full error is with approx exp()




(define (run P0 P1 n [exp exp])
  (let* ((q  (exp (/ (- P1 P0) n)))
         (p0 (exp P0))
         (p1 (for/fold ((p p0))
                       ((i (in-range n)))
               (* p q))))
    (values p1 (exp P1))))

box> (run +1i +.9i 100)
0.6216099682706652+0.7833269096274834i
0.6216099682706644+0.7833269096274834i

From this it seems that numerical error due to accumulation can
probably be largely ignored.  It pretty much stays in the lower
significant bits.  The real problem is in the approximation of exp.

We need something that satisfies

   exp ((P1 - P2 / n)) ^ n  == exp(P1) / exp(P0) + e

up to an acceptable e.

(define (error P0 P1 n)
  (lambda (exp)
    (- (expt (exp (/ (- P1 P0) n)) n)
       (/ (exp P1) (exp P0)))))

box> (define e (error +1i +.9i 100))
box> (e exp)
3.774758283725532e-15-4.579669976578771e-16i
;; Test  
(define e (error +1i +.9i 100))



;; exp Taylor series
(define (exp-taylor order [acc 0])
  (lambda (x)
    (call-with-values
        (lambda ()
          (for/fold ((acc acc)
                     (xk 1)
                     (kfac 1))
              ((k (in-range (add1 order))))
            (values
             (+ acc (/ xk kfac))
             (* xk x)
             (* kfac (add1 k)))))
      (lambda (acc . _) acc))))



box> (e (exp-taylor 4))
-0.0029735148625679164-0.0017139087382160162i
box> (e (exp-taylor 6))
9.001294569177531e-05+5.099546987226422e-05i
box> (e (exp-taylor 8))
-1.4694914934887393e-06-8.294693956273358e-07i

That isn't too bad.

[1] entry://../math/20130309-104518




[Reply][About]
[<<][meta][>>][..]