3D Spline Interpolation 

Andreas Schramm 

Reinhold-Wuerth-University, Kuenzelsau, Germany

andreas.schramm@hs-heilbronn.de

Description:  

This worksheet gives an example of how to use the Maple spline function to create a 3 dimensional spline surface and a function R^2->R, based on discrete values given with respect to their axes in a matrix with the first row (without the first element) containing the x-values and the first column (without the first element) containing the y-values.  

The computet function and surface is continuous and differentiable. The relevant difference between least square method and spline interpolation is that no knowledge on the mathematical function describing the system behavior of the source which delivered the measured or computed values is needed. The example below results from a computation of a finite element method calculation. The next step is to get interpolated values for the purpose of minimizing calculation time of the finite element method calculation. 

Maple Version: 11.02 

Loading of toolboxes and defining matrix 

>	restart;with(plots):with(ListTools):with(linalg):

 

Warning, the name changecoords has been redefined

 

Warning, the protected names norm and trace have been redefined and unprotected

 

The above described test matrix is defined here. 

>

Value_matr:=matrix([[0., 51.750, 46., 40.250, 34.500, 28.750, 23., 17.250, 11.500, 5.750, 0.], [0., 15.223, 14.300, 15.930, 21.222, 29.982, 40.086, 49.454, 57.060, 62.316, 64.653], [-5.750, 15.223, 14.300, 15.930, 21.222, 29.982, 40.086, 49.454, 57.060, 62.316, 64.653], [-11.500, 32.264, 33.044, 35.314, 39.346, 44.597, 49.939, 54.684, 58.461, 61.035, 62.169], [-17.250, 37.816, 39.116, 41.455, 44.771, 48.497, 52.060, 55.179, 57.638, 59.303, 60.023], [-23., 40.480, 41.979, 44.206, 46.945, 49.755, 52.368, 54.638, 56.437, 57.644, 58.157], [-28.750, 41.961, 43.511, 45.544, 47.807, 50.015, 52.042, 53.794, 55.188, 56.122, 56.513], [-34.500, 42.830, 44.344, 46.154, 48.037, 49.831, 51.462, 52.865, 53.978, 54.728, 55.034], [-40.250, 43.327, 44.755, 46.349, 47.941, 49.433, 50.779, 51.927, 52.835, 53.440, 53.684], [-46., 43.580, 44.900, 46.299, 47.662, 48.930, 50.064, 51.020, 51.763, 52.242, 52.437], [-51.750, 43.658, 44.868, 46.101, 47.285, 48.378, 49.346, 50.153, 50.761, 51.133, 51.280]]):evalm(Value_matr);

 

(1.1)

 

The x values of the contained data is stored in the first row, without the first element. Note that in this case the x-values are in descending order. 

>	row(Value_matr,1);

 

(1.2)

 

The y values of the contained data is stored in the first column, without the first element. Note that in this case the y-values are in descending order. 

>	col(Value_matr,1);

 

(1.3)

 

Data processing 

For better accessibility of the data (for following plotting purposes) the data has to be processed in lists: 

First we have to extract the x-values and create a list which contains m times (m = number of columns of the matrix) the x-values.  

>	xlist:=[]:

 

>	for m from 1 to coldim(Value_matr)-1 do

 

>	for n from (rowdim(Value_matr)-(rowdim(Value_matr)-2)) to rowdim(Value_matr) do

 

>	xlist:=[xlist,Value_matr[1,n]]:

 

>

end do:

 

>

end do:

 

>	xlist:=Flatten(xlist);

 

(2.1)

 

We have to do nearly the same for the y-values: creating a list which contains m times (m = number of rows of the matrix) the y-values. Note that instead of combining the interval of the y-values m times, here the y-values themselfes have to be listed after each other. 

>	ylist:=[]:

 

>	for n from (coldim(Value_matr)-(coldim(Value_matr)-2)) to coldim(Value_matr) do

 

>	for m from 1 to rowdim(Value_matr)-1 do

 

>	ylist:=[ylist,Value_matr[n,1]]

 

>

end do:

 

>

end do:

 

>	ylist:=Flatten(ylist);

 

(2.2)

 

The values are also processed into a list: 

>	Valuelist:=[]:

 

>	for n from 2 to rowdim(Value_matr) do

 

>	for m from 2 to coldim(Value_matr) do

 

>	Valuelist:=[Valuelist,Value_matr[n,m]]:

 

>

end do:

 

>

end do:

 

>	Valuelist:=Flatten(Valuelist);

 

(2.3)

 

This operation determines the plotting intervals: 

>	y_max := max(seq(ylist[i],i=1..nops(ylist))); y_min := min(seq(ylist[i],i=1..nops(ylist)));

 

(2.4)

 

(2.4)

 

>	x_max := max(seq(xlist[i],i=1..nops(xlist))); x_min := min(seq(xlist[i],i=1..nops(xlist)));

 

(2.5)

 

(2.5)

 

After the data processing we can get to the core of the matter: 

Spline Interpolation 

For better accessibility we need to do some data processing: 

We extract the x-values eliminating the non-information containing zero of the first element in the first row to get the raw data of the x-values: 

>	Valuelist_Spline_x_temp:=[]:

 

>	for m from 2 to coldim(Value_matr) do

 

>	Valuelist_Spline_x_temp:=[Valuelist_Spline_x_temp,Value_matr[1,m]]:

 

>

end do:

 

>	Valuelist_Spline_x_temp:=Flatten(Valuelist_Spline_x_temp);

 

(3.1)

 

The first step is to build a set of curves calculated by spline interpolation with respect to variable x with parameter y: 

>	Value_matr_Spline_y:=Matrix(coldim(Value_matr)-1,2):

 

>	for n from 2 to rowdim(Value_matr) do

 

>	Value_matr_Spline_y[n-1,1]:=Value_matr[n,1]:

 

>	Valuelist_Spline:=[]:

 

>	for m from 2 to coldim(Value_matr) do

 

>	Valuelist_Spline:=[Valuelist_Spline,Value_matr[n,m]]:

 

>

end do:

 

>	Valuelist_Spline:=Flatten(Valuelist_Spline):

 

>	Value_matr_Spline_y[n-1,2]:=spline(Reverse(Valuelist_Spline_x_temp),Reverse(Valuelist_Spline),x):

 

>

end do:

 

>

 

Display of the intermittent results 

For better understanding we can plot the set of curves like following: 

>	plotlist:=[]:

 

>	for n from 1 to rowdim(Value_matr_Spline_y) do

 

>	plotlist:=[plotlist,plot3d([x,Value_matr_Spline_y[n,1],Value_matr_Spline_y[n,2]],x=x_min..x_max,y=y_min..y_max)]:

 

>

end do:

 

>	display(Flatten(plotlist),axes=boxed,thickness=2,style=hidden,title="Resulting set of curves after 1st step of spline interpolation");

 

 

The second step is to interpolate the set of curves in the second dimension, which is y: 

>	ylist_splines:=[]:Splinelist_splines:=[]:

 

>	for n from 1 to rowdim(Value_matr_Spline_y) do

 

>	ylist_splines:=[ylist_splines,Value_matr_Spline_y[n,1]]:

 

>	Splinelist_splines:=[Splinelist_splines,Value_matr_Spline_y[n,2]]:

 

>

end do:

 

>	ylist_splines:=Flatten(ylist_splines):Splinelist_splines:=Flatten(Splinelist_splines):

 

>	Value_matr_Spline:=spline(Reverse(ylist_splines),Reverse(Splinelist_splines),y):

 

>	Value_matr_Spline:=unapply(Value_matr_Spline,x,y):

 

Display of the resulting spline surface 

The blue circles point out the discrete data given in the above defined matrix. 

>	data_points:=pointplot3d( {seq([xlist[i],ylist[i],Valuelist[i]], i=1..nops(Valuelist)) }, symbol=circle, color=blue):

 

>	spline_surface:=plot3d(Value_matr_Spline(x,y), x=x_min..x_max, y=y_min..y_max, style=wireframe):

 

>	display3d([spline_surface,data_points], labels=['x','y','"f(x,y)"'], axes=normal, orientation=[145,45],title="Resulting spline-surface");

 

 

Conclusion 

Spline interpolation serves its purpose if you want to interpolate between data points, without knowledge of the real function. The method will give you good results in the intervals for x and y. But be aware of extrapolating out of the boundaries of x and y. Because of the method of spline interpolation the absolute return values of the interpolation function will increase to infinity. In most cases it will fail to represent real behavior.  

>	spline_surface:=plot3d(Value_matr_Spline(x,y), x=x_min..1.5*x_max, y=1.5*y_min..y_max, style=wireframe):

 

>	display3d([spline_surface,data_points], labels=['x','y','"f(x,y)"'], axes=normal, orientation=[145,45],title="Resulting spline-surface with expanded x and y intervals");

 

 

>

 

>

 

Legal Notice: The copyright for this application is owned by the author(s). Neither Maplesoft nor the author are responsible for any errors contained within and are not liable for any damages resulting from the use of this material. This application is intended for non-commercial, non-profit use only. Contact the author for permission if you wish to use this application in for-profit activities.