\section{Interpolation of unbounded derivatives of $L_{\alpha}$ functions}
The purpose of this paper is to test the quality of approximation
of $f(x)$ and $f'(x)$ on [0,1] using the sinc series
  \begin{equation}
    \sum _{k=-N}^{N} c_{k}S(k,h)\circ \phi(x_{j})
  \end{equation}
{\sl but with the $c_{k}$ obtained  by setting up and solving the
collocation system } 
  \begin{equation}
    \sum _{k=-N}^{N} c_{k} \left [ S(k,h)\circ \phi(x_{j}) \right ]' =
    f'(x_{j}) 
    \label{eqn:derivative-collocation} 
  \end{equation}
instead of using simple interpolation of $f$ at sinc points.  The
advantage, if sucessful, comes in the solution of differential
equations for which derivative boundary conditions are specified even
though the derivatives are known to be unbounded at certain points.
For boundaries with unbounded normal derivatives, one can then use the
series 
  \begin{equation}
    \sum _{k=-N}^{N} c_{k} \left [ S(k,h)\circ \phi(x_{j}) \right ]' 
  \end{equation}
rather than 
  \begin{equation}
    \sum _{k=-N}^{N} c_{k} \left [ S(k,h)\circ \phi(x_{j}) \right ]
  \end{equation}
in a sinc-spline product for approximation of the derivative.

\subsection{The Code}
The code will consist of 
  \begin{enumerate}
    \item a {\sc maple} file to generate needed C++ subroutines,
    \item the main C++ program to set up and solve the system,
    \item a second C++ program to plot the series and derivative.
  \end{enumerate}
To smoothly integrate the three and the {\sc maple} output, a
Makefile will be used\footnote{ To 
extract needed pieces from this noweb source, the makefile (notice the
lowercase m) is used.}.

Starting with the main setup loop implementing
  \begin{equation}
    \sum _{k=-N}^{N} c_{k} \left [ S(k,h)\circ \phi(x_{j}) \right ]' =
    f'(x_{j}) 
  \end{equation}
where the $x_{j}$ are the sinc points, $j=-M..N$, we have
<<c setup: main loop>>=
{
  // 	c setup: main loop
  int k,j;
  for (j=-M;  j<=N;  j++) {
    for (k=-M;  k<=N;  k++){
      A_PUT_ENTRY(j+M, k+M, skh_pri(k, h, xj(j) ));
    }
    B_PUT_ENTRY(j+M, f_pri( xj(j) ));
  }     // for j
}
@
%
This code segment requires some definitions and functions.  
 First, the matrix and rhs vector, and the access macros.
<<c setup: matrix and rhs>>=
  // 	c setup: matrix and rhs
#define A_PUT_ENTRY(m_j, m_k, m_x)\
  if ((m_j < 0) || (m_j >= A->m)) {\
    cout << "row index out of range at line " << __LINE__ << " in file "\
      __FILE__ << ".  Aborting." << endl;\
	exit(1); }\
  if ((m_k < 0) || (m_k >= A->n)) {\
    cout << "column index out of range at line " << __LINE__ << " in file "\
      __FILE__ << ".  Aborting." << endl;\
	exit(1); }\
  A->me[m_j][m_k] = m_x
#define B_PUT_ENTRY(m_j, m_x)\
  if ((m_j < 0) || (m_j >= B->dim)) {\
    cout << "row index out of range at line " << __LINE__ << " in file "\
      __FILE__ << ".  Aborting." << endl;\
	exit(1); }\
  B->ve[m_j] = m_x	  

  MAT *A;	// define
  VEC *B;
  A = m_get(m,m);	// allocate memory
  B = v_get(m);
@ %def A_PUT_ENTRY B_PUT_ENTRY
%
For these definitions, we further require some include files:
<<c setup: include files>>=
#include <stdio.h>
#include <stdlib.h>
#include <iostream.h>
#include <math.h>
 
extern "C" {
#include	"matrix2.h"
#include 	"matlab.h"
}
@
%
Next, the global variables $h, N$ and $m$. $N,M$ and $alpha$\footnote{The enpoint behavior, and the grid spacing factor} are read from the file 
\verb|unbndd-1d-ipl-global-vars.in|, $h$ and  $m$ are calculated.
The global definition:
<<c setup: global variable declaration>>=
//	c setup: global variable declaration
double alpha, h, pi;
int M, N, m;
@
%
And the value assignment:
<<c setup: global variable assignment>>=
//	c setup: global variable assignment

  pi = 3.14159265358979323846264338328;

  FILE *fp;
  if ((fp = fopen("unbndd-1d-ipl-global-vars.in", "r")) == NULL)
    {
      printf("unable to open input file %s \n", "unbndd-1d-ipl-global-vars.in");
      exit(1);
    }
  fscanf(fp, "%lf %d %*[^\n]", &(alpha), &(N));
  fscanf(fp, "%d %*[^\n]",  &(M));
  fclose(fp);

  h = pi/sqrt(alpha*(double)(N));
  m = M+N+1;

@
% 
For other files, the global variables will be external; we can put the
declarations in a common header file, along with other needed things:
<<unbndd-1d-ipl-stinc.h>>=
#ifndef MAIN_PROG
extern double alpha, h;
extern int N, M, m;
#endif

#include <math.h>
@
%
 

%
The input file thus will look like
<<unbndd-1d-ipl-global-vars.in>>=
1.0   10	alpha N
11		M

@
%
Lastly, the functions \verb|x(j), skh_pri(k,h,x), f_pri(x)| and
associated ones.  These
will be generated  by {\sc maple} and put into the file  
\verb|unbndd-1d-ipl-setup-subs.cc|.  Their prototypes go into
\verb|unbndd-1d-ipl-setup-subs.h|.
<<unbndd-1d-ipl-maple-setup.map>>=
#
Digits := 16;
interface(errorbreak = 2);	# abort on error
read `sinc8.mas`;

<<maple subs setup: define some needed derivatives>>

<<maple subs setup: define functions>>

<<maple subs setup: generate C++ code and header files>>
@
And here are the invokations used in the Makefile to obtain the
generated C++ code:
<<Makefile: get C++ setup subs from MAPLE>>=

unbndd-1d-ipl-setup-subs.cc:    unbndd-1d-ipl-maple-setup.map unbndd-1d-ipl-setup-subs-1.h
	(echo "gc(0);"; echo "words(0);"; \
	cat unbndd-1d-ipl-maple-setup.map ; echo "quit;") | maple > $@.mout

unbndd-1d-ipl-setup-subs.h:    unbndd-1d-ipl-maple-setup.map  unbndd-1d-ipl-setup-subs-1.h
	(echo "gc(0);"; echo "words(0);"; \
	cat unbndd-1d-ipl-maple-setup.map ; echo "quit;") | maple > $@.mout
@
%
To keep the size of generated code  to a minimum, the above maple
routines do not expand $\phi, \psi$ or $x_{j}$.  These functions are
done separately by {\sc maple}  and are put into
\verb|unbndd-1d-ipl-setup-subs-1.cc|;  their prototypes go into
\verb|unbndd-1d-ipl-setup-subs-1.h|.
%
<<unbndd-1d-ipl-maple-setup-1.map>>=
Digits := 16;
interface(errorbreak = 2);	# abort on error
read `sinc8.mas`;

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

phi_p := x-> diff(phi(x),x):
psi_p := x-> diff(psi(x),x):

psi(j*h);
xj := unapply(", j);

cff:= proc(expr);		# factor, convert to float (no factor)
   convert((expr), float);
end:
if (macroc_loaded <> true) then;
   read `macroc`;
fi;
optimized := true;
precision := double;
autodeclare := double;

OUTPUT_FILE := `unbndd-1d-ipl-setup-subs-1.cc`;
PROTOTYPE_FILE := `unbndd-1d-ipl-setup-subs-1.h`;
WHICH_UNKNOWN := ` `;

prog := [
      [textC, `#include "unbndd-1d-ipl-stinc.h"`]
      ];
start_fun(prog);
prot := [textC, `// prototypes from unbndd-1d-ipl-maple-setup-1.map`];
start_proto(prot);

function_list := [
      [`xj`, `j`, `[[int, j]]`],
      [`phi`, `x`, `[[double, x]]`],
      [`psi`, `x`, `[[double, x]]`],
      [`psi_p`, `x`, `[[double, x]]`],
      [`phi_p`, `x`, `[[double, x]]`]
      ];
gen_func(function_list);
gen_proto(function_list);
interface(errorbreak = 1);  
@
To actually get the C++ code out of this {\sc maple} listing, the
following Makefile entry is used:
<<Makefile: get C++ setup subs 1 from MAPLE>>=
unbndd-1d-ipl-setup-subs-1.cc:    unbndd-1d-ipl-maple-setup-1.map
	(echo "gc(0);"; echo "words(0);"; \
	cat $? ; echo "quit;") | maple > $@.mout
unbndd-1d-ipl-setup-subs-1.h:    unbndd-1d-ipl-maple-setup-1.map
	(echo "gc(0);"; echo "words(0);"; \
	cat $? ; echo "quit;") | maple > $@.mout
@
To ensure inclusion of the  prototypes generated above:
<<unbndd-1d-ipl-stinc.h>>=
#include "unbndd-1d-ipl-setup-subs-1.h"
@
This will also be needed in the main C code.  Thus, add it:
<<c setup: include files>>=
#include "unbndd-1d-ipl-stinc.h"
@
%
The extra {\sc maple} commands keep messages from corrupting the
generated code.
%
The only functions left to make are the sinc functions and their
derivatives.  Because this requires some care for small argument
values\footnote{For $x$ near zero, sinc(x) is best evaluated using a
taylor series.}, these routines are done manually:
<<c setup: sinc functions/derivatives>>=
double sinc(double x)		
{
  if (fabs(x) < 0.004) 
    return 1.0-PI*PI*x*x/6+pow(PI,4.0)*pow(x,4.0)/120;
  else
    return sin(PI*x)/(PI*x);
} 

double sinc_p(double x) 	// first derivative
{
  if (fabs(x) < 0.004)
    return -PI*PI*x/3+pow(PI,4.0)*x*x*x/30-pow(PI,6.0)*pow(x,5.0)/840;
  else
    return cos(PI*x)/x-sin(PI*x)/PI/(x*x);
} 

@
%
To make these definitions available globally, include the prototypes
in 
<<unbndd-1d-ipl-stinc.h>>=
double sinc(double x);
double sinc_p(double x);
@

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
In the following, we define all  functions needed in the main system loop. 
<<maple subs setup: define functions>>=
skh := unapply(sinc((phi(x)-k*h)/h), (k,h,x));
skh_pri := unapply(diff(skh(k,h,x), x), (k,h,x));

f := unapply(x**(1/2)*(1-x), x);
f_pri := unapply(diff(f(x),x), x);

@
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Next, we generate the C++ code and headers;  the listing
is followed by the appropriate Makefile entry to invoke maple.
%
<<maple subs setup: generate C++ code and header files>>=
# 	Preliminaries:
#
OUTPUT_FILE := `unbndd-1d-ipl-setup-subs.cc`;
PROTOTYPE_FILE := `unbndd-1d-ipl-setup-subs.h`;
WHICH_UNKNOWN := ` `;
#
cff:= proc(expr);		# factor, convert to float (no factor)
   convert((expr), float);
end:
if (macroc_loaded <> true) then;
   read `macroc`;
fi;
optimized := true;
precision := double;
autodeclare := double;

#	Subroutines:
#	~~~~~~~~~~~~
#
# 	startup:
prog := [
      [textC, `#include "unbndd-1d-ipl-stinc.h"`]
      ];
start_fun(prog);
prot := [textC, `// prototypes from unbndd-1d-ipl-maple-setup.map`];
start_proto(prot);

function_list := [
      [`f`, `x`, `[[double, x]]`],
      [`f_pri`, `x`, `[[double, x]]`],
      [`skh`, `k,h,x`, `[[int, k], [double, h], [double, x]]`],
      [`skh_pri`, `k,h,x`, `[[int, k], [double, h], [double, x]]`]
      ];
gen_func(function_list);
gen_proto(function_list);
interface(errorbreak = 1);
@
Since the prototypes generated above will be needed by most files, add
them to 
<<unbndd-1d-ipl-stinc.h>>=
#include "unbndd-1d-ipl-setup-subs.h"
@
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
<<maple subs setup: define some needed derivatives>>=
# Defs for derivatives (psi_p(), etc.)
`diff/psi` := proc(a,x);
   psi_p(a)*diff(a,x);
end:
`diff/phi` := proc(a,x);
   phi_p(a)*diff(a,x);
end:
`diff/psi_p` := proc(a,x);
   psi_pp(a)*diff(a,x);
end:
`diff/phi_p` := proc(a,x);
   phi_pp(a)*diff(a,x);
end:
@
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
With the interpolation/collocation matrix thus set up, it only remains
to solve the system and save the answer vector obtained.  Thus,
<<c setup: solve system>>=
  PERM	*pivot;
  VEC *u_vec;
  pivot = px_get(m);
  LUfactor(A, pivot);
  u_vec = LUsolve(A, pivot, B, VNULL);
@
The vector obtained will be saved in the file unbndd-1d-ipl-answer.mat:
%
<<c setup: save answer vector>>=
{
  FILE *fp;
  if ( (fp=fopen("unbndd-1d-ipl-answer.mat","w")) == (FILE *) NULL ) {
    printf("Cannot open save file \"%s\"\n", "unbndd-1d-ipl-answer.mat");
    exit(1);
  }

  v_save(fp, u_vec, "u_vec");

  fclose(fp);

}
@
It may be useful to have the matrix and right hand side available in
matlab; for this, we use
<<c setup: save matrix and rhs>>=
{
  FILE *fp;
  if ( (fp=fopen("unbndd-1d-ipl-A_and_b.mat","w")) == (FILE *) NULL ) {
    printf("Cannot open save file \"%s\"\n", "unbndd-1d-ipl-A_and_B.mat");
    exit(1);
  }

  v_save(fp, B, "b");
  m_save(fp, A, "A");

  fclose(fp);
}
@
%
We thus have the following main program:
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
<<unbndd-1d-ipl-setup.cc>>=
#define MAIN_PROG
<<c setup: include files>>
<<c setup: global variable declaration>>
<<c setup: sinc functions/derivatives>>
int main()  {
<<c setup: global variable assignment>>
<<c setup: matrix and rhs>>
<<c setup: main loop>>
<<c setup: save matrix and rhs>>
<<c setup: solve system>>
<<c setup: save answer vector>>
}
@
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% 
To properly run the {\sc maple } code and compile the subsequent C++
code, a Makefile is put in order:
<<Makefile>>=
cc = gcc

CFLAGS = -g -I/root/meschach
LIB = /root/meschach/meschach.a -lm

.SUFFIXES:	.map .cc .h

<<Makefile: get C++ setup subs from MAPLE>>

<<Makefile: get C++ setup subs 1 from MAPLE>>

# 	main program:
D1 = unbndd-1d-ipl-setup.o unbndd-1d-ipl-setup-subs-1.o \
	unbndd-1d-ipl-setup-subs.o
unbndd-1d-ipl-setup.exe:	$(D1) 
	$(CXX) -o $@ $(LFLAGS) $(D1) $(LIB)

.cc.o:
	$(CXX) -c $*.cc $(CFLAGS) 

@
%
The above will set up and solve the system for equation
(\ref{eqn:derivative-collocation}).  It now remains to use the answer
vector obtained to calculate the interpolant for $f$ and $f'$ at
arbitrary points in the domain.  This will be most flexible if done in
two steps:  (a) making a list of plot points in {\sc matlab} (or a
clone thereof) and saving it in a file;  (b) running a C++ program to
read the list, calculate all functions, and save results in {\sc
matlab}-readable format for plotting etc.

A simple {\sc octave} code to save the plot points, run the plot
program and read the results is
%
<<Unbndd1dIplPlotInput.m>>=
%       save plot points as column vector
x = [0.01:0.01:0.99]';
save -mat-binary unbndd-1d-ipl-plot-data.mat x

%	do the plotting
system("unbndd-1d-ipl-plot.exe")

%	read output
load -force unbndd-1d-ipl-plot-results.mat
%	release 1.1.1 of octave has braindamage reading matlab files 
%	in row-major order, so adjust results:
results = reshape(results, size(results, 2), size(results, 1))';

@
%
The C++ code to do the actual work will have a structure very similar
to the setup program:
<<unbndd-1d-ipl-plot.cc>>=
#define MAIN_PROG
<<c setup: include files>>
<<c setup: global variable declaration>>
<<c setup: sinc functions/derivatives>>
int main()  {
<<c setup: global variable assignment>>
<<c plot: read plot points>>
<<c plot: main loop>>
<<c plot: save results >>
}
@
Thus, only three more code segments are required.  But first, the
Makefile entry:
<<Makefile>>=
#
# 	plotting program:
DPLOT = unbndd-1d-ipl-plot.o unbndd-1d-ipl-setup-subs-1.o \
	unbndd-1d-ipl-setup-subs.o
unbndd-1d-ipl-plot.exe:	$(DPLOT) 
	$(CXX) -o $@ $(LFLAGS) $(DPLOT) $(LIB)
#
@
%
Also, the names are somewhat long, so add an ``all'' target:
<<Makefile>>=
#
all:	unbndd-1d-ipl-plot.exe unbndd-1d-ipl-setup.exe

@
%
Now the remaining code chunks.
<<c plot: read plot points>>=

MAT *points;
char *name;

{
  FILE *fp;
  if ( (fp=fopen("unbndd-1d-ipl-plot-data.mat","r")) == (FILE *) NULL ) {
    printf("Cannot open input file \"%s\"\n", "unbndd-1d-ipl-plot-data.mat");
    exit(1);
  }

  points = m_load(fp, &name);

  if (points->n != 1) {
    cout << "Error: unbndd-1d-ipl-plot-data.mat does not contain "\
      "a column vector.  Aborting." << endl;
    exit(1); 
  }

  fclose(fp);
}

@ %def points name
%
For running the main loop, we need the previously computed answer vector.
<<c plot: main loop>>=
//	define and read the answer vector:
MAT *u_vec;
{
  FILE *fp;
  if ( (fp=fopen("unbndd-1d-ipl-answer.mat","r")) == (FILE *) NULL ) {
    printf("Cannot open input file \"%s\"\n", "unbndd-1d-ipl-answer.mat");
    exit(1);
  }

  u_vec = m_load(fp, &name);

  if (u_vec->n != 1) {
    cout << "Error: unbndd-1d-ipl-answer.mat does not contain "\
      "a column vector.  Aborting." << endl;
    exit(1); 
  }
  if (u_vec->m != m) {
    cout << "Error: size of answer vector in unbndd-1d-ipl-answer.mat" \
      " is not compatible with global variable m.  Aborting." << endl;
    exit(1);
  }

  fclose(fp);
}
@
Next, create a matrix with 5 columns to hold, in order, $x, f, $ skh,
$ f',$ skh\_pri.
%
<<c plot: main loop>>=
//
//	do plotting, put results in matrix:
MAT *results;
results = m_get(points->m, 5);
{

int i,j, k;
for (i=0;  i<points->m;  i++) { 	// cover plot points
  double df = 0;
  double dfp = 0;
  for (j=-M;  j<=N;  j++) {
    df  += skh(j, h, points->me[i][0] ) * u_vec->me[j+M][0];
    dfp += skh_pri(j, h, points->me[i][0] ) * u_vec->me[j+M][0] ;
  }	// for j
  results->me[i][0] = points->me[i][0];
  results->me[i][1] = f(points->me[i][0]);
  results->me[i][2] = df;
  results->me[i][3] = f_pri(points->me[i][0]);
  results->me[i][4] = dfp;
}	// for i
}
@
%
<<c plot: save results >>=
{
  FILE *fp;
  if ( (fp=fopen("unbndd-1d-ipl-plot-results.mat","w")) == (FILE *) NULL ) {
    printf("Cannot open save file \"%s\"\n", "unbndd-1d-ipl-plot-results.mat");
    exit(1);
  }

  m_save(fp, results, "results");

  fclose(fp);
}
@
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%               RESULTS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Some test runs and results}
For $f(x) = \sqrt{x}(1-x), \alpha = 1.0, $ we have the following results.
\begin{itemize}
  \item  First derivative matricies of odd order are
        singular, and the solution obtained via direct
        elimination is useless -- this would reflect in the
        2D problem as well.  Of course, with $N=M$, the
        matrix is always of odd order, so using $M=N+1$ (or
        a similar adjustment) may help the 2D case also.  
        
  \item The least squares solution obtained via the SVD for odd order
        matrices has very oscillatory  results near the endpoints for
        the derivative, and a slight oscillation and an offset for $f$
        itself. 

  \item Using even-order matricies, here with $M=N+1=11$, there is no
        visual difference between the interpolants and the functions
        for $x=0.01:0.01:0.99$.   Numerically, the relative error for
        $f$'s interpolant is less than 1\%, except for $x=0.99$, where
        it is 12\%.  For $f'$, the error is less than 5\%, uniformly.
\end{itemize}

For $\displaystyle f(x) = \sqrt{x}(1-x)^{4}, \alpha = 1.0, M=N+1=11,$
we get very good visual results, although the relative error
deteriorates severely near the fourth order zero.  The derivative
approximation has a small jump near $x=1$ as well.  \\ 
With $\alpha = 0.5$, the endpoint results improve somewhat, but the
interior is not nearly as good -- as expected.

For $\displaystyle f(x) = x^{1/8}(1-x)^{2}, \alpha = 0.5, M=N+1=11, $
there is an offset in the approximation of $f$, and oscillations at
the endpoints for $f'$. \\
Using $\alpha = 1.0$, the results are much worse;  worse yet for
$\alpha = 2.0$. 
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
\section{Indicies}
\subsection{Code Chunks}
\nowebchunks
\subsection{Identifiers}
\nowebindex

% 
% to extract code from this file, do a  
%       make unbndd-1d-ipl-all-code
%
% to compile the resulting mess, do  
%      make -f Makefile all
% 
% to run, do
%       unbndd-1d-ipl-setup.exe  	for every new parameter set
%
% use Unbndd1dIplPlotInput.m in octave for the plotting.
%