Chapter 3: Problem Solutions

Problems Home

 P3.1.1 P3.2.1 P3.3.1 P3.1.2 P3.2.2 P3.3.2 P3.1.3 P3.2.3 P3.3.3 P3.1.4 P3.2.4 P3.3.4 P3.1.5 P3.2.5 P3.3.5 P3.1.6 P3.3.6 P3.1.7 P3.3.7 P3.1.8 P3.3.8 P3.1.9 P3.3.9 P3.1.10 P3.3.10 P3.3.11 P3.3.12 P3.3.13 P3.3.14

P3.1.1

```% 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

```

P3.1.2

```% 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);
```

P3.1.3

```% Problem P3_1_3
%
% Applying pwLAdapt on sin(x) on [0,2pi].

z = linspace(0,2*pi);
plot(z,sin(z),x,y,x,y,'o')
```

P3.1.4

```% 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
plot(x,y,'.');
title(sprintf('hmin = 1/%3.0f  Adapt n = %2.0f',d,length(x)))
end
```

P3.1.5

```% 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(' ')
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);
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);
```

P3.1.6

```% 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)))
subplot(1,2,2)
plot(x,y,'.');
title(sprintf('delta = %8.4f  Adapt n= %2.0f',delta,length(x)))
end

% 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);
x = [ xLeft;xRight(2:length(xRight))];
y = [ yLeft;yRight(2:length(yRight))];
end
end
```

P3.1.7

```% Problem P3_1_7
%

close all
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
```

P3.1.8

```% 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]
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

% 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
eTaken = eTaken + eTakenL;
eMax = eMax - eTakenL;
if eMax == 0
x = xLeft;
y = yLeft;
return
end

eTaken = eTaken + eTakenR;
x = [ xLeft;xRight(2:length(xRight))];
y = [ yLeft;yRight(2:length(yRight))];
end
end
```

P3.1.9

```% 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
plot(x,y,'.',xunused,yunused,'o');
title('o = the saved f-evaluations')
```

P3.1.10

```z
```

P3.2.1

```% 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));

```

P3.2.2

```% Problem P3_2_2
%

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]

nL = length(x);
nC = length(x);
disp(sprintf(' %7.5f      %3.0f            %3.0f',delta,nL,nC))
end

% 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);

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);

```

P3.2.3

```% 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);
```

P3.2.4

```z
```

P3.2.5

```z
```

P3.3.1

```% 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')
```

P3.3.2

```% 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.')
```

P3.3.3

```% 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))
```

P3.3.4

```% 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
```

P3.3.5

```% 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
```

P3.3.6

```% 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);
```

P3.3.7

```% 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

```

P3.3.8

```% 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);
```

P3.3.9

```% 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]);

```

P3.3.10

```% 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);

```

P3.3.11

```z
```

P3.3.12

```z
```

P3.3.13

```z
```

P3.3.14

```z
```