CSc 2262 Numerical Methods Fall 2004[1]

Notes

History

Calculating devices over time – ancient, middle, modern

Babbage, analytical (including slide rule), digital

WWII – scientific computing

Supercomputers and grand challenge problems – weather, human genome,

oil exploration, …

Mathematics and Electrical Engineering roots – circuits, hardware;

algebra, trigonometry, calculus

Number Representation

Decimal numbers, binary digits, octal and hexadecimal digits – conversion

Floating point numbers – sign, exponent, fraction (IEEE standard)

characteristic – 8 bit, excess 127

sign bit

mantissa or significant (fraction) – n=23 bit single precision,

n=52 bit double precision

accuracy – single 2^-23 or 10^-7; double 2^-52 or 10^-16

Errors

Can only store n bits of fraction – what to do with 1/3 or p?

F(x) = version of x stored as floating point = x(1+e)

Chopping – drop bits if more than n è eÎ[-2^-n,0]

Rounding – if digit n+1=0 then chop else chop and add 1 if that digit=1

so eÎ[-2^-n,2^-n]

è worst chopping error twice rounding error, sign of f(x) same as sign of

x so no cancellation possible

*Consequences

Error = true value–approximate value=x_t-x_a

Relative error = error/true value =(x_t-x_a)/x_t=1x_t/x_t

Example: p=3.14159265…= x_t while 22/7=3.1428571…=x_a

Error = 0.00126 and Relative error = 0.000402

Error sources

Modeling errors N_t=number of people on earth in year t = N₀e^kt

Out of scope but model does not hold based upon real data

Mistakes (blunders) – programmers

Physical measurement – c = speed of light = 2.997925+e *10¹⁰ cm/sec

Machine representation and arithmetic errors – rounding, chopping, underflow,

Overflow

Mathematical approximation and truncation (digitization) - p

Loss of significance

*Noise

Function, fuzzy curve

Program segment example:

n=30; x=1.0/9.9;

for (i=1; i<=n;i++) {x=10.0*x-1.0; print(i,x);}

Propagation (interval arithmetic)

x_a+y_a |(x_t+y_t)|-(x_a+y_a)|£e_x+e_y

x_a/y_a (x_a-e_x)/(y_a+e_y)£x_t/y_t)|£ x_a+e_x)/(y_a-e_y)

relative error of (x*y) = relative error(x)+relative error(y) – relative

error(x)*relative error(y)

relative error ~ f’(x_t)(x_t-x_a)/f(x_t)= f’(x_t)/(x_t) x_t relative error(x_a)

Theorems

intermediate value theorem

f(x) continuous xÎ[a,b]

M=max_xf(x) and m=min_xf(x)

"vÎ[m,M] $cÎ[a,b] ' f(c)=v

mean value theorem

f(x) continuous and f’(x) continuous and differentiable xÎ[a,b]

$cÎ[a,b] ' [f(b-f(a)]=(b-a)*f’(c)

integral mean value theorem

f(x) continuous and w(x) nonnegative and integrable xÎ[a,b]

$cÎ[a,b] ' ò_a^bf(x)w(x)dx=f(c)ò_a^b w(x)dx

Taylor Series

f(x) continuous with n+1 derivatives xÎ[a,b]

P_n(x)=S_j=0ⁿf^j(x₀)(x-x₀)^j/j!

R_n(x) = remainder = E_n(x)= error =f(x)-P_n(x)=(x-x₀)ⁿ⁺¹fⁿ⁺¹(c)/(n+1)!

j!=j(j-1)(j-2)…(2)(1)=factorial and 0!=1 f^j(x₀)=jth derivative of f at x₀

Note: if x₀=0 this is called a Maclaurin series

Example: f(x)=e^x and x₀=0 xÎ (-¥,¥)

f^j(x₀)=e^x₀=e⁰=1 "j=0,1,…,n+1

P_n(x)=1+x+x²/2+…+xⁿ/n! and R_n(x)=xⁿ⁺¹/(n+1)!e^c for cÎ[0,x]

Evaluation of Polynomials (Horner method)

P_n(x)=S_j=0ⁿa_jx^j for real coefficients a_j Evaluate P_n(x) at x=z

Let b_n= a_n then b_n-1=a_n-1+zb_n, b_n-2=a_n-2+zb_n-1,…,b_k=a_k+zb_k+1,…, b₀=a₀+zb₁

then P_n(z)= b₀

Synthetic Division

Let Q_n-1(x)=S_j=1ⁿb_jx^j-1= quotient of P_n(x)/(x-z) with b₀ as remainder è

P_n(x)=b₀+(x-z) Q_n-1(x) (x)

Roots e=error tolerance

Find x such that f(x)=0

Example: Suppose that one puts in P_in dollars into a savings account at 100*r % interest compounded in each of N_in time periods, and after that withdraws P_out dollars from the account for N_out time periods until the account has no money left. After N_in time periods the account has S_in= P_in (1+r)[(1+r)^N_in-1]/r dollars and the amount after N_out withdrawals of P_out dollars = f® = (1+r)^N_out^-1S_in-P_out[(1+r)^N_out –1]/r = 0. We can solve for the r needed for this transaction.

Bracketing Methods f(x) continuous for xÎ[a,b]

Incremental Search multiple roots

set a_j=(a+(b-a)j/n

if f(a_j)*f(a_j+1)<0 then output [a_j,a_j+1] as small interval with root

Note: could then use other bracketing methods to get more exact estimates

Bisection f(a)*f(b)<0

set c=(ab)/2; if b-c£e then c=root and stop

elseif f(b)*f(c)£ 0 then set a=c elseif set b=c; restart

Note: b_n- a_n=2^n-1(b-a) and root is in [a_n,b_n] so must converge as n increases

è n³log(b-a)/e)/log(2)

False Position f(a)*f(b)<0

x₀,x₁Î[a,b]; x=x₁-f(x₁)*(x₁-x₀)/[f(x₁)-f(x₀)];

if f(x)*f(b)<0 then x₀=x else x₁=x;

if |x₁-x|<Î else restart

Open Methods (no bracketing)

Fixed Point Iteration f(x) continuous

find g(x) such that root of f(x) solves g(x)=x

one possibility for g is g(x)=bf(x)+x for some nonzero parameter b

find original estimate x₀; set x=g(x0);

if |x-x₀|<Î stop; else restart with x₀ = x

Secant

analogous to false position but without the brackets

find initial estimates x₀,x₁;

x=x₁-f(x₁)*(x₁-x₀)/[f(x₁)-f(x₀)];

if f(x)*f(x₁)<0 then x₀=x else x₁=x;

if |x₁-x|<Î else restart

Newton-Raphson

based on Taylor series f(x) continuous and differentiable

find x₀; x=x₀-f (x₀)/f’(x₀);

if |x-x₀|<Î stop; else restart with x₀ = x

Matlab

Can use fzero to find roots

Can use polyeval to evaluate a polynomial

Can use poly to construct a polynomial given its roots

Can use conv to multiply a polynomial by another polynomial

Can use deconv to divide one polynomial into another polynomial

Interpolating (or Approximating) Polynomials

Computers have numbers via binary representations and can do simple arithmetic (+,-,*,/)

Need to estimate complex functions (e.g., e^x or sin(x)) via these simple arithmetic

functions è estimating polynomials (using multiplication and addition)

Interpolation – find P_n(x) as an interpolating polynomial that has the same value as the function through some known points, i.e., f(x_j)=P_n(x_j). Note that the known points are x₀,x₁,…, x_n and y₀,y₁,…, y_n where y_j=f(x_j).

Taylor Series x₀Î[a,b]

P_n(x)=S_j=0ⁿf^j(x₀)(x-x₀)^j/j!

E_n(x)= error =f(x)-P_n(x)=(x-x₀)ⁿ⁺¹fⁿ⁺¹(c)/(n+1)!

Linear Interpolation x₀, x₁Î[a,b]

P₁(x)=[(x₁-x)y₀+(x-x₀)y₁]/(x₁-x₀)

Note: slope = (y₁-y₀)/(x₁-x₀) and y-intercept = (x₁y₀-x₀y₁)/(x₁-x₀)

Note that P₁(x_i)=yi for i=0,1

Example: f(x)=Öx x₀=1 and x₁=4 so y₀=1 and y₁=2

Slope = (2-1)/(4-1)=1/3 and Y-intercept=(4*1-1*2)/(4-1)=2/3

P₁(x)= (1/3)x+(2/3)

for x=2, f(x)=1.4142135 and P₁(x)=4/3=1.3333333

Example: f(x)=e^x x₀=0.82 and x₁=0;83 so y₀=2,2705 and y₁=2.293319

for x=2, f(x)=2.2841638 and P₁(x)=2.2841914

Lagrange Polynomials x₀, x₁,…,x_nÎ[a,b]

P_n(x)=S_j=0ⁿ y_jL_j(x) where L_j(x)=P_i=0,i¹jⁿ[(x-x_i)/(x_j-x_i)]

Note that L_j(x_i)=δ_ij=1 if i=j, 0 if i¹j (δ=Kronecker delta)è P_n(x_j)=y_j

Note that polynomial can be shown to be unique

E_n(x)= error of P_n(x) at or x = P_i=0ⁿ(x-x_i) fⁿ⁺¹(c)/(n+1)! for cÎ[a,b]

Example: n=2; f(x)=e^x x₀=0.82, x₁=0;83, and x₂=0.84, so y₀=2,2705,

y₁=2.293319, and y₂=2.316367

for x=0.826, f(x)=2.2841638 and P₂(x)=2.2841639

Newton Polynomials x₀, x₁,…,x_nÎ[a,b]

P_n(x)=S_j=0ⁿ D_j*P_i=0^j-1(x-x_i)

Divided Differences

D₀=f[x₀]=f(x₀);

D₁=f[x₀,x₁]=[f(x₁)-f(x₀)]/(x₁-x₀)={f[x₁]-f[x₀]}/(x₁-x₀)

~=f’(c) and c~{x₀+ x₁}/2

D₂=f[x₀,x₁,x₂]={f[x₁,x₂]-f[x₀,x₁]}/x₂-x₀} …

D_n=f[x₀, x₁,..., x_n]={f[x₁,…, x_n]-f[x₀,…, x_n-1]}/{x_n-x₀} ~= f⁽ⁿ⁾(c)/n!

Note that we can permute the order of x₀ variables

E_n(x)=error = P_i=0ⁿ (x-x_i)f⁽ⁿ⁺¹⁾(c)/(n+1)! for cÎ[min x_j,max x_j]

Note that we can calculate recursively:

d_k,0=f(x_k), d_k,j=[d_k,j-1-d_k-1,j-1]/(x_k-x_k-j) è D_k=d_k,k

Approximation – find P_n(x) as an approximating polynomial minimizing the error, at

least at some known points x₀,x₁,…, x_n and y₀,y₁,…, y_n where y_j=f(x_j).

Note: error is a function of P_n(x_j)-y_j at x_j

Chebyshev Polynomials x₀, x₁,…,x_nÎ[-1,1]

E_n(x)=P_i=0ⁿ (x-x_i)f⁽ⁿ⁺¹⁾(c)/(n+1)! Q(x)= P_i=0ⁿ (x-x_i)

Minimize Q(x) by picking x_k=cos[(2k+1)P/2n]

P_n(x)=S_j=0ⁿ a_jC_j(x) where

C_j(x)=cos(j cos^-1(x))

C₀(x)=1; C₁(x)=1; C_j(x)=2xC_j-1(x)-C_j-2(x)

a₀=[1/(n+1)]S_k=0ⁿy_kC₀(x_k)= [1/(n+1)]S_k=0ⁿy_k

a_j=[2/(n+1)]S_k=0ⁿy_kC_j(x_k)=[2/(n+1)]S_k=0ⁿy_kcos[jP(2k+1)/(2n+2)]

Transforming the interval

if xÎ[a,b] and a¹-1 and/or b¹1, consider tÎ[-1,1] by letting

t=2[(x-a)/(b-a)]-1 which means x=[(b-a)t/2]+[(a+b)/2]

Select t_k=cos[(2n+1-2k)P/(2n+2)] so x_k=[(b-a)t_k/2]+[(a+b)/2]

Legendre Polynomials x₀, x₁,…,x_nÎ[a,b]

Minimize E = error = Ö[1/(b-a)]ò_a^b[f(x)-Pn(x)]²dx

P₀(x) = 1; P_j(x)=i/(j!2^j)d^j[(x²-1)^j]/dx^j

P_n(x)=S_j=0ⁿ b_jLG_j(x)

LG0(x)=1; LG_j+1(x)=[(2j+1)/(j+1)]xLG_j(x)- [j/(j+1)]LG_j-1(x)

b_j=ò_a^b f(x)LG_j(x)/ ò_a^bLG_j(x)²dx

Splines (thin, flexible strings) x₀áx₁á…áx_nÎ[a,b]

Knot – x values where two splines meet

Let h_i=x_i+1-x_i

Linear

S_i(x)=a_i+b_i(x-x_i)

S_i(x_i)=y_i (continuity) è a_i= y_i è S_i(x)=y_i+b_i(x-x_i)

b_i=(y_i+1-y_i)/(x_i+1-x_i) =(y_i+1-y_i)/h_i = slope

è S_i(x)= y_i+[(y_i+1-y_i)/(x_i+1-x_i)](x-x_i)=y_i+[(y_i+1-y_i)/h_i](x-x_i)

Note that S_i(x_i+1)= S_i(x_i+1)=y_{i+1 (}equality at knots)

Quadratic

S_i(x)=a_i+b_i(x-x_i)+c_i(x-x_i)²

S_i(x)=y_i è a_i= y_i è S_i(x)=y_i+b_i(x-x_i)+c_i(x-x_i)²

S_i(x_i+1)=y_i+b_i(x_i+1-x_i)+c_i(x_i+1-x_i)²=S_i+1(x_i+1)=y_i+1+b_i+1(x_i+1-x_i+1)+c_i(x_i+1-x_i+1)²=y_i+1

è y_i+b_ih_i+c_ih_i²=y_{i+1 (}equality at knots)

S_i’(x)=b_i+2c_i(x-x_i) è S_i’(x_i+1)=b_i+2c_ih_i=S_i+1(x_i+1)=b_i+1 è b_i+2c_ih_i=b_i+1

(equality at knots – first derivative)

Cubic

S_i(x)=a_i+b_i(x-x_i)+c_i(x-x_i)²+d_i(x-x_i)³ è a_i= y_i (continuity)

è S_i(x)=y_i+b_i(x-x_i)+c_i(x-x_i)²+d(x-x_i)³

y_i+b_ih_i+c_ih_i²+d_ih_i³=y_i+1 and b_i+2c_ih_i+3d_ih_i²=b_i+1 (equality of knots for S and S’)

S’(x)= b_i+2c_i(x-x_i)+3d_i(x-x_i)² and S”(x)=2c_i+6d_i(x-x_i) è c_i+3d_ih_i=c_i+1

(equality of knots for S”)

Natural spline c₀=c_n=0 (additional conditions)

Clamped (f’ and f” specified for x₀ and x_n)

è 2h₀c₀+h₀c₁=3(y₁-y₀)/(x₁-x₀)-3f’(x₀)

2h_n-1c_n-1+2h_n-1c_n=3f’(xn)-3(y_n-y_n-1)/(x_n-x_{n -1})

Not-a-knot (continuity of 3^rd derivatives at 2^nd and next-to-last knots – first

internal knots no longer true knots; gives same result using ordinary cubic

polynomial for first four points)

è h₁c₀-(h₀+h₁)c₁+h₀c₂=0

h_N-1c_N-2-(h_N-2+h_N-1)c_N-1+h_N-2cn=0

Least Squares

Y=f(x) or f(x_k) = y_k+e_k or e_k=f(x_k)-y_k k=1,…,n

E_¥(f)=max_k|f(x_k)-y_k| = maximum error; E₁(f)=(1/n)S_k=1ⁿ|f(x_k)-y_k| = average error;

E₂(f)=Ö(1/n)S_k=1ⁿ|f(x_k)-y_k|² = root-mean-square error

Least squares minimizes E₂(f)

Linear f(x)=ax+b E₂(f)=Ö(1/n)S_k=1ⁿ|ax_k+b-y_k|²

normal equations obtained from dE₂(f) d/a=0 and dE₂(f) d/b=0

è aS_k=1ⁿx_k²+bS_k=1ⁿx_k=S_k=1ⁿx_ky_k

aS_k=1ⁿx_k+bn=S_k=1ⁿy_k

è a=(nS_k=1ⁿx_ky_k-S_k=1ⁿx_kS_k=1ⁿy_k )/[nS_k=1ⁿx_k²-(S_k=1ⁿx_k)²]

b=(1/n)[S_k=1ⁿy_k-aS_k=1ⁿx_k]

= g-ac where g=(1/n)[S_k=1ⁿy_k=average y value and

c=(1/n)[S_k=1ⁿx_k=average x value

Matlab – linregr

x = vector of n x values, y = vector of n y values

linregr(x,y,n)

Nonlinear – can try to linearize

y=ae^bx è ln(y)=ln(a)+bxè z=A+bx with z=ln( y) and A=ln(a)

y=ax^bè log(y)=log(a)+blog(x)

è z=A+bw with z=log(y), A=log(a), and w-log(x)

y=ax/(b+x) è (1/y)=(1/a)+(b/a)(1/x)

è z=A+Bw with z=1/y, A=1/a, B=b/a, and w=1/x

Matlab – polyfit, polyeval

x = vector of n x values, y = vector of n y values,

X=value of x to be tested

a=polyfit(x,y,n)

Y=polyeval(a,X)

---

[1] Grader is Ankur Suri, at asuri1@lsu.edu