<- previous    index    next ->

Lecture 5, polynomials


A polynomial of degree n has the largest exponent value n.
There are n+1 terms and n+1 coefficients.
A normalized polynomial has the coefficient of the
largest exponent equal to 1.0.

An nth order polynomial with coefficients c_0 through c_n.

  y = c_0 + c_1 x + c_2 x^2 +...+ c_n-1 x^n-1 + c_n x^n

Horners method for evaluating a polynomial

"y" is computed numerically using Horners method that needs n multiplications and n additions. Thus faster and more accurate. Starting at the highest power, set y = c_n * x Then combine the next coefficient y = (c_n-1 + y) * x Continue for c_n-2 until c_1 y = (c_1 + y) * x Then finish y = c_0 + y The code in C, Fortran 90, Java and Ada 95 is: /* peval.c Horners method for evaluating a polynomial */ double peval(int n, double x, double c[]) { /* an nth order polynomial has n+1 coefficients */ /* stored in c[0] through c[n] */ /* y = c[0] + c[1]*x + c[2]*x^2 +...+ c[n]*x^n */ int i; double y = c[n]*x; if(n<=0) return c[0]; for(i=n-1; i>0; i--) y = (c[i]+y)*x; return y+c[0]; } /* end peval */ ! peval.f90 Horners method for evaluating a polynomial function peval(n, x, c) result (y) ! an nth order polynomial has n+1 coefficients ! stored in c(0) through c(n) ! y = c(0) + c(1)*x + c(2)*x^2 +...+ c(n)*x^n implicit none integer, intent(in) :: n double precision, intent(in) :: x double precision, dimension(0:n), intent(in) :: c double precision :: y integer i if(n<=0) y=c(0); return y = c(n)*x do i=n-1, 1, -1 y = (c(i)+y)*x end do y = y+c(0) end function peval // peval.java Horners method for evaluating a polynomial double peval(int n, double x, double c[]) { // an nth order polynomial has n+1 coefficients // stored in c[0] through c[n] // y = c[0] + c[1]*x + c[2]*x^2 +...+ c[n]*x^n double y = c[n]*x; if(n<=0) return c[0]; for(int i=n-1; i>0; i--) y = (c[i]+y)*x; return y+c[0]; } // end peval -- peval.adb Horners method for evaluating a polynomial function peval(n : integer; x : long_float; c : vector) return long_float is -- an nth order polynomial has n+1 coefficients -- stored in c(0) through c(n) -- y := c(0) + c(1)*x + c(2)*x^2 +...+ c(n)*x^n y : long_float; begin if n<=0 then return c(0); end if; y := c(n)*x; for i in reverse 1..n-1 loop -- do i:=n-1, 1, -1 y := (c(i)+y)*x; end loop; y := y+c(0); return y; end peval; MatLab y = polyval(c, x) Python from numpy import polyval y = polyval(c, x)

Testing high degree polynomials

Although Horners method is fast and accurate on most polynomials, the following test programs in C, Fortran 90, Java and Ada95 show that evaluating polynomials of order 9 and 10 can have significant absolute error. The test programs generate a set of roots r_0, r_1, ... and computer the coefficients of a set of polynomials y = (x-r_0)*(x-r_1)*...*(x-r_n-1) becomes y = c_0 + c_1*x + c_2*x^2 +...+ c_n-1*x^n-1 + 1.0*x^n when x is any of the roots, y should be zero. Study one of the .out files and see how the absolute error increases with polynomial degree and increases with the magnitude of the root. (The source code was run on various types of computers and various operating systems. Results vary.) peval.c peval_c.out peval.f90 peval_f90.out peval.java peval_java.out peval.adb peval_ada.out test_peval.rb test_peval_rb.out Having the coefficients of a polynomial allows creating the derivative and integral of the polynomial.

Derivative of Polynomial

Given: y = c_0 + c_1 x + c_2 x^2 +...+ c_n-1 x^n-1 + c_n x^n dy/dx = c_1 + 2 c_2 x + 3 c_3 x^2 +...+ n-1 c_n-1 x^n-2 + n c_n x^n-1 The coefficients of dy/dx may be called cd and computed: for(i=0; i<n; i++) cd[i] = (double)(i+1)*c[i+1]; nd = n-1; The derivative may be computed an any point, example x=a, using dy_dx_at_a = peval(nd, a, cd);

Integral of Polynomial

Similarly, the integral y(x) dx has coefficients Given: y = c_0 + c_1 x + c_2 x^2 +...+ c_n-1 x^n-1 + c_n x^n int_y_dx = 0 + c_0 x + c_1/2 x^2 +...+ c_n-1/n x^n + c_n/(n+1) x^n+1 The coefficients of int_y_dy may be called ci and computed: ci[0] = 0.0; /* a reasonable choice */ ni = n+1; for(i=1; i<=ni; i++) ci[i] = c[i-1]/(double)i; The integral of the original polynomial y from a to b is computed: int_y_a_b = peval(ni, b, ci) - peval(ni, a, ci);

Polynomial Arithmetic

The sum, difference, product and quotient of two polynomials p and q are: (unused locations filled with zeros) nsum = max(np,nq); for(i=0; i<=nsum; i++) sum[i] = p[i] + q[i]; ndifference = max(np,nq); for(i=0; i<=ndifference; i++) difference[i] = p[i] - q[i]; nproduct = np+nq; for(i=0; i<=nproduct; i++) product[i] = 0.0; for(i=0; i<=np; i++) for(j=0; j<=nq; j++) product[i+j] += p[i]*q[j]; /* np > nd */ nquotient = np-nd; nremainder = np k = np; for(j=0; j<=np; j++) r[j] = p[j]; /* initial remainder */ for(i=nquotient; i>=0; i--) { quotient[i] = r[k]/d[nd] for(j=nd; j>=0; j--) r[k-nd+j] = r[k-nd+j] - quotient[i]*d[j] k--; }

Polynomial operations including add,sub,mul,div,roots,deriv,integral

poly.h poly.c test_poly.c test_poly_c.out poly.java test_poly.java test_poly_java.out In Python all available import numpy.matlib To see all functions print dir(numpy.matlib)

Polynomial series

Taylor series, for any differentiable function, f(x) (x-a) f'(a) (x-a)^2 f''(a) (x-a)^3 f'''(a) f(x) = f(a) + ----------- + -------------- + --------------- + ... 1! 2! 3! MacLaurin series, Taylor series with a=0 x f'(0) x^2 f''(0) x^3 f'''(0) f(x) = f(0) + ------- + ---------- + ----------- + ... 1! 2! 3! Taylor series, offset h f'(x) h^2 f''(x) h^3 f'''(x) f(x+h) = f(x) + ------- + ---------- + ----------- + ... 1! 2! 3! Example f(x) = e^x, thus f'(x) = e^x f''(x) = e^x f'''(x) = e^x f(0) = 1 f'(0) = 1 f''(0) = 1 f'''(0) = 1 substituting in the MacLaurin series x x^2 x^3 x^4 f(x) = 1 + -- + --- + --- + --- + ... 1! 2! 3! 4! Example f(x) = sin(x), f'(x) = cos(x) f''(x) = -sin(x) f'''(x) = -cos(x) f(0) = 0 f'(0) = 1 f''(0) = 0 f'''(0) = -1 substituting in the MacLaurin series x x^3 x^5 x^7 f(x) = -- - --- + --- - --- + ... 1! 3! 5! 7! Example f(x) = cos(x), f'(x) = -sin(x) f''(x) = -cos(x) f'''(x) = sin(x) f(0) = 1 f'(0) = 0 f''(0) = -1 f'''(0) = 0 substituting in the MacLaurin series x^2 x^4 x^6 f(x) = 1 - --- + --- - --- + ... 2! 4! 6! Many MacLaurin series expansions are shown on Mathworld Taylor and MacLaurin series can be expanded for two and more variables. Two dimensional Taylor series expansion is here

Orthogonal Polynomials

For use in integration and fitting data by a polynomial, orthogonal polynomials are used. The definition of orthogonal is based on having a set of polynomials from degree 0 to degree n that are denoted by p_0(x), p_1(x),...,p_n-1(x), p_n(n). A weighting function is allowed, w(x), that depends on x but not on the degree of the polynomial. The set of polynomials is defined as orthogonal over the interval x=a to x=b when the following two conditions are satisfied: / 0 if i != j integral from x=a to x=b of w(x) p_i(x) p_j(x) dx = | \ not 0 if i=j sum from i=1 to i=n c_i *p_i(x) = 0 for all x only when all c_i are zero. Legendre Polynomials, Chebyshev Polynomials, Laguerre Polynomials and Lagrange Polynomials are covered below.

Legendre Polynomials

P_n(x) are defined as P_0(x) = 1 P_1(x) = x P_2(x) = 1/2 (3 x^2 - 1) P_3(x) = 1/2 (5 x^3 - 3 x) P_4(x) = 1/8 (35 x^4 - 30 x^2 + 3) the general recursion formula for P_n(x) is P_n(x) = 2n-1/n x P_n-1(x) - n-1/n P_n-2(x) The weight function is w(x) = 1 The interval is a = -1 to b = 1 Explicit expression P_n(x) = 1/2^n sum m=0 to n/2 (-1)^m (2n-2m)! / [(n-m)! m! (n-2m)!] x^(n-2m) P_n(x) = (-1)^n/(2^n n!) d^n/dx^n ((1-x^2)^n) Legendre polynomials are best known for use in Gaussian integration. These approximate a least square fit. integrate f(x) for x= -1 to 1 For n point integration, determine w_i and x_i i=1,n where x_i is the i th root of P_n(x) and w_i = (x_i) integral x=-1 to 1 of f(x) = sum i=1 to n of w_i f(x_i) Change interval from t in [a, b] to x in [-1, 1] x = (b+a+(b-a)t)/2 see example test program gauleg.c and results gauleg_c.out

Chebyshev polynomials

T_n(x) are defined as T_0(x) = 1 T_1(x) = x T_2(x) = 2 x^2 - 1 T_3(x) = 4 x^3 - 3 x T_4(x) = 8 x^4 - 8 x^2 + 1 the general recursion formula for T_n(x) is T_n(x) = 2 x T_n-1(x) - T_n-2(x) The weight function, w(x) = 1/sqrt(1-x^2) The interval is a = -1 to b = 1 Explicit expression T_n(x) = n/2 sum m=0 to n/2 (-1)^m (n-m-1)!/(m! (n-2m)!) (2x)^(n-2m) T_n(x) = cos(n arccos(x)) Chebyshev polynomials are best known for approximating a function while minimizing the maximum error, covered in a latter lecture. Fit f(x) with approximate function F(x) c_i = 2/pi integral x=-1 to 1 f(x) T_i(x) /sqrt(1-x^2) dx (difficult near x=-1 and x=1) F(x) = c_0 /2 + sum i=1 to n c_i T_i(x) Change interval from t in [a, b] to x in [-1, 1] x = (b+a+(b-a)t)/2 see example test program chebyshev.c and results chebyshev_c.out

Laguerre polynomials

L_n(x) are defined as L_0(x) = 1 L_1(x) = -x + 1 L_2(x) = x^2 - 4 x + 2 L_3(x) = -x^3 +9 x^2 - 18 x + 6 the general recursion formula for L_n(x) is L_n(x) = (2n-x-1) L_n-1(x) - (n-1)^2 L_n-2(x) The weight function, w(x) = e^(-x) The interval is a = 0 to b = infinity Explicit expression L_n(x) = sum m=0 to n (-1)^m n!/(m! m! (n-m)!) x^m L_n(x) = 1/(n! e^-x) d^n/dx^n (x^n e^(-x)) Laguerre polynomials are best known for use in integrating functions where the upper limit is infinity.

Lagrange Polynomials

These will be used later in Galerkin Finite Element Method for solving partial differential equations. Lagrange polynomials L_n(x) are defined as L_1,0(x) = 1 - x L_1,1(x) = x L_2,0(x) = 1 -3 x +2 x^2 L_2,1(x) = 4 x -4 x^2 L_2,2(x) = - x +2 x^2 L_3,0(x) = 1 -5.5 x +9 x^2 -4.5 x^3 L_3,1(x) = 9 x -22.5 x^2 +13.5 x^3 L_3,2(x) = -4.5 x +18 x^2 -13.5 x^3 L_3,3(x) = x -4.5 x^2 +4.5 x^3 L_4,1(x) = 16 x -69.33 x^2 + 96 x^3 - 42.33 x^4 L_5,2(x) = -25 x +222.92 x^2 -614.58 x^3 +677.08 x^4 -260.42 x^5 For 0 ≤ x ≤ 1 and equal spaced points in between. The general recursion formula for L_n+1,j(x) is unavailable Explicit expression L_n,k(x) = product i=0 to n, i/=k of (x-x_i)/(x_k-x_i) Lagrange polynomials are best known for solving differential equations with equal and unequal grid spacing. An nth order Lagrange polynomial will exactly fit n data points. Fit f(x) with approximate function F(x) F(x) = sum k=0 to n of f(x_k) L_n,k(x)

Change of coordinates for integration

Change interval from t in [a, b] to x in [0, 1] x = (t-a)/(b-a) Change interval from t in [a, b] to x in [-1, 1] x = 2(t-a)/(b-a)-1

Example plots

see example test program lagrange.c and results lagrange_c.out A family of Lagrange Polynomials can be constructed to be 1.0 at x_i and zero at every x_j where i ≠ j Using the notation above, np=11 and each color is a different k. Each of the 11 points has one color at 1.0 and all other colors at 0.0 . Lagrange Phi shape functions First derivative of Lagrange Phi shape functions

polynomial functions in python3

test_poly.py3 source code of functions test_poly_py3.out output test_poly2.py3 source code of functions test_poly2_py3.out output

python3 divide integers produces floating point

test.py3 source code of functions test_py3.out output Additional utility functions: cs455_l37.html
    <- previous    index    next ->

Other links

Go to top