August 11, 1996  test-2d-ipl-code.nw  1

1 Tests of 2d interpolation for L_{α} functions

1.1 The Code

To see how accurately the unbounded derivatives of L_{α} functions can be approximated in 2D, we use the following maple code.

1a ⟨test-2d-ipl.map 1a⟩≡

read `sinc8.mas`;       # some standard sinc subroutines
Digits:=14;

⟨procedure for relative error between two matrices 4b⟩

⟨define function and derivatives 1b⟩ ⟨conformal maps and interpolant 2a⟩ ⟨set global parameters 2b⟩ ⟨produce interpolant’s coefficient grid 2c⟩

⟨procedure: get interpolant value 2d⟩ ⟨procedure: get interpolant’s dx value 3a⟩ ⟨procedure: get interpolant’s dy value 3b⟩

⟨get interpolant and function values on a grid 5⟩ ⟨write values, interpolant, and error into file 3c⟩

⟨get derivative interpolant and derivative values on a grid 6⟩ ⟨append derivatives, interpolant, and error to file 4a⟩

1b ⟨define function and derivatives 1b⟩≡

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

2a ⟨conformal maps and interpolant 2a⟩≡

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

2b ⟨set global parameters 2b⟩≡

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.

2c ⟨produce interpolant’s coefficient grid 2c⟩≡

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:

2d ⟨procedure: get interpolant value 2d⟩≡

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;

3a ⟨procedure: get interpolant’s dx value 3a⟩≡

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;

3b ⟨procedure: get interpolant’s dy value 3b⟩≡

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 written to the output file, but derivatives are subsequently appended to the same file.

3c ⟨write values, interpolant, and error into file 3c⟩≡

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

4a ⟨append derivatives, interpolant, and error to file 4a⟩≡

#       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.

4b ⟨procedure for relative error between two matrices 4b⟩≡

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;

5 ⟨get interpolant and function values on a grid 5⟩≡

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;

6 ⟨get derivative interpolant and derivative values on a grid 6⟩≡

#       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.

August 11, 1996  test-2d-ipl-test_runs.nw  7

1.2 Some test runs and results

For the modifications

f(x, y) = \frac{\sqrt{x} (1 − x) y^{1/10} (1 − y)^{2}}{3}  (1)

and N = 10 made to the listed code, we have the following relative errors: For f,

[ \begin{matrix}0.03479 & 0.02597 & 0.02024 & 0.02052 & 0.01858 \\ −0.007716 & −0.003380 & −0.007902 & −0.003937 & −0.006782 \\ −0.01327 & −0.008304 & −0.01731 & −0.02282 & −0.02138 \\ 0.01171 & 0.004684 & 0.005910 & 0.0001523 & −0.0004227 \\ 0.001112 & −0.009258 & −0.007214 & −0.01015 & −0.01108 \end{matrix} ]

for ∂f/∂x:

[ \begin{matrix}0.02671 & −0.008716 & 0.02433 & 0.009778 & −0.01909 \\ −0.03521 & −0.07032 & 0.03061 & −0.006469 & −0.07817 \\ 0.004837 & −0.04405 & −0.06126 & −0.04526 & 0.05540 \\ 0.02411 & 0.001753 & −0.01669 & −0.04427 & 0.04858 \\ −0.003292 & −0.01564 & −0.01042 & −0.04122 & 0.01453 \end{matrix} ]

for ∂f/∂y:

[ \begin{matrix}0.07000 & 0.2547 & 0.07463 & −0.2004 & −0.1303 \\ 4.226 & 2.744 & 4.528 & 4.940 & 4.182 \\ −0.6513 & −0.4244 & −0.6238 & −0.5545 & −0.5229 \\ −0.03197 & 0.05759 & −0.01336 & −0.04826 & −0.03390 \\ 0.1819 & 0.1774 & 0.2088 & 0.1813 & 0.1747 \end{matrix} ]

on the grid with x-values

[ \begin{matrix}0.01 & 0.05750 & 0.1050 & 0.1525 & 0.2000 \\ 0.01 & 0.05750 & 0.1050 & 0.1525 & 0.2000 \\ 0.01 & 0.05750 & 0.1050 & 0.1525 & 0.2000 \\ 0.01 & 0.05750 & 0.1050 & 0.1525 & 0.2000 \\ 0.01 & 0.05750 & 0.1050 & 0.1525 & 0.2000 \end{matrix} ]

August 11, 1996  test-2d-ipl-test_runs.nw  8

and y-values

[ \begin{matrix}0.01 & 0.01 & 0.01 & 0.01 & 0.01 \\ 0.05750 & 0.05750 & 0.05750 & 0.05750 & 0.05750 \\ 0.1050 & 0.1050 & 0.1050 & 0.1050 & 0.1050 \\ 0.1525 & 0.1525 & 0.1525 & 0.1525 & 0.1525 \\ 0.2000 & 0.2000 & 0.2000 & 0.2000 & 0.2000 \end{matrix} ]

With the “nicer” function

f(x, y) = \frac{\sqrt{x} (1 − x) y^{3/5} (1 − y)^{2}}{3}  (2)

and N = 10, using the same (x, y) grid, we have the relative errors

[ \begin{matrix}0.0147 & 0.00644 & 0 & 0.00147 & 0 \\ −0.00757 & 0 & −0.00646 & 0 & −0.00524 \\ −0.0102 & −0.00641 & −0.0150 & −0.0175 & −0.0202 \\ −0.00261 & −0.0114 & −0.00445 & −0.0156 & −0.0144 \\ 0.00249 & −0.00545 & −0.00425 & −0.00744 & −0.0103 \end{matrix} ]

[ \begin{matrix}−0.0100 & −0.0281 & 0.00918 & 0 & −0.0380 \\ −0.0309 & −0.0695 & 0.0319 & −0.00540 & −0.0754 \\ 0.00896 & −0.0420 & −0.0589 & −0.0438 & 0.0583 \\ 0.00532 & −0.00748 & −0.0293 & −0.0520 & 0.0375 \\ −0.00254 & −0.0143 & −0.00932 & −0.0408 & 0.0138 \end{matrix} ]

and

[ \begin{matrix}0 & 0.0259 & −0.00289 & −0.0379 & −0.0304 \\ −0.0569 & 0.00499 & −0.0700 & −0.0954 & −0.0693 \\ 0.0378 & −0.0295 & 0.0273 & 0 & −0.00349 \\ 0.0994 & −0.0181 & 0.0622 & 0.0867 & 0.0665 \\ 0.102 & 0.0185 & −0.0764 & −0.00744 & −0.00688 \end{matrix} ]

for f, ∂f/∂x and ∂f/∂y, respectively.

With the non-separable function

f(x, y) = ^{4} \sqrt{x^{2} + y^{2}}  (3)

August 11, 1996  test-2d-ipl-test_runs.nw  9

and N = 10, again on the same grid, we get the following results. Notice that this function is not in L_{α} on the unit square, so poor results can be expected.

[ \begin{matrix}−0.0841 & 0 & 0.00616 & −0.0153 & 0.00894 \\ −0.0124 & 0 & −0.00289 & −0.0198 & −0.00438 \\ −0.00616 & −0.00289 & −0.0104 & −0.0372 & −0.0168 \\ −0.0179 & −0.0248 & −0.0209 & −0.0345 & −0.0179 \\ 0.0134 & −0.00658 & −0.00421 & −0.0160 & −0.00752 \end{matrix} ]  (4)

[ \begin{matrix}−0.0875 & 0.324 & −0.235 & 0.0705 & 0.254 \\ −1.80 & 0.231 & −0.205 & 0.0173 & 0.161 \\ −0.103 & 0.388 & −0.396 & −0.0199 & 0.271 \\ 2.29 & 0.538 & −0.317 & −0.00657 & 0.274 \\ −4.27 & 0.367 & −0.362 & −0.0546 & 0.265 \end{matrix} ]  (5)

[ \begin{matrix}−0.138 & −0.110 & −2.13 & −13.6 & −16.1 \\ 0.289 & 0.377 & 0.147 & −0.357 & −0.499 \\ −0.215 & −0.268 & −0.210 & −0.162 & −0.200 \\ 0.118 & 0.0173 & 0.0606 & 0.130 & 0.101 \\ 0.302 & 0.207 & 0.170 & 0.179 & 0.158 \end{matrix} ]  (6)

for f, ∂f/∂x and ∂f/∂y, respectively.

Next, an L_{α} function with a radial singularity at (0, 0).

f(x, y) = −y^{3/2} (−1 + y) (x − 1) + ^{4} \sqrt{x^{2} + y^{2}} (1 − x) (1 − y) y  (7)

With N = 10 and the above grid, we have the relative errors of

[ \begin{matrix}0.0162 & 0 & −0.00502 & −0.00410 & −0.00363 \\ 0.0613 & 0 & −0.00777 & −0.00398 & −0.00640 \\ −0.189 & −0.0103 & −0.0136 & −0.0236 & −0.0176 \\ −0.339 & −0.00622 & −0.00869 & −0.0210 & −0.00871 \\ −0.271 & −0.00739 & −0.00498 & −0.0149 & −0.00923 \end{matrix} ]  (8)

[ \begin{matrix}0.00691 & −0.0170 & −0.00881 & 0.0128 & −0.00371 \\ 0.00528 & −0.0246 & 0.00359 & −0.00451 & −0.0265 \\ −0.0814 & 0.0151 & −0.0434 & −0.0257 & 0.0143 \\ −0.141 & 0.0427 & −0.0436 & −0.0203 & 0.0177 \\ −0.122 & 0.0339 & −0.0305 & −0.0136 & 0.0126 \end{matrix} ]  (9)

[ \begin{matrix}0.00662 & 0.0194 & −0.0130 & −0.0349 & −0.0257 \\ 0.190 & −0.0975 & −0.108 & −0.0931 & −0.0527 \\ 0.785 & 0.0236 & −0.0421 & −1.31 & −0.0367 \\ −0.214 & −0.00650 & −0.0553 & −0.120 & −0.211 \\ −1.09 & −0.0171 & −0.0172 & −0.0546 & −0.0402 \end{matrix} ]  (10)

for f, ∂f/∂x and ∂f/∂y, respectively. The (4,1) and (5,1) entries in (9) correspond to small exact values; similarly, the poor entries of (10) also correspond to small exact answers.

As yet another nice test, consider

f(x, y) = \sqrt{x} ^{3} \sqrt{1 − x}y^{1/10} (1 − y)^{2} sin(2 π x)e^{y}  (11)

with N = 10 and the rest as coded. For the function interpolation, we have the relative errors

[ \begin{matrix}0.0385 & 0.0251 & 0.0168 & 0.0212 & 0.0203 \\ −0.00452 & −0.00680 & −0.00738 & −0.00935 & −0.00358 \\ −0.0360 & −0.00851 & −0.0147 & −0.0233 & −0.0178 \\ −0.0115 & 0.00174 & 0 & −0.00476 & 0 \\ −0.0240 & −0.0127 & −0.00784 & −0.00993 & −0.0114 \end{matrix} ]  (12)

Along the diagonal from (0,0) to (1,1), the relative error is similar, with a maximum of 5%.

Without the random error, we obtain

[ \begin{matrix}0.0667 & 0.0289 & 0.0251 & 0.0317 & 0.0324 \\ 0.0362 & −0.00341 & −0.00738 & 0 & 0 \\ 0.0360 & −0.00340 & 0 & 0 & 0 \\ 0.0508 & 0.0104 & 0.00753 & 0.0143 & 0.0146 \\ 0.0337 & −0.00543 & −0.00784 & 0 & 0 \end{matrix} ]  (13)

The results along the diagonal are again comparable, with a maximum relative error of 5%. By increasing N to 13, we obtain the relative error of

[ \begin{matrix}0.00770 & 0.0193 & 0.00836 & 0.0212 & 0.0203 \\ −0.0203 & −0.0119 & −0.0148 & −0.0141 & −0.00716 \\ 0 & 0.0119 & 0.00735 & 0.00465 & 0.0107 \\ −0.0161 & −0.00522 & −0.0151 & −0.00952 & −0.00364 \\ −0.0168 & −0.00907 & −0.00784 & −0.00497 & −0.00380 \end{matrix} ]  (14)

Along the diagonal, the maximum relative error is still 5%.

If we define ρ by

ρ = r^{\frac{1}{2}} sin(θ) = \frac{y}{(x^{2} + y^{2})^{1/4}}  (15)

with r the radial distance from the origin, and θ measured from the x-axis as usual, then the function

f(x, y) = ρ(1 − x)(1 − y) − \sqrt{y}(1 − x)(1 − y)  (16)

is in L_{α} on the unit square. With N = 10 and the above grid, we have the relative errors of

[ \begin{matrix}0.0449 & 0.0110 & 0.00163 & −0.00160 & 0 \\ −0.00602 & −0.00295 & −0.00644 & −0.00514 & −0.00583 \\ −0.196 & −0.0115 & −0.0145 & −0.0231 & −0.0203 \\ 0.723 & 0.0216 & 0.00365 & −0.0134 & −0.0119 \\ −0.357 & −0.0106 & −0.00529 & −0.0152 & −0.0110 \end{matrix} ]  (17)

[ \begin{matrix}0 & 0 & −0.138 & −0.0471 & −0.0336 \\ −0.0185 & −0.0148 & −0.00225 & 0 & −0.0356 \\ −0.0546 & 0 & −0.0380 & −0.0265 & 0.0276 \\ −0.0129 & 0.0317 & −0.0394 & −0.0269 & 0.0198 \\ −0.118 & 0.0282 & −0.0277 & −0.0139 & 0.0164 \end{matrix} ]  (18)

[ \begin{matrix}−0.0231 & 0.107 & −0.00605 & −0.0703 & −0.0516 \\ 0.784 & 0.0404 & 0.108 & 0.228 & 0.443 \\ −1.64 & −0.0688 & −0.0650 & −0.0425 & −0.0436 \\ 0.576 & 0.0197 & −0.0133 & −0.0341 & −0.0394 \\ 4.28 & 0.140 & 0.0575 & 0.00482 & 0 \end{matrix} ]  (19)

for f, ∂f/∂x and ∂f/∂y, respectively. The 0.443 entry in the y-partial corresponds to a small entry.

For the same f and N, but without introducing the random error, the relative errors of

[ \begin{matrix}0.0449 & 0.0146 & 0.0114 & 0.0112 & 0.0114 \\ −0.0181 & 0 & 0 & 0 & 0 \\ −0.0911 & −0.00576 & −0.00242 & 0 & −0.00135 \\ 0.974 & 0.0402 & 0.0146 & 0.00672 & 0.00683 \\ −0.163 & −0.00453 & 0 & −0.00305 & −0.00220 \end{matrix} ]  (20)

[ \begin{matrix}0.00489 & 0 & −0.0199 & 0.0637 & 0.0122 \\ −0.00308 & 0.00269 & −0.00225 & 0 & 0 \\ −0.00781 & −0.00404 & 0.00422 & 0 & 0 \\ 0.0703 & 0.00633 & −0.00262 & 0 & 0 \\ −0.0273 & 0.00471 & −0.00346 & 0 & 0.00410 \end{matrix} ]  (21)

[ \begin{matrix}−0.00771 & 0.0273 & 0.0121 & 0.00502 & 0.00939 \\ 0.809 & 0.0827 & 0.0841 & 0.136 & 0.352 \\ −1.81 & −0.101 & −0.0542 & −0.0504 & −0.0591 \\ 0.734 & 0.0344 & 0.0133 & 0.0137 & 0.00657 \\ 4.79 & 0.199 & 0.0719 & 0.0434 & 0.0410 \end{matrix} ]  (22)

are obtained for f, ∂f/∂x and ∂f/∂y, respectively. Notice the large errors for small x, in the first columns of the matrices. Along the diagonal, for x = y = 0.1 : 0.1 : 0.9, the largest relative error is 0.5%. For x = 0.02, y = 0.1 : 0.1 : 0.9 (near the left boundary), we get

[0.065, 0.024, 0.511, −1.32, 0.00302, 3.32, −3.32, −0.815, −8.07]

for the values of the relative error. For x = 0.1 : 0.1 : 0.9, y = 0.98 (the top boundary), the r.e. is

[1.29, 11.7, −1.13, −0.473, −0.267, −0.179, −0.137, −0.107, −0.0893]

while for y = 0.1 : 0.1 : 0.9, x = 0.98 (the right boundary) it is much improved:

[0, 0, 0.00283, 0, 0, 0.00579, 0, −0.00137, −0.00598]

For x = 0.1 : 0.1 : 0.9, y = 0.02 (the bottom boundary), we also have the much better error

[0.00896, 0.00779, 0.00723, 0.00665, 0.00658, 0.00659, 0.00633, 0.00634, 0.00622]

For the above f, N, without random error and with h = \frac{π}{\sqrt{N}} we have the following errors near the boundaries: On the top,

[−0.0108, −0.151, −0.0946, −0.0339, 0.000233, 0.0123, 0.00746, −0.00496, 0.00899]

(23)

on the bottom,

[0.00197, −0.000588, 0.000203, 0.000743, 0.000122, −0.000729, −0.000223, 0.00203, −0.00152]

(24)

on the left,

[−0.0453, 0.170, −0.0583, −0.39, −0.0761, 0.828, 0.556, −2.692, 3.44] (25)

and on the right

[0.0102, −0.00506, −0.0203, −0.00692, −0.00441, −0.0237, −0.0172, −0.0321, −0.0210]

(26)

The graphs comparing the exact answer and the interpolant are shown in the figures.