Math'375 Interpolation

1
Math.375 Interpolation

• Vageli Coutsias, Fall 05

2
(No Transcript)
3
function plot_temp() L,K,T av_temp plot(L,T(,
1),'g',L,T(,2),'r',L,T(,3),'b',L,T(,4),'k')
title(' The temperature as a function of latitude
') xlabel(' Latitude (in degrees) ') ylabel('
Temperature ') legend('K.67','K1.5','K2.0','K
3.0') function L,K,T av_temp QS Table
3.1 L 65-10-55 K .67, 1.5,2.0,3.0 T
is 13x4 T -3.10, 3.52, 6.05, 9.30...
-3.22, 3.62, 6.02, 9.30... -3.30, 3.65,
5.92, 9.17... -3.32, 3.52, 5.70, 8.82...
-3.17, 3.47, 5.30, 8.10... -3.07, 3.25,
5.02, 7.52... -3.02, 3.15, 4.95, 7.30...
-3.02, 3.15, 4.97, 7.35... -3.12, 3.20,
5.07, 7.62... -3.20, 3.27, 5.35, 8.22...
-3.35, 3.52, 5.62, 8.80... -3.37, 3.70,
5.95, 9.25... -3.25, 3.70, 6.10, 9.50
4
Topics

• Polynomial Interpolation in 1d Lagrange
  Chebyshev
• POLYFIT POLYVAL
• 2D Matrices and Color Images
• Trigonometric Interpolation
• Sound

5
function runge1() close all a -1 b 1 m
17 Equally Spaced, m17 x linspace(a,b,m)
y rung(x, 16) x0 linspace(-1,1,100)
y0 rung(x0,16) P,S polyfit(x,y,17) pVal
polyval(P,x0) plot(x,y,'ro',x0,pVal,'b',x0,y0
,'r--') axis(-1,1,-2,2) function y
rung(x,a) y 1./(1ax.2)
6
P,S polyfit(x,y,17) pVal polyval(P,x0)
7
least squares fit degre 8 fit based on 17
pts. P,S polyfit(x,y,8) pVal polyval(P,x0)
8
for i 1m1 Chebyshev
Nodes x(i,1) .5(ba) .5(b-a)cos(pi(i-1)/
m) end y 1./(116x.2)
 
P,S polyfit(x,y,17) pVal polyval(P,x0)
9
Machine assignment 2.7 (a) Write a subroutine
that produces the value of the interp.
polynomial
at any real t, where is a
given integer, are n1 distinct
nodes, and f is any function available in the
form of a function subroutine. Use Newtons
interpolation formula and generate the divided
differences (from the bottom up) in place in
single array of dim n1 that originally contains
the function values. (b) Run your routine on the
function using
and n228 (Runge) Plot the
polynomials against the exact function.
10
function c InterpNRecur(x,y) n
length(x) c zeros(n,1) c(1) y(1)
if n gt 1 c(2n) InterpNRecur(x(2n),
((y(2n)-y(1))./(x(2n)-x(1))) end
11
function c InterpN(x,y) n
length(x) for k 1n-1 y(k1n)
(y(k1n)-y(k))/(x(k1n)-x(k)) end
c y
12
function pVal HornerN(c,x,z) n
length(c) pVal c(n)ones(size(z)) for
kn-1-11 pVal (z-x(k)).pVal c(k)
end
13
Ma505 Gautschi M.2.7 for m14 n2m
x linspace(-5,5,n1) y 1./(1x.2)
x0 linspace(-5,5,1000) y0 1./(1x0.2)
a InterpNRecur(x,y) pVal
HornerN(a,x,x0) subplot(2,2,m)
plot(x0,y0,x0,pVal,'--',x,y,'') axis(-5 5
-1 1) end
14
(No Transcript)
15
LAGRANGE polynomials
16
Lagrangian Polynomial
Interpolating Polynomial is a UNIQUE polynomial
of degree n that agrees with f(x) at the n1 nodes
17
(No Transcript)
18
Problem 2.24 Prepare plots of the Lebesgue
function for interpolation,
for with the
nodes given by Compute on
a grid obtained by dividing each interval
into 20 equal
subintervals. Plot
in case (a) and in case (b).

19
The Lebesgue function and constant
for a grid are defined
as with the Lagrange
polynomial Below we give a matlab script that
computes the Lebesgue functions (and constants)
for the given grids.
20
Math 505 Gautchi 2.24 Lebesgue
functions fprintf(' n equispaced
Chebyshev \n') for ii02 n52ii
xa(1n1,1)0 xb(1n1,1)0 for i 1n1
xa(i,1) -12(i-1)/n xb(n2-i,1)
cos(pi(2(i-1)1)/(2n2)) end for i1n
for j120 ga(20(i-1)j,1)((21-j)
xa(i)(j-1)xa(i1))/20
gb(20(i-1)j,1)((21-j)xb(i)(j-1)xb(i1))/20
end end ga(20n1,1)xa(n1)
gb(20n1,1)xb(n1) for k 120n1
compute Lebesgue functions
na(,k)ga(k)-xa() nb(,k)gb(k)-xb()
numerators end for i 1n1
denominators for j1n1
da(i,j) xa(i)-xa(j) db(i,j)
xb(i)-xb(j) end end
21
for k 120n1 for i 1n1
la(i,k) 1 lb(i,k) 1 for j1n1
if ji
la(i,k)la(i,k)na(j,k)/da(i,j)
lb(i,k)lb(i,k)nb(j,k)/db(i,j)
end end end lama(k)1
lamb(k)1 for i1n1
lama(k)lama(k)abs(la(i,k))
lamb(k)lamb(k)abs(lb(i,k)) end end
figure subplot(2,1,1) plot(ga,log(lama))
subplot(2,1,2) plot(gb,lamb) lamax
lama(1) lambx lamb(1) compute Lebesgue
constants for k 220n1 if lama(k) gt
lamax lamax lama(k) end if lamb(k)
gt lambx lambx lamb(k) end end
fprintf('9i 22.14f 22.14f \n',n,lamax,lambx)
clear end
22
a
b
n5
n10
Lebesgue constants n equispaced
Chebyshev 5 4.10496
2.68356 10 30.8907 3.06859 20
10987.5 3.47908 The plots
give the Lebesgue fncs (log of same for grid (a)).
a
b
n20
23
Fourier series
24
(No Transcript)
25
(No Transcript)
26
Complex Form
27
Band Limited Functions
Complex Form
28
Real to Complex conversion
29
Equispaced Sampling
30
Discrete Orthogonality
31
Coefficients are computed exactly from discrete
sampling
32
Real form even or odd sample sizes
33
f
t
34
The Fast Fourier transform

• fft compute Fourier coefficients
• from sampled function values
• ifft reconstruct sampled values
• from (complex) coefficients
• interpft construct trigonometric
• interpolant

35
FFT Discrete Fourier transform. FFT(X) is
the discrete Fourier transform (DFT) of vector X.
For length N input vector x, the DFT is
a length N vector X, with elements
N X(k) sum
x(n)exp(-j2pi(k-1)(n-1)/N), 1 lt k lt N.
n1 The inverse DFT
(computed by IFFT) is given by
N x(n) (1/N) sum X(k)exp(
j2pi(k-1)(n-1)/N), 1 lt n lt N.
k1 See also IFFT, FFT2, IFFT2,
FFTSHIFT.
36
INTERPFT 1-D interpolation using FFT method.
Y INTERPFT(X,N) returns a vector Y of
length N obtained by interpolation in the
Fourier transform of X. Assume x(t) is a
periodic function of t with period p,
sampled at equally spaced points, X(i)
x(T(i)) where T(i) (i-1)p/M, i 1M, M
length(X). Then y(t) is another
periodic function with the same period and
Y(j) y(T(j)) where T(j) (j-1)p/N, j
1N, N length(Y). If N is an integer
multiple of M, then Y(1N/MN) X.
37
function F CSInterp(f) n length(f) m
n/2 tau (pi/m)(0n-1)' P for j 0m,
P P cos(jtau) end for j 1m-1, P
P sin(jtau) end y P\f F
struct('a',y(1m1),'b',y(m2n))
38
THE FOURIER COEFFICIENTS
F.ay(1),,y(m1)
F
F.by(m2),,y(2m)
A CELL array its elements are STRUCTURES
(here arrays) F struct('a',y(1m1),'b',y(m2
n))
39
function Fvals CSeval(F,T,tvals)
Fvals zeros(length(tvals),1) tau
(2pi/T)tvals for j 0length(F.a)-1
Fvals Fvals F.a(j1)cos(jtau)
end for j 1length(F.b) Fvals
Fvals F.b(j )sin(jtau) end
40
Script File Pallas A linspace(0,360,13)' D
408 89 -66 10 338 807 1238 1511 1583 1462 1183
804 408' Avals linspace(0,360,200)' F
CSInterp(D(112)) Fvals CSeval(F,360,Avals) pl
ot(Avals,Fvals,A,D,'o') axis(-10 370 -200
1700) set(gca,'xTick',linspace(0,360,13)) xlabel(
'Ascension (Degrees)') ylabel('Declination
(minutes)')
41
The Gibbs Phenomenon .
42
Script File square Trigonometric interpolant
of a square wave. M300 A linspace(0,360,2M1)'
D(2M)0D(M22M)1 D(1).5D(M1).5D(2M1
).5 DD' Avals linspace(0,360,1000)' F
CSInterp(D(12M)) Fvals CSeval(F,360,Avals)
plot(Avals,Fvals)
43
. can be seen at all discontinuities
44
Script File sawtooth_wave. clear M512 A
linspace(0,360,2M1)' D(2M)(2M)/MD(M22M)
-1(2M)/M D(1)0D(M1)0D(2M1)0 DD' Aval
s linspace(0,360,1024)' F
CSInterp(D(12M)) Fvals CSeval(F,360,Avals) p
lot(Avals,Fvals)
45
The overshoot is present at any resolution
46
although the interpolant is correct at the nodes
as dx-gt0
47
Compare with equispaced polynomial interpolation
48
Input-output structures

• xinput(string)
• titleinput(Title for plot,string)
• x,yginput(n)
• disp(a 3x3 magic square)
• disp(magic(3))
• fprintf(5.2f \n5.2f\n,pi2, exp(1))

49
Writing to and reading from a file

• fid fopen('output','w') alinspace(1,100,10
  0)
• fprintf(fid,'g degC g degF\n',a,321.8a)
• fidfopen('output','r')
• xfscanf(fid,'g degC g degF')
• xreshape(x,(length(x)/2),2)

50
---------------------------------------- I. A1
2 34 5 67 8 9 yx' transpose xzeros(1,n) l
ength(x) xlinspace(a,b,n) ax(5-12) SineTab
le (Help SineTable) disp(sprintf('2.0f 3.0f
',k, a(k))) -----------------------------------
-----
51
--------------------------------------------------
-------------- II. VECTOR OPERATIONS (1)
vector scale, add, subtract (scalarvector)
210 20 30 20 40 60 (2) POINTWISE vector
multiply, divide, exponentiate (pointwise vector
vector) 2 3 4.10 20 30 20 60
120 (for pointwise operations, dimensions must
be identical!) (3) Setting ranges for axes,
Superimposing plots plot(x,y,x,exp(x),'o')
axis(-40 0 0 0.1) using "hold on/off"
extended parameter lists (4) Subplots
subplot(m,n,k) (5) Vector linear combos
matrixvector (6) Tensor Products (x is array, A
is array of arrays of same shape as x)
Af1(x) f2(x) ... fn(x) A(k,) k-th row
A(,i) i-th column (7) loops loop until
condition or for fixed no. of times while
any(abs(term) gt epsabs(y)) for term 1n
end (8) size(A) gives array dimensions (array,
can be used to index other arrays) length(x)
gives vector dimension ---------------------------
-------------------------------------
52
--------------------------------------------------
-------------- II' (Miscellaneous
topics) rem(x,n) remainder if x 1 (lt lt
gt gt ) ( and or not xor exclusive or)
else endif s3 s1 s2 concatenation of
strings, treated like arrays xmax, imax
max(x) (max val, index) rat max(x)/x(1) x(1)
input('Enter initial positive integer') density
sum(xltx(1))/x(1)
53
Summary

• Interpolation
• vanderMonde, Newton, Least Squares
• 2d interpolation linear and cubic
• Color figures and uint8
• Smoothing of images

54
References

• find m-files at
• www.math.unm.edu/vageli
• Higham Higham,
• Matlab Handbook
• C.F.van Loan,
• Intro to Sci.Comp. (2-3)