CODEWITHSUNDEEP
fortrannumerical-methodsfortran90lapackpolynomials
maciek
maciek

Reputation: 2117

Insufficient numerical precision of LAPACK procedures (Fortran 90)?

I have written a simple Fortran code which carries out polynomial interpolation of n+1 points of R^2. It solves a system of linear equations (I am creating Vandermonde matrix) with two LAPACK procedures (everything is double precision in my code):
At first, it factorises the matrix: http://sites.science.oregonstate.edu/~landaur/nacphy/lapack/routines/dgetrf.html
Later, it solves the system: http://sites.science.oregonstate.edu/~landaur/nacphy/lapack/routines/dgetrs.html

the program works fine for a few tested cases of data produced from polynomials of degrees: 0,1,2,3,4. However, when I provide 11 points from a polynomial P(x) = x^10 it infers completely wrong coefficients.

Input (x..., y...):

1.0
1.001
1.002
1.003
1.004
1.005
1.006
1.007
1.008
1.009
1.01
1.0
1.01004512021
1.02018096336
1.03040825707
1.04072773401
1.05114013204
1.06164619412
1.07224666847
1.08294230847
1.0937338728
1.10462212541

Output (a_n,...,a_0):

-4.6992230177E+004
2.2042918738E+005
-3.2949938635E+005
5.0740528522E+004
2.4764893257E+005
-3.1846974845E+004
-1.7195378863E+005
-1.4751075818E+005
4.1766709741E+005
-2.6476448046E+005
5.6082872757E+004

Am I hitting numerical stability problems? Or did I do something wrong?

Edit: I attach the code for the interpolation procedure (watch out, we actually have n points not n+1).

module InterpolatorModule
contains
subroutine interpolate(n, x, y, a)
  implicit none
  integer :: n, i, j, info
  integer, dimension (:), allocatable :: ipiv
  real(8), dimension (:), pointer :: x, y, a
  real(8), dimension(:,:), allocatable :: Mat_X, Mat_B

  ! create the Vandermonde matrix:
  allocate ( Mat_X(n,n) )
  do i = 1,n
    do j = 1,n
      Mat_X(i,j) = x(i) ** ( n - j )
    end do
  end do

  ! reshape ipnut data into a matrix form:
  allocate ( Mat_B(n,1) )
  do i = 1,n
    Mat_B(i,1) = y(i)
  end do

  ! prepare an array for the pivot indices from DGETRF:
  allocate(ipiv(n))

  call DGETRF(n, n, Mat_X, n, ipiv, info)
  call DGETRS("N", n, 1, Mat_X, n, ipiv, Mat_B, n, info)

  ! save the results into the argument array
  do i = 1,n
    a(i) = Mat_B(i,1)
  end do

  deallocate(ipiv)
  deallocate ( Mat_X )
  deallocate ( Mat_B )

end subroutine interpolate
end module InterpolatorModule

Edit: the main program:

program MyInterpolation
use InterpolatorModule
implicit none

  integer :: line, n
  real(8), dimension (:), pointer :: p_x, p_y, p_a
  real(8), dimension(:), target, allocatable :: in1_x, in1_y, in1_a 
  character(len=23) :: str

  open (1, file="test/in.txt", status="old", action="read", form="formatted")

    n = 11 ! this value is not known at compulation time, simplification for SO

    ! allocate memory for the input data
    allocate ( in1_x(n) )
    allocate ( in1_y(n) )
    allocate ( in1_a(n) )

    ! read the x_i coordinates:
    do line = 1,n
      read(1,*) in1_x(line)
    end do
    ! read the y_i coordinates:
    do line = 1,n
      read(1,*) in1_y(line)
    end do
  close(1)
  ! assign pointers to the arrays:
  p_x => in1_x
  p_y => in1_y
  p_a => in1_a
  ! call the interpolating procedure:
  call interpolate(n, p_x, p_y, p_a)

  ! print out calculated coefficients:
  do line = 1,n
    write(str,'(ES23.10 E3)') in1_a(line)
    write(*,'(a)') adjustl(trim(str))
  end do

  ! free the allocated memory
  deallocate (in1_x)
  deallocate (in1_y)
  deallocate (in1_a)

  ! ===========================================================================

end program MyInterpolation

Upvotes: 2

Views: 434

Answers (1)

njuffa
njuffa

Reputation: 26185

It is widely recognized – and has been pointed out in comments – that the Vandermonde matrix is frequently ill-conditioned. For example a standard textbook on numerical analysis, D. Kincaid and W. Cheney, "Numerical Analysis, 2nd ed.", Brooks/Cole Publishing 1996, states on page 336:

The coefficient matrix in Equation (12) is called the Vandermonde matrix. It is non-singular because the system has a unique solution for any choice of y0, y1, ..., yn [...]. The determinant of the Vandermonde matrix is therefore nonzero for distinct nodes x0, x1, ..., xn. [...] However, the Vandermonde matrix is often ill conditioned, and the coefficients ai may therefore be inaccurately determined by solving System (12). (See Gautschi [1984]). [...] Therefore, this approach is not recommended.

Inspection of the matrix returned by GETRF already informs us that we are in trouble for polynomials of fairly low degrees. At degree 4 the magnitude of the smallest relevant element after LU decomposition is on the order of 10-12, and at degree 5 the smallest such element is 10-14. Based on the fact that all elements of the original matrix were near unity, and knowing the basic steps of the LU decomposition process, it is clear that the tiny elements in the GETRF result come about by subtractive cancellation.

The magnitude of the elements getting close to the double precision epsilon (≈ 10-16) tells us that most of the accuracy has been lost. For polynomials of degrees 6 and higher the code is basically operating on pure noise.

We can further confirm this by comparing the computed coefficients of the interpolation polynomial with a more robust reference. For simplicity I used Wolfram Alpha for this comparison. For a degree 4 polynomial, the Fortran code computes the coefficients accurate to about three decimal digits, and for a polynomial of degree 5, this shrinks to a single correct decimal digit.

In terms of a simple and more robust algorithm to generate the coefficients of the interpolation polynomial, I found the following:

J. N. Lyness and C.B. Moler, "Van Der Monde Systems and Numerical Differentiation." Numerische Mathematik 8, 458-464 (1966)

I translated the Algol code in the paper to ISO-C99 and it seems to deliver reasonable results up to degree 8 by comparison with Wolfram Alpha. Wolfram Alpha gives up for degrees greater than 8, and I do not have another reference handy. Even with this more robust algorithm there seems to be significant loss of precision for higher degrees, with only 6 decimal digits matching between Wolfram Alpha and the Luness/Moler algorithm.

#include <stdio.h>
#include <stdlib.h>

/* J. N. Lyness and C.B. Moler, "Van Der Monde Systems and Numerical 
   Differentiation." Numerische Mathematik 8, 458-464 (1966)
*/
void update (int k, int p, double *x, double fxk, double *C)
{
    int d, s, m, n;
    double xk, xkd;
    xk = x[k];
    m = k*(k+3)/2;
    C[m] = fxk;
    for (d = 1; d <= k; d++) {
        xkd = xk - x[k-d];
        for (s = 0; s <= ((d > p) ? p : d); s++) {
            m = m - 1;
            n = m + d;
            if (s == 0) {
                C[m] = C[n] + xk * (C[n-k-1] - C[n]) / xkd;
            } else if (s == d) {
                C[m] = (C[n+1] - C[n-k]) / xkd;
            } else {
                C[m] = C[n] + (xk * (C[n-k-1] - C[n]) + (C[n+1] - C[n-k]))/ xkd;
            }
            if (d > p) {
                m = m - d + p;
            }
        }
    }        
}

/* Solve (k+1) x (k+1) system Vc = f, where V[i,j] = x[i]**j  */
void vandal (int k, double *x, double *f, double *c)
{
    double C [k*(k+3)/2+1];
    for (int i = 0; i <= k; i++) {
        update (i, k, x, f[i], C);
    }
    for (int i = 0; i <= k; i++) {
        c[i] = C[k-i];
    }
}

#define k 10
int main (void)
{
    double x[k+1] = {1.000, 1.001, 1.002, 1.003, 1.004, 1.005, 
                     1.006, 1.007, 1.008, 1.009, 1.010};
    double f[k+1] = {1.00000000000, 1.01004512021, 1.02018096336, 
                     1.03040825707, 1.04072773401, 1.05114013204, 
                     1.06164619412, 1.07224666847, 1.08294230847, 
                     1.09373387280, 1.10462212541};
    double c[k+1] = {0};
    
    vandal (k, x, f, c);

    for (int i = 0; i <= k; i++) {
        printf ("c[%2d] = % 23.16e\n", i, c[i]);
    }
    return EXIT_SUCCESS;
}

Upvotes: 5

Related Questions