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 x-coordinates of the function, ys is an array of length n representing the y-coordinates 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 = 1e-25 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

  1. The ratint function takes three arguments: xs, ys, and x.

  2. The length of the input arrays xs and ys are obtained and verified to be equal. If they are not equal, an error is thrown.

  3. Two arrays c and d are created as copies of the ys array.

  4. A variable hh is initialized to the absolute difference between x and the first element of xs. This will be used to keep track of the smallest distance between x and an element of xs.

  5. Three variables y, dy, and ns are initialized to 0. y will be used to store the approximation of the function at x, dy will be used to store the derivative of the approximation at x, and ns will be used to keep track of the index of the closest element of xs to x.

  6. A variable eps is initialized to a small number. This will be used to prevent division by zero.

  7. A for loop is used to iterate over the xs array.

  8. For each xs[i], the absolute difference between x and xs[i] is computed and stored in a variable h.

  9. If h is equal to 0, then ys[i] is the exact value of the function at x, so y is set to ys[i], dy is set to 0, and an exception is raised to exit the loop.

  10. If h is less than hh, then ns is set to i, hh is set to h, and the c[i] and d[i] arrays are updated to ys[i]. This keeps track of the closest element of xs to x.

  11. After the loop, y is set to ys[ns], and ns is decremented by 1.

  12. Two nested for loops are used to iterate over the c and d arrays.

  13. For each iteration, m represents the current iteration of the outer loop, and i represents the current iteration of the inner loop.

  14. A variable w is computed as the difference between c[i] and d[i-1].

  15. A variable h is computed as the difference between xs[i+m-1] and x.

  16. A variable t is computed as (xs[i-1] - x) * d[i-1] / h.

  17. A variable dd is computed as t - c[i].

  18. If dd is equal to 0, then the function has a pole, and an exception is thrown.

  19. A variable dd is computed as w / dd.

  20. The d[i-1] array is updated to c[i] * dd, and the c[i-1] array is updated to t * dd.

  21. After the loops, y and dy are returned as the approximation and derivative of the function at x.

  22. 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 n-1 times, and the inner loop iterates n-m times, where m is the current iteration of the outer loop. This gives a total of (n-1) + (n-2) + ... + 1 iterations for the outer loop, which is equal to n(n-1)/2. The inner loop has a similar number of iterations, giving a total of n(n-1)^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.