Chapter 3: Problem Solutions
% Problem P3_1_1
%
% Compare Locate with MyLocate
n = input('Enter n (the length of the partition): ');
m = input('Enter m (the number of evaluation points): ');
a = rand(n-1,1);
b = rand(n-1,1);
x = linspace(0,1,n)';
zvals = linspace(0,1,m)';
flops(0)
Lvals = pwLEval(a,b,x,zvals);
f1 = flops;
flops(0)
Lvals = MypwLEval(a,b,x,zvals);
f2 = flops;
clc
disp(sprintf('x = linspace(0,1,%2.0f), zvals = linspace(0,1,%2.0f)',n,m))
disp(' ')
disp(sprintf('Flops using Locate: %5.0f ',f1))
disp(sprintf('Flops using MyLocate: %5.0f ',f2))
disp(' ')
disp(' ')
x = sort(rand(n,1));
zvals = sort(rand(m,1));
flops(0)
Lvals = pwLEval(a,b,x,zvals);
f1 = flops;
flops(0)
Lvals = MypwLEval(a,b,x,zvals);
f2 = flops;
disp(sprintf('x = sort(rand(%2.0f,1)), zvals = sort(rand(%2.0f,1))',n,m))
disp(' ')
disp(sprintf('Flops using Locate: %5.0f ',f1))
disp(sprintf('Flops using MyLocate: %5.0f ',f2))
function LVals = MypwLEval(a,b,x,zvals)
% x is a column n-vector with x(1) < ... < x(n).
% zvals is a column m-vector with each component in [x(1),x(n)].
% a,b are column (n-1)-vectors
%
% LVals is a column m-vector with the property that
% LVals(j) = L(zvals(j)) for j=1:m where
% L(z) = a(i) + b(i)(z-x(i)) for x(i)<=z<=x(i+1).
m = length(zvals);
LVals = zeros(m,1);
g = 1;
for j=1:m
i = MyLocate(x,zvals(j),g);
LVals(j) = a(i) + b(i)*(zvals(j)-x(i));
g = i;
end
% Problem P3_1_2
%
% Location in uniform partitions
x = linspace(-2,2,5);
clc
disp('i x(i) z x(i+1)')
disp('-------------------------')
for z = linspace(-2,2,9)
i = LocateUniform(x(1),x(5),5,z);
disp(sprintf('%1.0f %3.0f %4.1f %3.0f',i,x(i),z,x(i+1)))
end
function i = LocateUniform(alpha,beta,n,z)
% alpha<=z<=beta and n a positive integer >=2
% x(i)<=z<=x(i+1) where x = linspace(alpha,beta,n)
h = (beta-alpha)/(n-1);
i = max(ceil((z-alpha)/h),1);
% Problem P3_1_3
%
% Applying pwLAdapt on sin(x) on [0,2pi].
[x,y] = pwLAdapt('sin',0,0,2*pi,0,.01,.01);
z = linspace(0,2*pi);
plot(z,sin(z),x,y,x,y,'o')
% Problem P3.1.4
%
% Looks at pwLAdapt on [0,1] using the function
%
% humps(x) = 1/((x-.3)^2 + .01) + 1/((x-.9)^2+.04)
%
% with delta = 0.
close all
delta = 0;
for d = [128 129];
hmin = 1/d;
figure
[x,y] = pwlAdapt('humps',0,humps(0),1,humps(1),delta,hmin);
plot(x,y,'.');
title(sprintf('hmin = 1/%3.0f Adapt n = %2.0f',d,length(x)))
end
% Problem P3_1_5
%
% Examine the effect of second derivative on pwLAdapt
% f(x) = exp(-x)sin(kx)
% f'(x) = -exp(-x)sin(kx) + 10exp(-x)cos(kx)
% f''(x) = -2exp(-x)cos(x)
clc
disp('f(x) = exp(-.3x)sin(6.1x)')
disp('Second derivative bound on [a,b] is 40exp(-.3a) assuming a>=0.')
delta = .0001;
disp(' ')
disp(' [x,y] = pwLAdapt(''MyF315'',a,fa,b,fb,.0001,0)')
disp(' ')
disp(' a b bound min(diff(x)) length(x)')
disp('------------------------------------------------------------')
for k=0:8
a = k*pi; fa = feval('MyF315',a);
b = a+pi; fb = feval('MyF315',b);
[x,y] = pwlAdapt('MyF315',a,fa,b,fb,delta,0);
n = length(x);
bound = 40*exp(-.3*a);
xdiffmin = min(diff(x));
disp(sprintf(' %6.2f % 6.2f %9.6f %8.4f %5.0f',a,b,bound,xdiffmin,n))
end
function y = MyF315(x)
y = exp(-.3*x).*sin(6.1*x);
% Problem P3.1.6
%
% Compares pwLStatic and pwLAdapt1 on [0,1] using the function
%
% humps(x) = 1/((x-.3)^2 + .01) + 1/((x-.9)^2+.04)
%
close all
% Second derivative estimate based on divided differences
z = linspace(0,1,101);
humpvals = humps(z);
M2 = max(abs(diff(humpvals,2)/(.01)^2));
for delta = [1 .5 .1 .05 .01]
figure
[x,y] = pwlStatic('humps',M2,0,1,delta);
subplot(1,2,1)
plot(x,y,'.');
title(sprintf('delta = %8.4f Static n= %2.0f',delta,length(x)))
[x,y] = pwlAdapt1('humps',0,humps(0),1/3,humps(1/3),2/3,humps(2/3),1,humps(1),delta,.001);
subplot(1,2,2)
plot(x,y,'.');
title(sprintf('delta = %8.4f Adapt n= %2.0f',delta,length(x)))
end
function [x,y] = pwLAdapt1(fname,xL,fL,xL1,fL1,xR1,fR1,xR,fR,delta,hmin)
% fname is a string that names an available function of the form y = f(u).
% xL,fL are scalars where fL = f(xL)
% xL1,fL1 are scalars where fL1 = f(xL1)
% xR,fR are scalars where fR = f(xR)
% xR1,fR1 are scalars where fR1 = f(xR1)
%
% xL1 = (2L + R)/3 and xR1 = L + 2R)/3
% delta is a positive real
% hmin is a positive real
%
% x is a column n-vector with the property that
% xL = x(1) < ... < x(n) = xR.
% Each subinterval [L,R] is either <= hmin in length or has the property
% that |f(m) - L(m)| <= delta for m = (2L+R)/3 or (L+2R)/3 and where L(x)
% is the line that connects (xL,fL) and (xR,fR).
%
% y is a column n-vector with the property that y(i) = f(x(i)).
if (xR-xL) <= hmin
% Subinterval is acceptable
x = [xL;xR];
y = [fL;fR];
else
if abs((2*fL+fR)/3 - fL1) <= delta & abs((fL+2*fR)/3 - fR1) <= delta
% Subinterval accepted.
x = [xL;xR];
y = [fL;fR];
else
% Produce left and right partitions, then synthesize.
mid = (xL+xR)/2;
fmid = feval(fname,mid);
zL = (5*xL+xR)/6;
fzL = feval(fname,zL);
zR = (xL+5*xR)/6;
fzR = feval(fname,zR);
[xLeft,yLeft] = pwLAdapt1(fname,xL,fL,zL,fzL,xL1,fL1,mid,fmid,delta,hmin);
[xRight,yRight] = pwLAdapt1(fname,mid,fmid,xR1,fR1,zR,fzR,xR,fR,delta,hmin);
x = [ xLeft;xRight(2:length(xRight))];
y = [ yLeft;yRight(2:length(yRight))];
end
end
% Problem P3_1_7
%
% Applying pwLAdapt to sqrt
close all
[x,y] = pwLAdapt('sqrt',0,0,1,1,.001,.01);
z = linspace(0,1);
plot(z,sqrt(z),x,y,x,y,'o')
axis([-.1 1.1 -.2 1.2])
hold on
plot(x,zeros(size(x)),'.')
hold off
% Problem P3_1_8
%
% Examine pwLAdapt with feval max on [0,1] using the function
%
% humps(x) = 1/((x-.3)^2 + .01) + 1/((x-.9)^2+.04)
%
close all
eMax = 50;
for delta = [1 .5 .1 .05 .01]
[x,y,eTaken] = pwlAdapt2('humps',0,humps(0),1,humps(1),delta,.001,eMax);
figure
plot(x,y,'.');
axis([0 1 0 100])
title(sprintf('delta = %8.4f x(%1.0f) = %5.3f',delta,length(x),x(length(x))))
pause
end
function [x,y,eTaken] = pwLAdapt2(fname,xL,fL,xR,fR,delta,hmin,eMax)
% fname is a string that names an available function of the form y = f(u).
% xL is real and fL = f(xL)
% xR is real and fR = f(xR)
% delta is a positive real
% hmin is a positive real
% eMax is a positive integer representing the maximum number of fevals allowed.
%
% x is a column n-vector with the property that
% xL = x(1) < ... < x(n) = xR. Each subinterval
% is either <= hmin in length or has the property
% that at its midpoint m, |f(m) - L(m)| <= delta
% where L(x) is the line that connects (xR,fR).
%
% y is a column n-vector with the property that y(i) = f(x(i)).
% eTaken is the number of f-evaluations required.
if (xR-xL) <= hmin
% Subinterval is acceptable
x = [xL;xR];
y = [fL;fR];
eTaken = 0
else
mid = (xL+xR)/2;
fmid = feval(fname,mid);
eTaken = 1;
eMax = eMax-1;
if (abs(((fL+fR)/2) - fmid) <= delta )
% Subinterval accepted.
x = [xL;xR];
y = [fL;fR];
else
% Produce left and right partitions, then synthesize.
if eMax ==0
x = xL;
y = fL;
return
end
[xLeft,yLeft,eTakenL] = pwLAdapt2(fname,xL,fL,mid,fmid,delta,hmin,eMax);
eTaken = eTaken + eTakenL;
eMax = eMax - eTakenL;
if eMax == 0
x = xLeft;
y = yLeft;
return
end
[xRight,yRight,eTakenR] = pwLAdapt2(fname,mid,fmid,xR,fR,delta,hmin,eMax);
eTaken = eTaken + eTakenR;
x = [ xLeft;xRight(2:length(xRight))];
y = [ yLeft;yRight(2:length(yRight))];
end
end
% Problem P3_1_9
%
% Examine pwLAdapt with saved f-evals on [0,1] using the function
%
% humps(x) = 1/((x-.3)^2 + .01) + 1/((x-.9)^2+.04)
%
close all
[x,y,xunused,yunused] = pwlAdapt3('humps',0,humps(0),1,humps(1),1,.001,[],[]);
plot(x,y,'.',xunused,yunused,'o');
title('o = the saved f-evaluations')
z
% Problem P3_2_1
%
% Apply pwcStatic to exp(-2x)sin(10pi*x) on [0,1].
close all
alpha = 0;
beta = 1;
delta = .001;
z = linspace(0,1)';
fz = feval('MyF321',z);
for M4 = [100 1000 10000 100000 1000000]
[a,b,c,d,x] = pwCStatic('MyF321','dMyF321',M4,alpha,beta,delta);
Cvals = pwCEval(a,b,c,d,x,z);
err = max(abs(Cvals-fz));
figure
plot(z,fz,z,Cvals,x,feval('MyF321',x),'o')
title(sprintf('M4 = %3.0e n = %2.0f delta = %2.0e error = %2.0e',M4,length(x),delta,err))
end
function [a,b,c,d,x] = pwCStatic(fname,fpname,M4,alpha,beta,delta)
% fname is a string that names an available function f(x).
% Assume that f can take vector arguments.
% M4 is an upper bound for|f''''(x)| on [alpha,beta].
% alpha and beta are scalars with alpha<beta
% delta is a positive scalar
%
% x is a column n-vector that defines the partition.
% a,b,c,d is a column (n-1)-vectors with the property that C(x(i)) = f(x(i)),
% i=1:n. The piecewise cubic Hermite interpolant C(x) of this
% data satisfies |C(x) - f(x)| <= delta, where
% x = linspace(alpha,beta,n).
n = max(2,ceil(1 + ((beta-alpha)*M4/(384*delta))^.25) );
x = linspace(alpha,beta,n)';
y = feval(fname,x);
s = feval(fpname,x);
a = y(1:n-1);
b = s(1:n-1);
Dx = (beta-alpha)/(n-1);
Dy = diff(y);
yp = Dy/Dx;
c = (yp - s(1:n-1))/Dx;
d = (s(1:n-1) - 2*yp + s(2:n))/Dx^2;
function y = dMyF321(x)
x2 = 10*pi*x;
y = exp(-2*x).*( -2*sin(x2) + (10*pi)*cos(x2));
% Problem P3_2_2
%
% Compare pwLAdapt and pwCAdapt on f(x) = sqrt(x) on [.001,9]
hmin = .001;
L = .001; fL = sqrt(L); sL = .5/sqrt(L);
R = 9.000; fR = sqrt(R); sR = .5/sqrt(R);
clc
disp(' pwL pwC')
disp(' delta length(x) length(x)')
disp('--------------------------------------')
for delta = [.1 .01 .001 .0001 .00001]
[x,y] = pwLAdapt('sqrt',L,fL,R,fR,delta,hmin);
nL = length(x);
[x,y,s] = pwCAdapt('sqrt','dMyF322',L,fL,sL,R,fR,sR,delta,hmin);
nC = length(x);
disp(sprintf(' %7.5f %3.0f %3.0f',delta,nL,nC))
end
function [x,y,s] = pwCAdapt(fname,fpname,xL,fL,sL,xR,fR,sR,delta,hmin)
% fname is a string that names an available function of the form y = f(u).
% fpname is a string that names a function that is the derivative of f.
% xL is real and fL = f(xL) and sL = f'(xL).
% xR is real and fR = f(xR) and sR = f'(xR).
% delta is a positive real
% hmin is a positive real
%
% x is a column n-vector with the property that
% xL = x(1) < ... < x(n) = xR. Each subinterval
% is either <= hmin in length or has the property
% that at its midpoint m, |f(m) - q(m)| <= delta
% where q(x) is the cubic hermite interpolant of
% (xL,fL,sL) and (xR,fR,sR).
%
% y is a column n-vector with the property that
% y(i) = f(x(i)).
%
% sc is a column n-vector with the property that
% s(i) = f'(x(i)).
if (xR-xL) <= hmin
% Subinterval is acceptable
x = [xL;xR];
y = [fL;fR];
s = [sL;sR];
else
mid = (xL+xR)/2;
fmid = feval(fname,mid);
% Compute the cubic hermite interpolant and evaluate at the midpoint:
[alfa,beta,gamma,eta] = HCubic(xL,fL,sL,xR,fR,sR);
qeval = pwCEval(alfa,beta,gamma,eta,[xL;xR],mid);
if abs(qeval - fmid) <= delta
% Subinterval accepted.
x = [xL;xR];
y = [fL;fR];
s = [sL;sR];
else
% Produce left and right partitions and synthesize.
smid = feval(fpname,mid);
[xLeft,yLeft,sLeft] = pwCAdapt(fname,fpname,xL,fL,sL,mid,fmid,smid,delta,hmin);
[xRight,yRight,sRight] = pwCAdapt(fname,fpname,mid,fmid,smid,xR,fR,sR,delta,hmin);
x = [ xLeft;xRight(2:length(xRight))];
y = [ yLeft;yRight(2:length(yRight))];
s = [ sLeft;sRight(2:length(sRight))];
end
end
function y = dMyF322(x)
y = .5./sqrt(x);
% Problem P3_2_3
%
% Variable subinterval piecewise cubic hermite interpolation.
close all
[x,y] = pwCHExp(0,5,.0001);
[a,b,c,d] = pwC(x,y,-y);
z = linspace(0,5)';
Cvals = pwCEval(a,b,c,d,x,z);
plot(z,exp(-z),z,Cvals,x,y,'o')
title('Piecewise cubic hermite interpolant of exp(-x)')
function [x,y] = pwCHExp(a,b,tol)
% a,b are real scalars with a < b.
% tol is a positive real scalar.
%
% x,y is a column n-vectors where a = x(1) < x(2) < ... < x(n) = b
% and y(i) = exp(-x(i)), i=1:n.
% The partition is chosen so that if C(z) is the piecewise cubic
% hermite interpolant of exp(-z) on this partition,
% then |C(z) - exp(-z)| <= tol for all z in [a,b].
x = a;
y = exp(-a);
n = 1;
while x(n) < b
[R,fR] = stretch(x(n),y(n),tol);
if R < b
x = [x;R];
y = [y;fR];
else
x = [x;b];
y = [y;exp(-b)];
end
n = n+1;
end
function [R,fR] = stretch(L,fL,tol);
% L,fL are scalars that satisfy fL = exp(-L)
% tol is a positive real.
%
% R,fR are scalars that satisfy fR = exp(-R) with the
% property that if q(z) is the cubic hermite interpolant
% of exp(-z) at z=L and z=R, then |q(z) - exp(-z)| <= tol on [L,R].
h = sqrt(sqrt((384*tol/fL)));
R = L + h;
fR = exp(-R);
z
z
% Problem P3_3_1
%
% Not-a-knot spline interpolation of f(x) = x^3.
x = sort(rand(4,1));
y = x.^3;
[a,b,c,d] = CubicSpline(x,y);
z = linspace(x(1),x(4));
Cvals = pwCEval(a,b,c,d,x,z);
z0 = linspace(x(1),x(4),20);
plot(z,Cvals,z0,z0.^3,'o')
title('Not-a-knot spline interpolant of x^3')
% Problem P3_3_2
%
% 4-point not-a-knot spline interpolant
close all
x = sort(rand(4,1));
y = rand(4,1);
[a,b,c,d] = CubicSpline(x,y);
z = linspace(x(1),x(4))';
Cvals = pwCEval(a,b,c,d,x,z);
a = InterpV(x,y);
z0 = linspace(x(1),x(4),20)';
pvals = HornerV(a,z0);
plot(z,Cvals,x,y,'*',z0,pvals,'o')
title('n=4 not-a-knot is just cubic interpolation.')
% Problem P3_3_3
%
% Natural spline interpolant of x^3.
close all
x = [-3;-1;1;3];
y = x.^3;
[a,b,c,d] = CubicSpline(x,y,2,0,0);
S0 = pwCEval(a,b,c,d,x,0);
z = linspace(-3,3);
Cvals = pwCEval(a,b,c,d,x,z);
plot(z,Cvals,x,y,'o')
title(sprintf('S(0) = %8.4e',S0))
% Problem P3_3_4
%
% Illustrate spline-under-tension idea. As sigma increases,
% the approximant looks more and more like a piecewise linear
% function.
close all
n = 10;
x = linspace(0,1,n)';
y = exp(-2*x).*sin(10*pi*x);
s = -2*exp(-2*x).*sin(10*pi*x) + 10*pi*exp(-2*x).*cos(10*pi*x);
z = linspace(0,1)';
fz = exp(-2*z).*sin(10*pi*z);
plot(z,fz,x,y,'o')
sigma = 100;
[a,b,c,d] = Fit(x,y,s,sigma);
hold on
for i=1:n-1
xvals = linspace(x(i),x(i+1),100)';
u = exp(sigma*(xvals-x(i)));
yvals = a(i) + b(i)*(xvals-x(i)) + c(i)*u + d(i)./u;
plot(xvals,yvals,'.')
end
hold off
title('n=10, f(x) = exp(-2x)sin(10*pi*x), sigma = 100')
for sigma = [ 10 100 250];
[a,b,c,d] = Fit(x,y,s,sigma);
figure
hold on
for i=1:n-1
xvals = linspace(x(i),x(i+1),100)';
u = exp(sigma*(xvals-x(i)));
yvals = a(i) + b(i)*(xvals-x(i)) + c(i)*u + d(i)./u;
plot(xvals,yvals,x,y,x,y,'o')
end
title(sprintf('PWL fit and the tension spline fit with sigma = %6.1f',sigma))
hold off
end
function [a,b,c,d] = Fit(x,y,s,sigma)
% x,y,s are column n-vectors
% sigma is a positive scalar
%
% a,b,c,d are column (n-1)-vectors so that if
%
% g(z) = a(i) + b(i)(z-x(i)) + c(i)exp(sigma(z-x(i))) + d(i)exp(-sigma(z-x(i)))
%
% then g(x(i))= y(i), g'(x(i)) = s(i), g(x(i+1)) = y(i+1), and g'(x(i+1)) = s(i+1)
% for i=1:n.
n = length(x);
a = zeros(n-1,1);
b = zeros(n-1,1);
c = zeros(n-1,1);
d = zeros(n-1,1);
for i=1:n-1
del = x(i+1)-x(i);
f = exp(sigma*del);
M = [1 0 1 1; ...
1 del f 1/f; ...
0 1 sigma -sigma; ...
0 1 sigma*f -sigma/f];
rhs = [y(i);y(i+1);s(i);s(i+1)];
u = M\rhs;
a(i) = u(1);
b(i) = u(2);
c(i) = u(3);
d(i) = u(4);
end
% Problem P3_3_5
%
% Illustrates CubicSpline1 with various end conditions,
% some not so good.
close all
xvals = linspace(0,4*pi,100);
yvals = sin(xvals);
for n = 4:3:16
x = linspace(0,4*pi,n)';
y = sin(x);
[a,b,c,d] = CubicSpline1(x,y,1);
svals = pwCEval(a,b,c,d,x,xvals);
plot(xvals,yvals,xvals,svals)
title(sprintf('Spline interpolant of sin(x) using qL'' and qR''. (%2.0f subintervals)',n-1))
pause
end
for n = 4:3:16
x = linspace(0,4*pi,n)';
y = sin(x);
[a,b,c,d] = CubicSpline1(x,y,2);
svals = pwCEval(a,b,c,d,x,xvals);
plot(xvals,yvals,xvals,svals)
title(sprintf('Spline interpolant of sin(x) using qL'''' and qR''''. (%2.0f subintervals)',n-1))
pause
end
for n = 4:3:16
x = linspace(0,4*pi,n)';
y = sin(x);
[a,b,c,d] = CubicSpline1(x,y,1);
svals = pwCEval(a,b,c,d,x,xvals);
plot(xvals,yvals,xvals,svals)
title(sprintf('Spline interp. of sin(x) using derivatives of end interpolants. (%2.0f subintervals)',n-1))
pause
end
for n = 4:3:16
x = linspace(0,4*pi,n)';
y = sin(x);
[a,b,c,d] = CubicSpline1(x,y,1,1,1);
svals = pwCEval(a,b,c,d,x,xvals);
plot(xvals,yvals,xvals,svals)
title(sprintf('''Good'' Complete spline interpolant of sin(x). (%2.0f subintervals)',n-1))
pause
end
for n = 4:3:16
x = linspace(0,4*pi,n)';
y = sin(x);
[a,b,c,d] = CubicSpline1(x,y,1,-1,-1);
svals = pwCEval(a,b,c,d,x,xvals);
plot(xvals,yvals,xvals,svals)
title(sprintf('''Bad'' Complete spline interpolant of sin(x). (%2.0f subintervals)',n-1))
pause
end
for n = 4:3:16
x = linspace(0,4*pi,n)';
y = sin(x);
[a,b,c,d] = CubicSpline1(x,y,2,0,0);
svals = pwCEval(a,b,c,d,x,xvals);
plot(xvals,yvals,xvals,svals)
title(sprintf('Natural spline interpolant of sin(x). (%2.0f subintervals)',n-1))
pause
end
for n = 4:3:16
x = linspace(0,4*pi,n)';
y = sin(x);
[a,b,c,d] = CubicSpline1(x,y);
svals = pwCEval(a,b,c,d,x,xvals);
plot(xvals,yvals,xvals,svals)
title(sprintf('Not-a-Knot spline interpolant of sin(x). (%2.0f subintervals)',n-1))
pause
end
function [a,b,c,d] = CubicSpline1(x,y,derivative,muL,muR)
% x,y are column n-vectors. n >= 4 and x(1) < ... x(n)
% derivative is the order of the spline's derivative that are
% used in the end conditions (= 1 or 2)
% muL,muR are the values of the specified derivative at the
% left and right endpoint.
%
% a,b,c,d are column (n-1)-vectors that define the spline.
% On [x(i),x(i+1)], the spline S(z) is specified by the cubic
%
% a(i) + b(i)(z-x(i)) + c(i)(z-x(i))^2 + d(i)(z-x(i))^2(z-x(i+1).
%
% Usage:
% [a,b,c,d] = CubicSpline1(x,y,1,muL,muR) S'(x(1)) = muL, S'(x(n)) = muR
% [a,b,c,d] = CubicSpline1(x,y,1) S'(x(1)) = qL'(x(1)), S'(x(n)) = qR'(x(n))
% [a,b,c,d] = CubicSpline1(x,y,2,muL,muR) S''(x(1)) = muL, S''(x(n)) = muR
% [a,b,c,d] = CubicSpline1(x,y,2) S''(x(1)) = qL''(x(1)), S''(x(n)) = qR''(x(n))
% [a,b,c,d] = CubicSpline1(x,y) S'''(z) continuous at x(2) and x(n-1)
% Here qL = cubic interpolant of (x(i),y(i)), i=1:4
% qR = cubic interpolant of (x(i),y(i)), i=n-3:n
if nargin == 3
% need qL and qR
n = length(x);
cL = InterpN(x(1:4),y(1:4));
cR = InterpN(x(n:-1:n-3),y(n:-1:n-3));
end
if nargin == 2
[a,b,c,d] = CubicSpline(x,y);
else
if derivative == 1
if nargin == 3
% Differentiate the Newton interpolants and evaluate.
muL = cL(2) + cL(3)*(x(1)-x(2)) + cL(4)*((x(1)-x(2)) *(x(1)-x(3)));
muR = cR(2) + cR(3)*(x(n)-x(n-1)) + cR(4)*((x(n)-x(n-1))*(x(n)-x(n-2)));
end
[a,b,c,d] = CubicSpline(x,y,1,muL,muR);
else
if nargin ==3
% Differentiate the Newton interpolants twice and evaluate.
muL = 2*cL(3) + 2*cL(4)*((x(1)-x(2)) + (x(1)-x(3)));
muR = 2*cR(3) + 2*cR(4)*((x(n)-x(n-1)) + (x(n)-x(n-2)));
end
[a,b,c,d] = CubicSpline(x,y,2,muL,muR);
end
end
% Problem P3_3_6
%
% Integral of the splines second derivative.
clc
disp(' n Natural Not-a-Knot Bad Complete')
disp('--------------------------------------------------')
for n = [10 20 40 80 160 320 ]
x = linspace(0,1,n)';
y = 1./(1+25*x.^2);
[a,b,c,d] = CubicSpline(x,y,2,0,0);
e1 = Energy(ConvertS(a,b,c,d,x));
[a,b,c,d] = CubicSpline(x,y);
e2 = Energy(ConvertS(a,b,c,d,x));
[a,b,c,d] = CubicSpline(x,y,1,30,-15);
e3 = Energy(ConvertS(a,b,c,d,x));
disp(sprintf(' %3.0f %8.5e % 8.5e %8.5e',n,e1,e2,e3))
end
function e = Energy(S)
% S is the pp representation of a spline
% e is the integral from x(1) to x(n) of the
% square of the second derivative of the spline
[x,rho,L,k] = unmkpp(S);
% The second derivative of the ith cubic is 2*rho(i,2) + 6*rho(i,1)(x-x(i)).
% The integral from x(i) to x(i+1) of [2*rho(i,2) + 6*rho(i,1)(x-x(i))]^2
% is
% 4* rho(i,2)^2 * del + 12*rho(i,2)*rho(i,1)*del^2 + 12*rho(i,1)^2*del^3
del = (x(2:L+1)-x(1:L))';
r1 = 4*rho(:,2).^2;
r2 = 12*rho(:,2).*rho(:,1);
r3 = 12*rho(:,1).^2;
e = sum(((del.*r3 + r2).*del + r1).*del);
% Problem P3_3_7
%
% Max 3rd derivative continuity.
clc
n = input('Enter the number of points to interpolate: ');
x = linspace(0,pi,n)';
y = exp(-x).*sin(x);
S = spline(x,y);
d3 = MaxJump(S)
function d3 = MaxJump(S);
% S is the pp representation of a cubic spline.
% d3 is the absolute value of the largest 3rd derivative
% jump at any knot.
[x,rho,L,k] = unmkpp(S);
% 3rd derivative of the i-th local cubic is 6*rho(i,1)
if L==1
d3 = 0;
else
d3 = 6*max(abs(rho(2:L,1)-rho(1:L-1,1)));
end
% Problem P3_3_8 % % Different representations for piecewise cubics. x = [0 .1 .3 .5 .6 1]'; y = sin(2*pi*x); % Generate not-a-knot spline 2 ways: [a,b,c,d] = CubicSpline(x,y); S1 = Spline(x,y); [x1,rho1,L1,k1] = unmkpp(S1); % Convertthe book spline tothe pp representation and compare % with S1. S2 = ConvertS(a,b,c,d,x); [x2,rho2,L2,k2] = unmkpp(S2); x_error = x1-x2 rho_error = rho1-rho2 L_error = L1-L2 k_error = k1-k2 function S = ConvertS(a,b,c,d,x) % x is a column n-vector with x(1) < ... < x(n) % a,b,c, and d are column n-1 vectors. % % S is the pp-representation of the piecewise cubic whose % ith local cubic is defined by % % a(i) + b(i)*(x-x(i)) + c(i)*(x-x(i))^2 + d(i)*(x-x(i))^2*(x-x(i+1)) % % Note that the i-th local cubic is also prescribed by % a(i) + b(i)*(x-x(i)) + (c(i)+d(i)(x(i)-x(i+1))*(x-x(i))^2 + d(i)*(x-x(i))^3 rho = [d c-d.*diff(x) b a]; S = mkpp(x,rho);
% Problem P3_3_9 % % Illustrate a 1-knot natural spline. y = [-1 ; 3 ; -2]; [a,b,c,d] = SmallSpline(y); zL = linspace(-1,0)'; CL = ((d(1)*zL + c(1)).*zL+b(1)).*zL+a(1); zR = linspace(0,1)'; CR = ((d(2)*zR + c(2)).*zR+b(2)).*zR+a(2); plot([zL;zR],[CL;CR],[-1;0;1],y,'o') function [a,b,c,d] = SmallSpline(y) % y is 3-vector. % % a,b,c,d are column 2-vectors with the property that if % % S(x) = a(1) + b(1)x + c(1)x^2 + d(1)x^3 on [-1,0] % and % S(x) = a(2) + b(2)x + c(2)x^2 + d(2)x^3 on [0,1] % then % (a) S(-1) = y(1), S(0) = y(2), S(1) = y(3), % (b) S''(-1) = S''(1) = 0 % (c) S, S', and S'' are continuous on [-1,1] % rhs = zeros(8,1); % S(-1) = y(1) M = [ 1 -1 1 -1 0 0 0 0 ]; rhs(1) = y(1); % S(0) = y(2) M = [M ; 1 0 0 0 0 0 0 0]; rhs(2) = y(2); % S(1) = y(3) M = [M; 0 0 0 0 1 1 1 1]; rhs(3) = y(3); % S''(-1) = 0 M = [M; 0 0 2 -6 0 0 0 0]; rhs(4) = 0; % S''(1) = 0 M = [M; 0 0 0 0 0 0 2 6]; rhs(5) = 0; % S continuous at x = 0 M = [M; 1 0 0 0 -1 0 0 0 ]; rhs(6) = 0; % S' continuous at x = 0 M = [M; 0 1 0 0 0 -1 0 0 ]; rhs(7) = 0; % S'' continuous at x = 0 M = [M; 0 0 1 0 0 0 -1 0 ]; rhs(8) = 0; u = M\rhs; a = u([1 5]); b = u([2 6]); c = u([3 7]); d = u([4 8]);
% Problem P3_3_10
%
% Spline interpolation in the plane.
close all
theta = linspace(0,2*pi)';
xvals = 5*cos(theta);
yvals = 3*sin(theta);
for n = 5:4:13
theta = linspace(0,2*pi,n)';
x = 5*cos(theta);
y = 3*sin(theta);
[xi,yi] = SplineInPlane(x,y,100);
figure
plot(xvals,yvals,xi,yi,x,y,'o')
axis('equal')
pause
end
function [xi,yi] = SplineInPlane(x,y,m)
% x and y are columnn-vectors.
% xi and yi are pp-representations of not-a-knot splines with the
% property that they interpolate x(i) and y(i) respectively
% at their knots.
n = length(x);
d = sqrt( (x(2:n)-x(1:n-1)).^2 + (y(2:n)-y(1:n-1)).^2 );
t = [0; cumsum(d)];
ti = linspace(0,t(n),m)';
xi = spline(t,x,ti);
yi = spline(t,y,ti);
z
z
z
z