Rational Interpolation
Introduction
Rational interpolation is a numerical method used to approximate a function using a rational function. The algorithm is used in various fields including engineering, physics, and finance.
Implementation
The ratint
algorithm is implemented in OCaml as shown below. The function takes three arguments: xs
, ys
, and x
. xs
is an array of length n
representing the xcoordinates of the function, ys
is an array of length n
representing the ycoordinates of the function, and x
is the point at which the function is to be approximated. The function returns a tuple containing the approximation of the function at x
and the derivative of the approximation at x
.
let ratint xs ys x =
let n = Array.length xs in
let m = Array.length ys in
let error () =
let s =
Printf.sprintf
"ratint requires that xs and ys have the same length, but xs length is %i \
whereas ys length is %i"
n
m
in
Owl_exception.INVALID_ARGUMENT s
in
Owl_exception.verify (m = n) error;
let c = Array.copy ys in
let d = Array.copy ys in
let hh = ref (abs_float (x . xs.(0))) in
let y = ref 0. in
let dy = ref 0. in
let ns = ref 0 in
let eps = 1e25 in
try
for i = 0 to n do
let h = abs_float (x . xs.(i)) in
if h = 0.
then (
y := ys.(i);
dy := 0.;
raise Owl_exception.FOUND)
else if h < !hh
then (
ns := i;
hh := h;
c.(i) < ys.(i);
d.(i) < ys.(i) +. eps)
done;
y := ys.(!ns);
ns := !ns  1;
for m = 1 to n  1 do
for i = 1 to n  m do
let w = c.(i) . d.(i  1) in
let h = xs.(i + m  1) . x in
let t = (xs.(i  1) . x) *. d.(i  1) /. h in
let dd = t . c.(i) in
if dd = 0. then failwith "Has a pole";
let dd = w /. dd in
d.(i  1) < c.(i) *. dd;
c.(i  1) < t *. dd
done
done;
!y, !dy
with
 Owl_exception.FOUND > !y, !dy
 e > raise e

The
ratint
function takes three arguments:xs
,ys
, andx
. 
The length of the input arrays
xs
andys
are obtained and verified to be equal. If they are not equal, an error is thrown. 
Two arrays
c
andd
are created as copies of theys
array. 
A variable
hh
is initialized to the absolute difference betweenx
and the first element ofxs
. This will be used to keep track of the smallest distance betweenx
and an element ofxs
. 
Three variables
y
,dy
, andns
are initialized to 0.y
will be used to store the approximation of the function atx
,dy
will be used to store the derivative of the approximation atx
, andns
will be used to keep track of the index of the closest element ofxs
tox
. 
A variable
eps
is initialized to a small number. This will be used to prevent division by zero. 
A for loop is used to iterate over the
xs
array. 
For each
xs[i]
, the absolute difference betweenx
andxs[i]
is computed and stored in a variableh
. 
If
h
is equal to 0, thenys[i]
is the exact value of the function atx
, soy
is set toys[i]
,dy
is set to 0, and an exception is raised to exit the loop. 
If
h
is less thanhh
, thenns
is set toi
,hh
is set toh
, and thec[i]
andd[i]
arrays are updated toys[i]
. This keeps track of the closest element ofxs
tox
. 
After the loop,
y
is set toys[ns]
, andns
is decremented by 1. 
Two nested for loops are used to iterate over the
c
andd
arrays. 
For each iteration,
m
represents the current iteration of the outer loop, andi
represents the current iteration of the inner loop. 
A variable
w
is computed as the difference betweenc[i]
andd[i1]
. 
A variable
h
is computed as the difference betweenxs[i+m1]
andx
. 
A variable
t
is computed as(xs[i1]  x) * d[i1] / h
. 
A variable
dd
is computed ast  c[i]
. 
If
dd
is equal to 0, then the function has a pole, and an exception is thrown. 
A variable
dd
is computed asw / dd
. 
The
d[i1]
array is updated toc[i] * dd
, and thec[i1]
array is updated tot * dd
. 
After the loops,
y
anddy
are returned as the approximation and derivative of the function atx
. 
If any exceptions are raised during the computation, they are rethrown.
Complexity Analysis
The ratint
algorithm has a time complexity of O(n^2), where n is the length of the input arrays xs
and ys
.
This is because the algorithm has two nested loops, each iterating over the c
and d
arrays. The outer loop iterates n1
times, and the inner loop iterates nm
times, where m
is the current iteration of the outer loop. This gives a total of (n1) + (n2) + ... + 1
iterations for the outer loop, which is equal to n(n1)/2
. The inner loop has a similar number of iterations, giving a total of n(n1)^2/2
iterations for both loops combined.
Therefore, the time complexity of the algorithm is O(n^2). This means that as the length of the input arrays xs
and ys
increases, the time taken by the algorithm to compute the approximation of the function at x
and its derivative also increases quadratically.