<- previous index next ->
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
"y" is computed numerically using Horners method that
needs n multiplications and n additions.
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)
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
Having the coefficients of a polynomial allows creating
the derivative and integral of the 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);
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);
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 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, 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!
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.
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 cos^(-1) 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 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
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 grid spacing.
Fit f(x) with approximate function F(x)
F(x) = sum k=0 to n of f(x_k) L_n,k(x)
Change interval from t in [a, b] to x in [0, 1]
x = (a+(b-a)t)
see example test program lagrange.c
and results lagrange_c.out
<- previous index next ->