Math 4401 Numerical Methods

Notebook Author/Course Instructor: Barry McQuarrie
Date: Sept 2016
Lecture topic: Cubic Splines and Approximating Points

Natural Cubic Splines

Brute Force Method of Solution

This is a brute force method of solution, that solves the system of equations that arises. It is quite cumbersome, and is why the algorithm that follows is so useful. Notice that I am using Mathematica’s Solve command here, since this is ultimately not how we will implement this. It was kind of finicky to get all the underscores and Mathematica syntax correct, but it works!

Here is a data set that we want to be able to interpolate on:

Clear[data]
data={{0.9,1.3},{1.3,1.5},{1.9,1.85},{2.1,2.1},{2.6,2.6},{3.0,2.7},{3.9,2.4},{4.4,2.15}};

dataplot=ListPlot[data,PlotStyle->{Red,PointSize[0.02]},Frame->True,ImageSize->{500,Automatic},BaseStyle->20]

Clear[a,b,c,d,x,xp,yp,eqlist,varlist,SS,SSd,SSdd,S]

n=Length[data];

(*data points*)
Do[
xp[i]=data[[i,1]];
yp[i]=data[[i,2]];
,{i,1,n}]

(*define the cubic spline and derivatives in each section*)
SS[i_,x_]:=a[i]+b[i](x-xp[i])+c[i](((x-xp[i]))^2)+d[i](((x-xp[i]))^3)
SSd[i_,x_]:=b[i]+2c[i](x-xp[i])+3d[i](((x-xp[i]))^2)
SSdd[i_,x_]:=2c[i]+6d[i](x-xp[i])

(*construct a list of the variables we need to solve for*)
varlist={};
Do[varlist=Append[varlist,{a[i],b[i],c[i],d[i]}]
,{i,1,n-1}]
varlist=Flatten[varlist];

(*construct the equations that must be satisfied*)
eqlist={};
Do[
eqlist=Append[eqlist,SS[i,xp[i]]==yp[i]];
eqlist=Append[eqlist,SS[i,xp[i+1]]==yp[i+1]];
If[i≠1,
{eqlist=Append[eqlist,SSd[i-1,xp[i]]==SSd[i,xp[i]]],
eqlist=Append[eqlist,SSdd[i-1,xp[i]]==SSdd[i,xp[i]]]}
];
,{i,1,n-1}]

(*free cubic spline*)
eqlist=Append[eqlist,SSdd[1,xp[1]]==0];
eqlist=Append[eqlist,SSdd[n-1,xp[n]]==0];

(*solve the equations*)
sol=Solve[eqlist,varlist];

(*assign the variables to the cubic functions we defined*)
Do[
SS[i,x_]=SS[i,x]/.Flatten[sol];
SSd[i,x_]=SSd[i,x]/.Flatten[sol];
SSdd[i,x_]=SSdd[i,x]/.Flatten[sol];
,{i,1,n-1}]

(*knit all the cubic functions together into the cubic spline and its derivatives*)
S[x_]:=Do[ If[xp[i]≤x<xp[i+1],Return[SS[i,x]]],{i,1,n-1}]
Sd[x_]:=Do[ If[xp[i]≤x<xp[i+1],Return[SSd[i,x]]],{i,1,n-1}]
Sdd[x_]:=Do[ If[xp[i]≤x<xp[i+1],Return[SSdd[i,x]]],{i,1,n-1}]

Here are the plots of the cubic functions we have created:

splineplots={};
Do[
splineplots=Append[splineplots,Plot[{SS[i,x]},{x,xp[i],xp[i+1]},PlotRange->{{xp[1],xp[n]},{0,3}},Frame->True,Axes->False,PlotStyle->{Hue[i/n],Thick}]];
,{i,1,n-1}]

Show[dataplot,splineplots,ImageSize->{600,Automatic},BaseStyle->20]

We can work with the entire cubic spline using the function S[x]:

Plot[S[x],{x,xp[1],xp[n]},PlotStyle->{Blue,Thick},Frame->True,FrameLabel->{"x","Cubic Spline"},PlotRange->{{xp[1],xp[n]},{0,3}},ImageSize->{600,Automatic}]

The first derivative is smooth:

Plot[Sd[x],{x,xp[1],xp[n]},FrameLabel->{"x","First Derivative"},Frame->True,Axes->False,PlotStyle->{Thick,Blue},ImageSize->{600,Automatic}]

The second derivative is continuous:

Plot[Sdd[x],{x,xp[1],xp[n]},FrameLabel->{"x","Second Derivative"},Frame->True,Axes->False,PlotStyle->{Thick,Blue},ImageSize->{600,Automatic}]

The Natural Cubic Spline Conditions: Solve for c, then define a,b,d.

Let’s examine the structure of the equations we are solving, and see if we can arrive at a better algorithm than just brute force.

Clear[i,f,a,b,c,d,x,xp,yp,eqlist,varlist,SS,SSd,SSdd,S]

SS[i_,x_]=a[i]+b[i](x-xp[i])+c[i](((x-xp[i]))^2)+d[i](((x-xp[i]))^3);
SSd[i_,x_]=b[i]+2c[i](x-xp[i])+3d[i](((x-xp[i]))^2);
SSdd[i_,x_]=2c[i]+6d[i](x-xp[i]);

Here are the four conditions (we will add the free cubic spline conditions later):

eq1=SS[i,xp[i]]==yp[i]

eq2=SS[i,xp[i+1]]==yp[i+1]

eq3=SSd[i-1,xp[i]]==SSd[i,xp[i]]

eq4=SSdd[i-1,xp[i]]==SSdd[i,xp[i]]

If we solve Eq (1) for a, we get the following:

Solve[eq1,a[i]]

Simplify

a[i]=yp[i]

where 
i=1,...,n-1

If we solve Eq (4) for d, we get the following:

Solve[eq4,d[i-1]]
%/.i->i+1

Simplify

d[i]=((c[i+1]-c[i])/(3 h[i]))

where 
h[i]=xp[i+1]-xp[i]
i=1,...,n-1
(Note: The i=n-1 is correct if we define c[n](defined below, but at the moment it is not in the list of our unknown parameters).

If we solve Eq (2) for b, we get the following:

Solve[eq2,b[i]]

b[i]=((Δ[i])/(h[i]))-c[i] h[i]-d[i] (((h[i]))^2)
=((Δ[i])/(h[i]))-c[i] h[i]-((c[i+1]-c[i])/(3 h[i])) (((h[i]))^2)
=((Δ[i])/(h[i]))-c[i] h[i]-((c[i+1]-c[i])/3) h[i]
=((Δ[i])/(h[i]))-((h[i])/3)( 2c[i]+c[i+1])

where 
h[i]=xp[i+1]-xp[i]
Δ[i]=yp[i+1]-yp[i]
i=1,...,n-1

So the system can be reduced to determing the values of c[i] using Eq (3), and then getting the other parameters from the relationships:

a[i]=yp[i]

b[i]=((Δ[i])/(h[i]))-((h[i])/3)( 2c[i]+c[i+1])

d[i]=((c[i+1]-c[i])/(3 h[i]))

where 
h[i]=xp[i+1]-xp[i]
Δ[i]=yp[i+1]-yp[i]
i=1,...,n-1

Let’s get the equation that c[i] satisfies:

eq3

eq3=b[-1+i]+2 c[-1+i] (h[i-1])+3 d[-1+i] (((h[i-1]))^2)==b[i]

eq3/.{b[i]->((Δ[i])/(h[i]))-((h[i])/3)( 2c[i]+c[i+1]),b[i-1]->((Δ[i-1])/(h[i-1]))-((h[i-1])/3)( 2c[i-1]+c[i]),d[i-1]->((c[i]-c[i-1])/(3 h[i-1]))}
Collect[%,c[i]]

c[i-1] h[i-1]+c[i](2 h[i-1]+2 h[i])+ c[i+1] h[i]==3(((Δ[i])/(h[i]))-((Δ[i-1])/(h[i-1])))

Now, we can use the free spline boundary conditions to get the final equations:

SSdd[1,xp[1]]==0
SSdd[n-1,xp[n]]==0

Simplify the latter (and this is where the c[n] gets defined):

2 c[-1+n]+6 ((c[n]-c[n-1])/(3 h[n-1])) (h[n-1])==0
2 c[-1+n]+2 c[n]-2c[n-1]==0
2 c[n]==0

So, the equations that c satisfy are (this is equation 3.24 in the text):

c[1]=0
c[i-1] h[i-1]+c[i](2 h[i-1]+2 h[i])+ c[i+1] h[i]==3(((Δ[i])/(h[i]))-((Δ[i-1])/(h[i-1]))), i=2,3,4,...n-1
c[n]=0

So we can solve for the quantities in a free cubic spline through the n points ((Subscript[x, i]),(Subscript[y, i])) i=1,2,...,n by the following algorithm:

For i=1 to n-1,
a[i]=yp[i]
h[i]=xp[i+1]-xp[i]
Δ[i]=yp[i+1]-yp[i]

Solve the system of equations for c[i],
c[1]=0
c[i-1] h[i-1]+c[i](2 h[i-1]+2 h[i])+ c[i+1] h[i]==3(((Δ[i])/(h[i]))-((Δ[i-1])/(h[i-1]))), i=2,3,4,...n-1
c[n]=0

For i=1 to n-1,
b[i]=((Δ[i])/(h[i]))-((h[i])/3)( 2c[i]+c[i+1])

d[i]=((c[i+1]-c[i])/(3 h[i]))

S[i,x]=a[i]+b[i](x-xp[i])+c[i](((x-xp[i]))^2)+d[i](((x-xp[i]))^3)

Now that we have this process laid out, we need to focus on how to solve the tridiagonal system involving the c[i].

LU Decomposition to solve for the parameters c

Solving Linear Systems of Equations

The system of equations we need to solve is a tridiagonal system.

We can solve a linear system of equations for an unknown vector x in a variety of ways.

A x = B
A is an n × n matrix,
x, B  are n vectors.

One way is through Gaussian elimination (from linear algebra) with backward substitution. There is also Cramer's rule, and through computation of (A^(-1)). Another way is through decomposing the matrix A into the product of a lower triangular matrix and an upper triangular matrix. Once this has been accomplished, we can write

A x = B
L U x = B  
U x = y, L y=B.

We can solve L y =B by simple backsubstitution, and then solve U x = y again by back substitution. All the difficulties of solving linear systems of equations is moved into the problem of factoring the matrix A.

For cubic splines, the matrix A is a tridiagonal matrix (meaning, it has entries only along the diagonal, and above and below the diagonal).

The process of LU decomposition attempts to turn a matrix into the product of a lower triangular matrix and an upper triangular matrix. The procedure is outlined here.

Crout's Method (One’s on Diagonal of Upper Triangular Matrix) Applied to a Tridiagonal Matrix

Here is Crout's Method to LU decompose a matrix. We will use n=6 just to recall the structure of the process.

Clear[i,j,n]
n=6;

Here we define a general tridiagonal matrix and see what we can learn.

A=Normal[SparseArray[{i_,j_}/;Abs[i-j]≤1->a[i,j],{n,n}]];
MatrixForm[A]

The lower diagonal has zeros now:

L =Normal[SparseArray[{i_,j_}/;i≥j&&Abs[i-j]≤1->l[i,j],{n,n}]];
MatrixForm[L]

Here we have added the additional n=6 conditions that the lower diagonal matrix has 1's on the diagonal.

U =Normal[SparseArray[{i_,j_}/;i≤j&&Abs[i-j]≤1->u[i,j],{n,n}]];
Do[U[[i,i]]=1,{i,1,n}]
MatrixForm[U]

So now we have only 3n-2 undertermined quantities (2n-1 from lower, n-1 from upper). It is useful to look at the equations that have to be solved to create the LU factorization of a tridiagonal matrix, since many are automatically satisfied.

LU=L.U;

eqlist={};
varlist={};
Do[eqlist=Append[eqlist,LU[[i,j]]==A[[i,j]]],{i,1,n},{j,1,n}]
Do[varlist=Append[varlist,L[[i,i]]],{i,1,n}]
Do[varlist=Append[varlist,L[[i,i-1]]],{i,2,n}]
Do[varlist=Append[varlist,U[[i,i+1]]],{i,1,n-1}]

TableForm[eqlist]
varlist
Length[eqlist]
Length[varlist]

So we have 3n-2 equations in 3n-2 unknowns. As you can see, many of the equations are automatically satisfied. The seeds of a cubic spline algorithm are contained in the list of equations that have to be solved via Crout’s algorithm to determine the LU factorization.

Notice that the equations can be solved as follows to create the LU form (all other elements are zero):

First,
l[1,1]=a[1,1]
u[1,2]=a[1,2]/l[1,1]

Then, for i=2,3,...,n-1,
l[i,i-1]=a[i,i-1] 
l[i,i]=a[i,i]-l[i,i-1]u[i-1,i]
u[i,i+1]=a[i,i+1]/l[i,i]

Finally,
l[n,n-1]=a[n,n-1]
l[n,n]=a[n,n]-l[n,n-1]u[n-1,n]

Now, we have the LU decomposition for the tridiagonal matrix A. We use this as follows to solve for the vector c:

A c = B
Write as L U c = B  (use above recursion relations, so (Subscript[l, i,j]) and (Subscript[u, i,j]) are known)
Solve L y=B for y (use backsubstitution),
	(Subscript[y, 1])=((Subscript[B, 1])/(Subscript[l, 1,1]))
	(Subscript[y, i])=(((Subscript[B, i])-(Subscript[y, i-1])(Subscript[l, i,i-1]))/(Subscript[l, i,i])) for i=2,...,n
Solve U c = y for c (use backsubstitution),
	(Subscript[c, n])=(Subscript[y, n])
	(Subscript[c, i])=(Subscript[y, i])-(Subscript[u, i,i+1])(Subscript[c, i+1]) for i=n-1,1

For the cubic splines algorithm, this gives us the parameters c.

Natural Cubic Spline: Data Set with Algorithm

The Algorithm

To determine the cubic spline for n data points

For i=1 to n-1,
	(Subscript[a, i])=(Subscript[y, i]) 
	(Subscript[h, i])=(Subscript[x, i+1])-(Subscript[x, i])
	(Subscript[Δ, i])=(Subscript[y, i+1])-(Subscript[y, i])
 
LU Decomposition,
	(Subscript[L, 1,1])=1
	(Subscript[U, 1,2])=0
For i=2 to n-1,
	(Subscript[L, i,i-1])=(Subscript[h, i-1])
	(Subscript[L, i,i])=2(Subscript[h, i-1])+2(Subscript[h, i])-(Subscript[L, i,i-1])(Subscript[U, i-1,i]) 
	(Subscript[U, i,i+1])=((Subscript[h, i])/(Subscript[L, i,i])) 
then,
	(Subscript[L, n,n-1])=0
	(Subscript[L, n,n])=1-(Subscript[L, n,n-1])(Subscript[U, n-1,n])

Solve L y=B for y,
	(Subscript[y, 1])=0
For i=2 to n-1,
	(Subscript[y, i])=((3(((Subscript[Δ, i])/(Subscript[h, i]))-((Subscript[Δ, i-1])/(Subscript[h, i-1])))-(Subscript[y, i-1])(Subscript[L, i,i-1]))/(Subscript[L, i,i])) 
then,
	(Subscript[y, n])=((-(Subscript[y, n-1])(Subscript[L, n,n-1]))/(Subscript[L, n,n])) 

Solve U c = y for c,
	(Subscript[c, n])=(Subscript[y, n])
For i=n-1 to 1 by -1,
	(Subscript[c, i])=(Subscript[y, i])+(Subscript[U, i,i+1])(Subscript[c, i+1])

Define the other cubic spline parameters,
For i = 1 to n-1,
	(Subscript[b, i])=((Subscript[Δ, i])/(Subscript[h, i]))-((Subscript[h, i])/3)((Subscript[c, i+1])+2(Subscript[c, i])) 
	(Subscript[d, i])=(((Subscript[c, i+1])-(Subscript[c, i]))/(3(Subscript[h, i])))
	(Subscript[S, i])(x) = (Subscript[a, i])+(Subscript[b, i])(x-(Subscript[x, i]))+(Subscript[c, i])(((x-(Subscript[x, i])))^2)+(Subscript[d, i])(((x-(Subscript[x, i])))^3)

Define cubic spline by stitching together the result in each interval.

Applying the Algorithm

Here is a data set to try this out on:

Clear[data,L,U]
data={{0.9,1.3},{1.3,1.5},{1.9,1.85},{2.1,2.1},{2.6,2.6},{3.0,2.7},{3.9,2.4},{4.4,2.15}};

dataplot=ListPlot[data,PlotStyle->{Red,PointSize[0.02]},Frame->True,ImageSize->{500,Automatic}]

Clear[a,b,c,d,SS,xp,yp,L,U]
n=Length[data];

(*Define data points*)
Do[
xp[i]=data[[i,1]];
yp[i]=data[[i,2]];
,{i,1,n}]

(*Define a, h, delta*)
Do[
a[i]=yp[i];
h[i]=xp[i+1]-xp[i];
delta[i]=yp[i+1]-yp[i];
,{i,1,n-1}];

(*Perform the LU decomposition*)
L[1,1]=1;
U[1,2]=0;
Do[
L[i,i-1]=h[i-1];
L[i,i]=2h[i-1]+2h[i]-L[i,i-1]U[i-1,i];
U[i,i+1]=((h[i])/(L[i,i]));
,{i,2,n-1}];
L[n,n-1]=0;
L[n,n]=1-L[n,n-1]U[n-1,n];

(*Solve Ly=B for y*)
y[1]=0;
Do[
y[i]=((3(((delta[i])/(h[i]))-((delta[i-1])/(h[i-1])))-y[i-1]L[i,i-1])/(L[i,i]));
,{i,2,n-1}];
y[n]=((-y[n-1]L[n,n-1])/(L[n,n]));

(*Solve Uc=y for c*)
c[n]=y[n];
Do[
c[i]=y[i]-U[i,i+1]c[i+1];
,{i,n-1,1,-1}];

(*Define b and d*)
Do[
b[i]=((delta[i])/(h[i]))-((h[i])/3)(2c[i]+c[i+1]);
d[i]=((c[i+1]-c[i])/(3h[i]));
,{i,1,n-1}];

(*Create the cubic spline*)
SS[i_,x_]:=a[i]+b[i](x-xp[i])+c[i](((x-xp[i]))^2)+d[i](((x-xp[i]))^3);
S[x_]:=Do[ If[xp[i]≤x≤xp[i+1],Return[SS[i,x]]],{i,1,n-1}]

(*Output the cubic spline data*)
cubicsplinedata=Table[{i,xp[i],a[i],b[i],c[i],d[i]},{i,1,n-1}];
TableForm[SetPrecision[cubicsplinedata,4],TableSpacing->1,TableHeadings->{None,{"i","(Subscript[x, i])","(Subscript[a, i])","(Subscript[b, i])","(Subscript[c, i])","(Subscript[d, i])"}}]

Note that I defined the S[x] with equality on both sides of the interval--this is OK since the values are the same at the points, and allows us to use S[x] to get the values at both the left and right endpoints of the interval.

Here is a plot

csplot=Plot[S[x],{x,xp[1],xp[n]},PlotRange->{{0,6},{1,3}},FrameLabel->{"x","y"},Frame->True,Axes->False,PlotStyle->{{Thick,Blue}}];

Show[csplot,dataplot,ImageSize->{500,Automatic}]

If we want the derivative, the way things are defined computationally we can't just ask for the derivative. It's simple enough to get the derivative, however:

SSp[i_,x_]:=b[i]+2*c[i](x-xp[i])+3*d[i](((x-xp[i]))^2)
Sp[x_]:=Do[ If[xp[i]≤x≤xp[i+1],Return[SSp[i,x]]],{i,1,n-1}]

csdplot=Plot[Sp[x],{x,xp[1],xp[n]},FrameLabel->{"x","First Derivative"},Frame->True,Axes->False,PlotStyle->{Thick,Blue},ImageSize->{500,Automatic}]

Let's look at the second derivative. It should be continuous but not smooth.

SSpp[i_,x_]:=2*c[i]+6*d[i](x-xp[i])
Spp[x_]:=Do[ If[xp[i]≤x≤xp[i+1],Return[SSpp[i,x]]],{i,1,n-1}]

csddplot=Plot[Spp[x],{x,xp[1],xp[n]},FrameLabel->{"x","Second Derivative"},
Frame->True,Axes->False,PlotStyle->{Thick,Blue},ImageSize->{500,Automatic}]

Notice how the second derivative at the endpoints of the region is zero--this tells us this is a natural cubic spline.

Higher derivatives will be discontinuous.

SSppp[i_,x_]:=6*d[i]
Sppp[x_]:=Do[ If[xp[i]≤x≤xp[i+1],Return[SSppp[i,x]]],{i,1,n-1}]

csdddplot=Plot[Sppp[x],{x,xp[1],xp[n]},FrameLabel->{"x","Third Derivative"},Frame->True,Axes->False,PlotStyle->{Thick,Blue}, ImageSize->{500,Automatic}]

Alternate Algorithm You May See

Here is the algorithm to compute the cubic splines, rewritten in a slighty different form.

Clear[n,x,a,h,l,mu,z,b,c,d,i]

n=Length[data];

Do[
xp[i]=data[[i,1]];
a[i]=data[[i,2]],
{i,1,n}]

For[i=1,i≤n-1,i=i+1,h[i]=xp[i+1]-xp[i]];

For[i=2,i≤n-1,i=i+1,alpha[i]=(3/(h[i]))*(a[i+1]-a[i])-(3/(h[i-1]))*(a[i]-a[i-1])];

l[1]=1;
mu[1]=0;
z[1]=0;

For[i=2,i≤n-1,i=i+1,
{l[i]=2*(xp[i+1]-xp[i-1])-h[i-1]mu[i-1];
mu[i]=((h[i])/(l[i]));
z[i]=((alpha[i]-h[i-1]z[i-1])/(l[i]))}];

l[n]=1;
z[n]=0;
c[n]=0;

For[j=n-1,j≥ 1,j=j-1,
{c[j]=z[j]-mu[j]c[j+1];
b[j]=((a[j+1]-a[j])/(h[j]))-h[j]*((c[j+1]+2c[j])/3);
d[j]=((c[j+1]-c[j])/(3h[j]))}];

cubicsplinedata=Table[{j,xp[j],a[j],b[j],c[j],d[j]},{j,1,n-1}];

TableForm[SetPrecision[cubicsplinedata,4],TableSpacing->1,TableHeadings->{None,{"i","(Subscript[x, i])","(Subscript[a, i])","(Subscript[b, i])","(Subscript[c, i])","(Subscript[d, i])"}}]

Comparison: Newton's Interpolatory Divided Difference Formula and Cubic Splines

Let's compare the cubic spline method with Newton's Interpolating Polynomial.

Newton's Algorithm

Here is Newton's interpolating polynomial, defined through the divided difference representation:

Clear[f,divdif,len,list,prod,LP]

len=Length[data]-1;

Do[divdif[{Part[data,n,1]}]=Part[data,n,2],{n,1,len+1}];

divdif[{p_,q___,r_}]:=divdif[{p,q,r}]=((divdif[{q,r}]-divdif[{p,q}])/(r-p));

list[0]={Part[data,1,1]};
Do[list[k]=Append[list[k-1],Part[data,k+1,1]],{k,1,len}];
prod[0]=(x-Part[data,1,1]);
Do[prod[k]=prod[k-1]*(x-Part[data,k+1,1]),{k,1,len}];

This is Equation 3.2 in the text:

LP[x_]= divdif[list[0]]+(∑_(k=1)^len)divdif[list[k]]*prod[k-1];

lpplot=Plot[LP[x],{x,0,15},PlotRange->{{0,6},{1,3}},FrameLabel->{"x","y"},
Frame->True,Axes->False,PlotStyle->{Red,Thick}];

Show[lpplot,dataplot,ImageSize->{500,Automatic}]

Here are the derivatives. Since Newton's interpolating polynomial is a polynomial, getting derivatives is easy. I have used Mathematica to do them for simplicity of presentation, but they can be found by hand easily if needed from:

Expand[LP[x],x]

lpdplot=Plot[LP'[x],{x,0,6},FrameLabel->{"x","First Derivative"},
Frame->True,Axes->False,PlotStyle->{Thick,Red},ImageSize->{500,Automatic}]
lpddplot=Plot[LP''[x],{x,0,6},FrameLabel->{"x","Second Derivative"},Frame->True,Axes->False,PlotStyle->{Thick,Red},ImageSize->{500,Automatic}]
lpdddplot=Plot[LP'''[x],{x,0,6},FrameLabel->{"x","Third Derivative"},Frame->True,Axes->False,PlotStyle->{Thick,Red},ImageSize->{500,Automatic}]

Comparison Plots (Blue: cubic spline; Red: Lagrange Interpolating Polynomial)

Show[lpplot,csplot,dataplot,PlotRange->{{0,6},{1,3}},ImageSize->{500,Automatic}]

Notice the appearance of a small wiggle in the right for the Lagrange interpolating polynomial! Cubic spline avoids the wiggles.

Show[lpdplot,csdplot,PlotRange->{{0,6},{-1,2}},ImageSize->{500,Automatic}]

Show[lpddplot,csddplot,PlotRange->{{0,6},{-2.5,3.5}},ImageSize->{500,Automatic}]

Show[lpdddplot,csdddplot,PlotRange->{{0,6},{-20,10}},ImageSize->{500,Automatic}]

Remember, we have no derivative information for the function. One of these is not necessarily better than the other. They will get worse as the derivatives increase (third derivative estimate is going to be worse than the first derivative estimate).

References:

Heikki Ruskeepaa, Mathematica Navigator, Graphics and Methods of Applied Mathematics (1999) San Diego,Academic Press.
Richard Burden & J. Douglas Faires, Numerical Analysis, 7th ed. (2001) Pacific Grove, Brooks/Cole.