0

I am trying to use PPVAL function from the IMSL library, to evaluate a piecewise polynomial, similar to MATLAB's ppval function.

Sadly, I could not understand well the documentation on how to properly use the function. I tried to evaluate a simple polynomial. The polynomials p(x) = 1 and p(x) = x were ok. What I don't understand is when I try to get higher degree polynomials.

The program below tries to evaluate the polynomial p(x)=x^4:

program main
    include 'LINK_FNL_SHARED.h'
    USE PPVAL_INT
    implicit none

    double precision :: x(11), y(11)
    double precision :: BREAK(2),PPCOEF(5,1)
    integer :: ii

    break = [0.0d0,10.0d0] ! breakpoints in increasing order
    PPCOEF(:,1) = [0.0d0,0.0d0,0.0d0,0.0d0,1.0d0] ! polynomial coefficients
    x = dble([0,1,2,3,4,5,6,7,8,9,10])

    do ii=1,11
        y(ii) = PPVAL(x(ii),break,PPCOEF)
    end do

    print *, 'x = ', x
    print *, 'y = ', y
    pause
end program main

But the function returns the polynomial p(x) = x^4/4!. For degrees 2,3 and 5, the return is always, p(x) = x^2/2!, p(x)=x^3/3!, p(x)=x^5/5!. Why does this factor appear in the ppval function? Why can't I just supply the polynomial coefficients and evaluate it?

Is there another simpler function to evaluate polynomials in Fortran, like MATLAB polyval?

Thales
  • 1,181
  • 9
  • 10
  • This is on purpose. See the definition of PPDER https://docs.roguewave.com/imsl/fortran/6.0/math/default.htm?turl=ppder.htm and replace derivative order (j) by zero... – aka.nice Nov 25 '19 at 10:05
  • `y(ii) = PPDER(x(ii),BREAK,PPCOEF, IDERIV=0)`. Still same behaviour. If I want a quadratic polynomial, `y=x^2`, `ppder` will return `y=x^2/2!`, and so on. Is there a workaround? I would like to simply give the coefficients of the polynomial and the function would return the polynomial evaluated at the x point. – Thales Nov 25 '19 at 10:46
  • Yes, I did not give you a solution, just a confirmation that the behavior you observed was documented. The simplest thing that could possibly work is to multiply coeff[ i ] by factorial( i - 1 ), assuming 1-based coeffs like in your example. – aka.nice Nov 25 '19 at 12:19

0 Answers0