Transcript Math.375 III- Interpolation

Math.375
Interpolation
Vageli Coutsias, Fall `05
7/21/2015
2
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','K=1.5','K=2.0','K=3.0')
function [L,K,T] = av_temp
% Q&S 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;...
7/21/2015
-3.37,
3.70, 5.95, 9.25;...
-3.25, 3.70, 6.10, 9.50];
3
Topics
• Polynomial Interpolation in 1d:
Lagrange & Chebyshev
• POLYFIT & POLYVAL
• 2D: Matrices and Color Images
• Trigonometric Interpolation
• Sound
function runge1()
close all
a = -1; b = 1; m = 17; %Equally Spaced, m=17
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 =7/21/2015
1./(1+a*x.^2);
5
[P,S] = polyfit(x,y,17);
pVal = polyval(P,x0);
7/21/2015
6
% least squares fit: degre 8 fit based on 17 pts.
[P,S] = polyfit(x,y,8);
pVal = polyval(P,x0);
7/21/2015
7
for i = 1:m+1
% Chebyshev Nodes
x(i,1) = .5*(b+a) +.5*(b-a)*cos(pi*(i-1)/m) ; end
y =1./(1+16*x.^2);
[P,S] = polyfit(x,y,17);
pVal = polyval(P,x0);
7/21/2015
8
Machine assignment 2.7
(a) Write a subroutine that produces the value of the interp.
polynomial
pn  f ; x0 , x1,, xn ; t  at any real t, where
n0
is a given integer, x i
are n+1 distinct nodes,
and f is any function available in the form of a function
subroutine. Use Newton’s interpolation formula and generate
the divided differences (from the bottom up) “in place” in
single array of dim n+1 that originally contains the function
values.
1
,5  t  5
(b) Run your routine on the function f t  
2
1 t
i
using xi  5  10 , i  0,1,  , n. and n=2:2:8 (Runge)
n
Plot the polynomials against the exact function.
7/21/2015
9
function c = InterpNRecur(x,y)
n = length(x); c = zeros(n,1);
c(1) = y(1);
if n > 1
c(2:n) = InterpNRecur(x(2:n), ((y(2:n)-y(1))./(x(2:n)-x(1)));
end
n 1
pn x   pn1 x   cn  x  xi 
i 0
7/21/2015
pn x0   y0  c0 ; pn xi   yi
10
function c = InterpN(x,y)
n = length(x);
for k = 1:n-1
y(k+1:n) = (y(k+1:n)-y(k))/(x(k+1:n)-x(k));
end
c = y;

x0 ,, xn1  f  x1 ,, xn  f
x0 ,, xn  f 
xn  x0 
x0 , x1,, xn y : cn ; x0 y  c0 : y0
7/21/2015
11
function pVal = HornerN(c,x,z)
n = length(c);
pVal = c(n)*ones(size(z));
for k=n-1:-1:1
pVal = (z-x(k)).*pVal + c(k);
end
px  c0  x  x0 c1  x  x1 c2  
7/21/2015
12
% Ma505; Gautschi M.2.7
for m=1:4
n=2*m;
x = linspace(-5,5,n+1);
y = 1./(1+x.^2);
x0 = linspace(-5,5,1000);
y0 = 1./(1+x0.^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
7/21/2015
13
7/21/2015
14
n
x  xj
j 0
xk  x j
k  x   
jk
, k  0,  , n
k  Pn , k x j    jk  
1, j  k
0, j  k
LAGRANGE polynomials
7/21/2015
15
Lagrangian Polynomial
Nodes: x 
n
i i 0
Values :  f xi 
n
i 1
Interpolating
Polynom ial:  n f xi   f xi , i  0,, n
Interpolating Polynomial is a UNIQUE
polynomial of degree n that agrees with f(x)
7/21/2015
16
at the n+1 nodes
n
 n f x    f xk k x 
k 0
Error

 n


En f  x   f  x    n f  x  
x

x
i

n  1! i 0
 n 1

max f
x
n 1
xI
max En f x  
h
xI
4n  1
f
7/21/2015
 n 1
17
Problem 2.24
Prepare plots of the Lebesgue function for interpolation,
n x ,1  x  1,
for
n  5,10,20
with the nodes
xi
given by:
2i
2i  1
a xi  1  , b xi  cos
 , i  0,1,, n
n
2n  2
Compute
interval
Plot
n x 
xi1, xi , i  1,2,, n,
log10 n x 
7/21/2015
on a grid obtained by dividing each
into 20 equal subintervals.
in case (a) and
n x 
in case (b).
18
The Lebesgue function
are defined as:
n
n x  and constant n x
n x    li x 
with
li x 
i 0
n x  n
the Lagrange polynomial:

for a grid
n
x  xj
j 0
xi  x j
li x   
j i
Below we give a matlab script that computes the Lebesgue
functions (and constants) for the given grids.
7/21/2015
19
% Math 505; Gautchi 2.24: Lebesgue functions
fprintf('
n
equispaced
Chebyshev \n')
for ii=0:2
n=5*2^ii; xa(1:n+1,1)=0; xb(1:n+1,1)=0;
for i= 1:n+1
xa(i,1) = -1+2*(i-1)/n;
xb(n+2-i,1) = cos(pi*(2*(i-1)+1)/(2*n+2));
end
for i=1:n
for j=1:20
ga(20*(i-1)+j,1)=((21-j)*xa(i)+(j-1)*xa(i+1))/20;
gb(20*(i-1)+j,1)=((21-j)*xb(i)+(j-1)*xb(i+1))/20;
end
end
ga(20*n+1,1)=xa(n+1);
gb(20*n+1,1)=xb(n+1);
for k = 1:20*n+1
%compute Lebesgue functions
na(:,k)=ga(k)-xa(:); nb(:,k)=gb(k)-xb(:);
%numerators
end
for i = 1:n+1
% denominators
for j=1:n+1
da(i,j) = xa(i)-xa(j); db(i,j) = xb(i)-xb(j);
end
end7/21/2015
20
for k = 1:20*n+1
for i = 1:n+1
la(i,k) = 1; lb(i,k) = 1;
for j=1:n+1
if j~=i
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 i=1:n+1
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 = 2:20*n+1
if lama(k) > lamax; lamax= lama(k); end
if lamb(k) > lambx; lambx= lamb(k); end
end
fprintf('%9i %22.14f %22.14f \n',n,lamax,lambx); clear
end 7/21/2015
21
a
b
n=5
n=10
Lebesgue constants
n
equispaced
5
4.10496
10
30.8907
20 10987.5
Chebyshev
2.68356
3.06859
3.47908
a
b
The plots give the Lebesgue fncs
(log of same for grid (a)).
7/21/2015
22
n=20
Fourier
series
a0 
f t     ak coskt   bk sin kt 
2 k 1
7/21/2015
2
a
 k 1
 cos kt 
    f t 
dt
 sin kt  
 bk   t 0
23
a0 m1 
 kn 
 kn  
n
f n     ak cos
  bk sin
    1 am
2 k 1 
 m 
 m 
n  0,,2m  1
2
1 2 m1
 kn 
ak   f t  coskt dt   f n cos

 t 0
m n 0
 m 
1
2
1
bk   f t sin kt dt 
 t 0
m
1
7/21/2015
 kn 
f n sin


 m 
n 0
2 m 1
24
C
1

1




 n 

cos


1
 m 




1 cos 2m  1  

m



1

 1n

 12 m1
0



 n 
sin 


 m 


 2m  1 
sin 
 
m






 nm  1  
sin 

m






n
2
m

1



sin 

m


0
F  a, b  a0 , a1 ,, am , b1 ,, bm1 
^
T
f   f0 , f1,, f 2m1 (, f 2m  f0 )
T
^
7/21/2015
^
1
CF  f  F C f
25
Complex Form
f x  
n 2 
c
e
k
ikx
k   n 2 
7/21/2015
26
Band Limited Functions
a0 M
f t     ak coskt   bk sin kt   aM 1 cosM  1x 
2 k 1
Complex Form
f x  
7/21/2015
M ( 1)
ce



ikx
k
k   M ( 1)
27
Real to Complex conversion
a k  ck  c k
bk  i ck  c k 
k  0,  , M
aM 1 

 cM 1  c M 1 

2 

7/21/2015
28
Equispaced Sampling
M 
f x  
ce

 
k  M 
f x j  
j 0
7/21/2015
j
k
M 
ce

 
k  M 
 f x e
n
ikx
imjh
ikjh
, x j  kh
n
M 
k

ce

 
j 0 k   M 
k
ikjh imjh
e
29
Discrete Orthogonality
ih k  m  n 1
0, k  m
1

e
ijh k  m 
e

 n  1 km  

ih k  m 
n  1, k  m
1 e
j 0
n
 f x e
n
j
j 0
 f x e
n
j 0
7/21/2015
imjh
j
imjh

M 
c e

 
k  M 

n
k
ikjh imjh
j 0
M 
c n  1

 
k  M 
k
e
mk
 cm n  1
30
Coefficients are computed exactly
from discrete sampling
 f x e
n
j
j 0
imjh

M 
c n  1

 
k  M 
k
mk
 cm n  1
1
imjh
cm 
f x j e

n  1 j 0
n
7/21/2015
31
Real form; even or odd sample sizes

 kj 
 kj  
f j  a0    ak cos
  bk sin
 
 m 
 m 
k 1 
j  0,,2m
m

 kj 
 kj  
j
f j  a0    ak cos
  bk sin
    1 am
 m 
 m 
k 1 
j  0,,2m  1
m 1
7/21/2015
32
n  0,,2m  1
f
fn
n
tn 
m
7/21/2015
t
33
The Fast Fourier transform
• fft: compute Fourier coefficients
from sampled function values
• ifft: reconstruct sampled values
from (complex) coefficients
• interpft: construct trigonometric
interpolant
7/21/2015
34
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(-j*2*pi*(k-1)*(n-1)/N), 1 <= k <= N.
n=1
The inverse DFT (computed by IFFT) is given by
N
x(n) = (1/N) sum X(k)*exp( j*2*pi*(k-1)*(n-1)/N), 1 <= n <= N.
k=1
7/21/2015
See also IFFT, FFT2, IFFT2, FFTSHIFT.
35
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 = 1:M, 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 = 1:N, N = length(Y).
If
N
is
an
integer
multiple
of
M,
then
7/21/2015
36
Y(1:N/M:N) = X.
function F =CSInterp(f)
n = length(f);
m = n/2;
tau = (pi/m)*(0:n-1)';
P = [];
for j = 0:m, P = [P cos(j*tau)]; end
for j = 1:m-1, P = [P sin(j*tau)]; end
y = P\f;
F = struct('a',y(1:m+1),'b',y(m+2:n));

 k 
 kn  
n
f n  a0    ak cos   bk sin
    1 am
m
 m 
k 1 
n  0,,2m  1
m 1
7/21/2015
37
THE FOURIER COEFFICIENTS
F
F.a=[y(1),…,y(m+1)]
F.b=[y(m+2),…,y(2m)]
A CELL array: it’s elements are STRUCTURES
(here arrays)
F = struct('a',y(1:m+1),'b',y(m+2:n));
7/21/2015
38
function Fvals = CSeval(F,T,tvals)
Fvals = zeros(length(tvals),1);
tau = (2*pi/T)*tvals;
for j = 0:length(F.a)-1
Fvals = Fvals + F.a(j+1)*cos(j*tau);
end
for j = 1:length(F.b)
Fvals = Fvals + F.b(j )*sin(j*tau);
end
m 1

 2kt 
 2kt  
f t   a0    ak cos
  bk sin 
 
 T 
 T 
k 1 
 2m t 
  am cos

 T 
7/21/2015
39
% 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(1:12));
Fvals = CSeval(F,360,Avals);
plot(Avals,Fvals,A,D,'o')
axis([-10 370 -200 1700])
set(gca,'xTick',linspace(0,360,13))
xlabel('Ascension (Degrees)')
ylabel('Declination (minutes)')
7/21/2015
40
The Gibbs
Phenomenon
….
7/21/2015
41
% Script File: square
% Trigonometric interpolant of a square wave.
M=300
A = linspace(0,360,2*M+1)';
D(2:M)=0;D(M+2:2*M)=1;
D(1)=.5;D(M+1)=.5;D(2*M+1)=.5;
D=D';
Avals = linspace(0,360,1000)';
F = CSInterp(D(1:2*M));
Fvals = CSeval(F,360,Avals);
plot(Avals,Fvals)
7/21/2015
42
….
can be seen at
all discontinuities
7/21/2015
43
% Script File: sawtooth_wave.
clear
M=512
A = linspace(0,360,2*M+1)';
D(2:M)=(2:M)/M;D(M+2:2*M)=-1+(2:M)/M;
D(1)=0;D(M+1)=0;D(2*M+1)=0;
D=D';
Avals = linspace(0,360,1024)';
F = CSInterp(D(1:2*M));
Fvals = CSeval(F,360,Avals);
plot(Avals,Fvals)
7/21/2015
44
The overshoot
is present at
any resolution…
7/21/2015
45
…although the
interpolant is
correct at the
nodes as dx->0
7/21/2015
46
1
y
,
2
1  .5 sin x
1  x  1
Compare with
equispaced
polynomial
interpolation
7/21/2015
47
Input-output structures
•
•
•
•
•
•
x=input(‘string’)
title=input(‘Title for plot:’,’string’)
[x,y]=ginput(n)
disp(‘a 3x3 magic square’)
disp(magic(3))
fprintf(‘%5.2f \n%5.2f\n’,pi^2, exp(1))
7/21/2015
48
Writing to and reading from a file
fid = fopen('output','w')
a=linspace(1,100,100);
fprintf(fid,'%g degC = %g
degF\n',a,32+1.8*a);
fid=fopen('output','r')
x=fscanf(fid,'%g degC = %g degF');
x=reshape(x,(length(x)/2),2)
7/21/2015
49
---------------------------------------I.
A=[1 2 3;4 5 6;7 8 9];
y=x'; %transpose
x=zeros(1,n);
length(x);
x=linspace(a,b,n);
a=x(5:-1:2);
%SineTable (Help SineTable)
disp(sprintf('%2.0f %3.0f ',k, a(k)));
----------------------------------------
---------------------------------------------------------------II.
VECTOR OPERATIONS:
(1) vector: scale, add, subtract (scalar*vector)
2*[10 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: matrix*vector
(6) Tensor Products (x is array, A is array of arrays of same shape as x)
A=[f1(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) > eps*abs(y))
for term = 1:n
end
(8) size(A) gives array dimensions (array, can be used to index other arrays)
length(x) gives vector dimension
----------------------------------------------------------------
----------------------------------------------------------------
II' (Miscellaneous topics)
rem(x,n): remainder
if x == 1
(< <= == >= > ~=)
(& 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(x<=x(1))/x(1);
Summary
•
•
•
•
•
Interpolation
vanderMonde, Newton, Least Squares
2d interpolation: linear and cubic
Color figures and uint8
Smoothing of images
References
• find m-files at
www.math.unm.edu/~vageli
• Higham & Higham,
Matlab Handbook
• C.F.van Loan,
Intro to Sci.Comp. (2-3)