numerical-routines/thiele.lisp

98 lines
4.1 KiB
Common Lisp

; Copyright (c) 2022 Thomas Baruchel
;
; Permission is hereby granted, free of charge, to any person obtaining a copy
; of this software and associated documentation files (the "Software"), to deal
; in the Software without restriction, including without limitation the rights
; to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
; copies of the Software, and to permit persons to whom the Software is
; furnished to do so, subject to the following conditions:
;
; The above copyright notice and this permission notice shall be included in
; all copies or substantial portions of the Software.
;
; THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
; IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
; FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
; AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
; LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
; OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
; SOFTWARE.
;
; Installation
; ------------
; The functions can be used with or without compilation:
; * without compilation:
; load("thiele.lisp")$
; * with compilation (must be compiled only once):
; :lisp (compile-file "thiele.lisp");
; look for the compiled file like "thiele.o" or "thiele.fasl"
; and from now on:
; load("thiele.fasl")$
; Thiele's interpolation formula
; ==============================
;
; In the following example, we use Thiele's interpolation formula
; for computing the inverse cos function cos^(-1) from 32 points
; given by the cos function itself.
;
; We create two vectors:
; * the vector 'a' contains regularly spaced points in [0..1.55]
; * the vector 'b' contains cos(x) for each x in vector 'a'
;
; We set the keepfloat flag to true in order to prevent Maxima
; to replace the values with rational approximations.
;
; Then we compute the inverse function by using thiele(b,a);
; we evaluate this interpolation at 0.5 expecting some value
; close to pi/3.
;
; (%i1) a:makelist(0.05*i,i,0,31)$
; (%i2) b:makelist(float(cos(0.05*i)),i,0,31)$
; (%i3) keepfloat:true$
; (%i4) subst(x=0.5, thiele(b,a))*3;
; (%o4) 3.14159265357928
;
; The function may be used with exact values, but it is highly sensitive
; to input values, and equally spaced abscissæ often yield some division
; by zero error.
;
; (%i5) thiele([1,2,5,6], [10,12,11,13]);
; 2
; 9 x + 29 x - 238
; (%o5) -----------------
; 8 x - 28
(defun $thiele (u v)
(loop for rh on
(loop for w on (cdr u) ; we build rho with axes transposed (and reversed)
; first, building column 2 and 1
for r = (list
(loop for a on u
for b on v
while (cdr a)
collect (div (sub (car a) (cadr a))
(sub (car b) (cadr b)))) v)
; then building starting from column 3
then (cons
(loop for a1 in u
for a2 in w
for y1 on (car r)
for y2 in (cdr (cadr r))
collect (add (div (sub a1 a2)
(sub (car y1) (cadr y1)))
y2)) r)
; finally get the row 1
; finally (return (loop for i in r ; i is already a 1-list
; for j = i then (cons (car i) j)
; finally (return j))))
finally (return (nreverse (mapcar #'car r))))
for k in (cdr u)
with a = 0
do (setf a (add a (div (sub 'x k)
(sub (caddr rh) (car rh)))))
finally (return (add (car v)
(div (sub 'x (car u))
(sub (cadr rh) a))))))