\section{Tests of 2d interpolation for $L_{\alpha}$ functions}
\subsection{The Code}
To see how accurately the unbounded derivatives of $L_{\alpha}$
functions can be approximated in 2D, we use the following  {\sc
maple} code.

<<test-2d-ipl.map>>=
read `sinc8.mas`;	# some standard sinc subroutines
Digits:=14;

<<procedure for relative error between two matrices>>

<<define function and derivatives>>
<<conformal maps and interpolant>>
<<set global parameters>>
<<produce interpolant's coefficient grid>>

<<procedure:  get interpolant value>>
<<procedure:  get interpolant's dx value>>
<<procedure:  get interpolant's dy value>>

<<get interpolant and function values on a grid>>
<<write values, interpolant, and error into file>>

<<get derivative interpolant and derivative values on a grid>>
<<append derivatives, interpolant, and error to file>>
@


<<define function and derivatives>>=
f := unapply(x**(1/2)*(1-x)**1/3*y**(1/10)*(1-y)**(2), (x,y));
fx := unapply(diff(f(x,y),x), (x,y));
fy := unapply(diff(f(x,y),y), (x,y));
@

<<conformal maps and interpolant>>=
phi := unapply(ln(x/(1-x)), x);		# from [0,1] to R
z = phi(u):
solve(", u):
psi := unapply(", z);			# from R to [0,1]

# 	the interpolant
si := unapply(sinc((phi(x) - k*h)/h)*sinc((phi(y) - j*h)/h), (k,x,j,y));
#	define needed partials of si:
si_x := unapply(diff(si(l,x,k,y), x), (l,x,k,y));
si_y := unapply(diff(si(l,x,k,y), y), (l,x,k,y));
@

<<set global parameters>>=
N := 2;		# number of terms in series.
h := evalf(Pi/sqrt(2*N));
m := 2*N+1;
@

For interpolation, the coefficients of the sinc product are formed
from function evaluations.  To get an idea of stability
w.r.t. perturbations, we include a random error.

<<produce interpolant's coefficient grid>>=
mran := evalf(1+rand(1..1000000)/50000000);	# at most 2% error
with(linalg):
si_coe := matrix(m,m, 0):
for i from -N to N do;		# y
   for j from -N to N do;	# x
      si_coe[i+N+1,j+N+1] := f(psi(j*h), psi(i*h))*mran();
   od;	# i
od;	# j
@
The si\_coe matrix computed above is kept global, and used to compute
the interpolants values at non-sinc grid points:

<<procedure:  get interpolant value>>=
get_int_val := proc(x,y)
   local k,l,val;
   val := 0;
   for k from -N to N do;
      for l from -N to N do;
	 val := val + evalf(si_coe[k+N+1,l+1+N]*si(l,x, k,y))
      od;
   od;
   val;
end;
@

<<procedure:  get interpolant's dx value>>=
get_int_dx := proc(x,y)
   local k,l,val;
   val := 0;
   for k from -N to N do;
      for l from -N to N do;
	 val := val + evalf(si_coe[k+N+1,l+1+N]*si_x(l,x, k,y))
      od;
   od;
   val;
end;
@

<<procedure:  get interpolant's dy value>>=
get_int_dy := proc(x,y)
   local k,l,val;
   val := 0;
   for k from -N to N do;
      for l from -N to N do;
	 val := val + evalf(si_coe[k+N+1,l+1+N]*si_y(l,x, k,y))
      od;
   od;
   val;
end;
@


The values are {\sl written} to the output file, but derivatives are
subsequently {\sl appended} to the same file.

<<write values, interpolant, and error into file>>=
# 	for humans to read:
tmp := `/tmp/test-2d-ipl.out.`.N;
writeto(tmp);
lprint(`For N=`, N, ` we have the following: `);
lprint(`The exact function values:`);
op(func_grid);
lprint(`The interpolant values:`);
op(ipl_grid);
lprint(`The relative error:`);
rel_err(func_grid, ipl_grid);
writeto(terminal);
# 	for maple to re-read:
tmp := `/tmp/test-2d-ipl.out.`.N . `.m`;
save func_grid, ipl_grid, tmp;
@


<<append derivatives, interpolant, and error to file>>=
# 	for humans to read:
tmp := `/tmp/test-2d-ipl.out.`.N;
appendto(tmp);

lprint(`The exact dx values:`);
op(dx_func_grid);
lprint(`The interpolant dx values:`);
op(dx_ipl_grid);
lprint(`The dx relative error:`);
rel_err(dx_func_grid, dx_ipl_grid);

lprint(`The exact dy values:`);
op(dy_func_grid);
lprint(`The interpolant dy values:`);
op(dy_ipl_grid);
lprint(`The dy relative error:`);
rel_err(dy_func_grid, dy_ipl_grid);

appendto(terminal);
#
# 	for maple to re-read:
tmp := `/tmp/test-2d-ipl.derivatives.out.`.N . `.m`;
save dx_func_grid, dx_ipl_grid, dy_func_grid, dy_ipl_grid, tmp;
@


In rel\_err, the first argument is taken as the exact answer.

<<procedure for relative error between two matrices>>=
rel_err := proc(a,b)
   local i,j, mmat;
   mmat := matrix(rowdim(a), coldim(a), 0);
   for i from 1 to rowdim(a) do;
      for j from 1 to coldim(a) do;
	 mmat[i,j] := (a[i,j] - b[i,j])/a[i,j];
      od;
   od;
   op(mmat);
end;

@

<<get interpolant and function values on a grid>>=
x_from := 0.01;
x_to := 0.2;
x_step := (x_to - x_from)/4;
y_from := 0.01;
y_to := 0.2;
y_step := (y_to - y_from)/4;
x_grid := matrix(5,5,0);
y_grid := matrix(5,5,0);
ipl_grid := matrix(5,5,0);
func_grid := matrix(5,5,0);
j:= 0;		# integer indicies
for x from x_from by x_step to x_to do;
   i := 0;    
   j := j+1;
   for y from y_from by y_step to y_to do;
      i := i+1;
      func_grid[i,j] := f(x,y);
      x_grid[i,j] := x;
      y_grid[i,j] := y;
      #		get interpolant value:
      ipl_grid[i,j] := get_int_val(x,y);
   od;
od;

@

<<get derivative interpolant and derivative values on a grid>>=
# 	use the same points that were used in the function comparison
dx_ipl_grid := matrix(5,5,0);
dx_func_grid := matrix(5,5,0);
dy_ipl_grid := matrix(5,5,0);
dy_func_grid := matrix(5,5,0);
#
j:= 0;		# integer indicies
for x from x_from by x_step to x_to do;
   i := 0;    
   j := j+1;
   for y from y_from by y_step to y_to do;
      i := i+1;
      dx_func_grid[i,j] := fx(x,y);
      dy_func_grid[i,j] := fy(x,y);
      x_grid[i,j] := x;
      y_grid[i,j] := y;
      dx_ipl_grid[i,j] := get_int_dx(x,y);
      dy_ipl_grid[i,j] := get_int_dy(x,y);
   od;
od;
@

Notice that for simplicity, the domain is the unit square and $k$ and
$h$ are equal for both regions.