\documentstyle[noweb]{article}\pagestyle{noweb}\noweboptions{}\begin{document}\nwfilename{unbndd-1d-ipl-code.nw}\nwbegindocs{0}\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 \nwenddocs{}\nwbegincode{1}\sublabel{NWunbL-c*sI-1}\nwmargintag{{\nwtagstyle{}\subpageref{NWunbL-c*sI-1}}}\moddef{c setup: main loop~{\nwtagstyle{}\subpageref{NWunbL-c*sI-1}}}\endmoddef \{ // 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 \} \nwidentuses{\\{{A{\char95}PUT{\char95}ENTRY}{A:unPUT:unENTRY}}\\{{B{\char95}PUT{\char95}ENTRY}{B:unPUT:unENTRY}}}\nwindexuse{A{\char95}PUT{\char95}ENTRY}{A:unPUT:unENTRY}{NWunbL-c*sI-1}\nwindexuse{B{\char95}PUT{\char95}ENTRY}{B:unPUT:unENTRY}{NWunbL-c*sI-1}\nwused{\\{NWunbL-unbM-1}}\nwendcode{}\nwbegindocs{2}\nwdocspar % This code segment requires some definitions and functions. First, the matrix and rhs vector, and the access macros. \nwenddocs{}\nwbegincode{3}\sublabel{NWunbL-c*sN-1}\nwmargintag{{\nwtagstyle{}\subpageref{NWunbL-c*sN-1}}}\moddef{c setup: matrix and rhs~{\nwtagstyle{}\subpageref{NWunbL-c*sN-1}}}\endmoddef // 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); \nwindexdefn{A{\char95}PUT{\char95}ENTRY}{A:unPUT:unENTRY}{NWunbL-c*sN-1}\nwindexdefn{B{\char95}PUT{\char95}ENTRY}{B:unPUT:unENTRY}{NWunbL-c*sN-1}\eatline \nwidentdefs{\\{{A{\char95}PUT{\char95}ENTRY}{A:unPUT:unENTRY}}\\{{B{\char95}PUT{\char95}ENTRY}{B:unPUT:unENTRY}}}\nwused{\\{NWunbL-unbM-1}}\nwendcode{}\nwbegindocs{4}% For these definitions, we further require some include files: \nwenddocs{}\nwbegincode{5}\sublabel{NWunbL-c*sM-1}\nwmargintag{{\nwtagstyle{}\subpageref{NWunbL-c*sM-1}}}\moddef{c setup: include files~{\nwtagstyle{}\subpageref{NWunbL-c*sM-1}}}\endmoddef #include #include #include #include extern "C" \{ #include "matrix2.h" #include "matlab.h" \} \nwalsodefined{\\{NWunbL-c*sM-2}}\nwused{\\{NWunbL-unbM-1}\\{NWunbL-unbL.3-1}}\nwendcode{}\nwbegindocs{6}\nwdocspar % 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: \nwenddocs{}\nwbegincode{7}\sublabel{NWunbL-c*sa-1}\nwmargintag{{\nwtagstyle{}\subpageref{NWunbL-c*sa-1}}}\moddef{c setup: global variable declaration~{\nwtagstyle{}\subpageref{NWunbL-c*sa-1}}}\endmoddef // c setup: global variable declaration double alpha, h, pi; int M, N, m; \nwused{\\{NWunbL-unbM-1}\\{NWunbL-unbL.3-1}}\nwendcode{}\nwbegindocs{8}\nwdocspar % And the value assignment: \nwenddocs{}\nwbegincode{9}\sublabel{NWunbL-c*sZ-1}\nwmargintag{{\nwtagstyle{}\subpageref{NWunbL-c*sZ-1}}}\moddef{c setup: global variable assignment~{\nwtagstyle{}\subpageref{NWunbL-c*sZ-1}}}\endmoddef // 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; \nwused{\\{NWunbL-unbM-1}\\{NWunbL-unbL.3-1}}\nwendcode{}\nwbegindocs{10}\nwdocspar % For other files, the global variables will be external; we can put the declarations in a common header file, along with other needed things: \nwenddocs{}\nwbegincode{11}\sublabel{NWunbL-unbL.2-1}\nwmargintag{{\nwtagstyle{}\subpageref{NWunbL-unbL.2-1}}}\moddef{unbndd-1d-ipl-stinc.h~{\nwtagstyle{}\subpageref{NWunbL-unbL.2-1}}}\endmoddef #ifndef MAIN_PROG extern double alpha, h; extern int N, M, m; #endif #include \nwalsodefined{\\{NWunbL-unbL.2-2}\\{NWunbL-unbL.2-3}\\{NWunbL-unbL.2-4}}\nwnotused{unbndd-1d-ipl-stinc.h}\nwendcode{}\nwbegindocs{12}\nwdocspar % % The input file thus will look like \nwenddocs{}\nwbegincode{13}\sublabel{NWunbL-unbS-1}\nwmargintag{{\nwtagstyle{}\subpageref{NWunbL-unbS-1}}}\moddef{unbndd-1d-ipl-global-vars.in~{\nwtagstyle{}\subpageref{NWunbL-unbS-1}}}\endmoddef 1.0 10 alpha N 11 M \nwnotused{unbndd-1d-ipl-global-vars.in}\nwendcode{}\nwbegindocs{14}\nwdocspar % 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|. \nwenddocs{}\nwbegincode{15}\sublabel{NWunbL-unbT-1}\nwmargintag{{\nwtagstyle{}\subpageref{NWunbL-unbT-1}}}\moddef{unbndd-1d-ipl-maple-setup.map~{\nwtagstyle{}\subpageref{NWunbL-unbT-1}}}\endmoddef # Digits := 16; interface(errorbreak = 2); # abort on error read `sinc8.mas`; \LA{}maple subs setup: define some needed derivatives~{\nwtagstyle{}\subpageref{NWunbL-mapm-1}}\RA{} \LA{}maple subs setup: define functions~{\nwtagstyle{}\subpageref{NWunbL-mapY-1}}\RA{} \LA{}maple subs setup: generate C++ code and header files~{\nwtagstyle{}\subpageref{NWunbL-mapq-1}}\RA{} \nwnotused{unbndd-1d-ipl-maple-setup.map}\nwendcode{}\nwbegindocs{16}\nwdocspar And here are the invokations used in the Makefile to obtain the generated C++ code: \nwenddocs{}\nwbegincode{17}\sublabel{NWunbL-Makd-1}\nwmargintag{{\nwtagstyle{}\subpageref{NWunbL-Makd-1}}}\moddef{Makefile: get C++ setup subs from MAPLE~{\nwtagstyle{}\subpageref{NWunbL-Makd-1}}}\endmoddef 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 \nwused{\\{NWunbL-Mak8-1}}\nwendcode{}\nwbegindocs{18}\nwdocspar % 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|. % \nwenddocs{}\nwbegincode{19}\sublabel{NWunbL-unbV-1}\nwmargintag{{\nwtagstyle{}\subpageref{NWunbL-unbV-1}}}\moddef{unbndd-1d-ipl-maple-setup-1.map~{\nwtagstyle{}\subpageref{NWunbL-unbV-1}}}\endmoddef 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); \nwnotused{unbndd-1d-ipl-maple-setup-1.map}\nwendcode{}\nwbegindocs{20}\nwdocspar To actually get the C++ code out of this {\sc maple} listing, the following Makefile entry is used: \nwenddocs{}\nwbegincode{21}\sublabel{NWunbL-Makf-1}\nwmargintag{{\nwtagstyle{}\subpageref{NWunbL-Makf-1}}}\moddef{Makefile: get C++ setup subs 1 from MAPLE~{\nwtagstyle{}\subpageref{NWunbL-Makf-1}}}\endmoddef 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 \nwused{\\{NWunbL-Mak8-1}}\nwendcode{}\nwbegindocs{22}\nwdocspar To ensure inclusion of the prototypes generated above: \nwenddocs{}\nwbegincode{23}\sublabel{NWunbL-unbL.2-2}\nwmargintag{{\nwtagstyle{}\subpageref{NWunbL-unbL.2-2}}}\moddef{unbndd-1d-ipl-stinc.h~{\nwtagstyle{}\subpageref{NWunbL-unbL.2-1}}}\plusendmoddef #include "unbndd-1d-ipl-setup-subs-1.h" \nwendcode{}\nwbegindocs{24}\nwdocspar This will also be needed in the main C code. Thus, add it: \nwenddocs{}\nwbegincode{25}\sublabel{NWunbL-c*sM-2}\nwmargintag{{\nwtagstyle{}\subpageref{NWunbL-c*sM-2}}}\moddef{c setup: include files~{\nwtagstyle{}\subpageref{NWunbL-c*sM-1}}}\plusendmoddef #include "unbndd-1d-ipl-stinc.h" \nwendcode{}\nwbegindocs{26}\nwdocspar % 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: \nwenddocs{}\nwbegincode{27}\sublabel{NWunbL-derZ-1}\nwmargintag{{\nwtagstyle{}\subpageref{NWunbL-derZ-1}}}\moddef{c setup: sinc functions/derivatives~{\nwtagstyle{}\subpageref{NWunbL-derZ-1}}}\endmoddef 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); \} \nwused{\\{NWunbL-unbM-1}\\{NWunbL-unbL.3-1}}\nwendcode{}\nwbegindocs{28}\nwdocspar % To make these definitions available globally, include the prototypes in \nwenddocs{}\nwbegincode{29}\sublabel{NWunbL-unbL.2-3}\nwmargintag{{\nwtagstyle{}\subpageref{NWunbL-unbL.2-3}}}\moddef{unbndd-1d-ipl-stinc.h~{\nwtagstyle{}\subpageref{NWunbL-unbL.2-1}}}\plusendmoddef double sinc(double x); double sinc_p(double x); \nwendcode{}\nwbegindocs{30}\nwdocspar %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% In the following, we define all functions needed in the main system loop. \nwenddocs{}\nwbegincode{31}\sublabel{NWunbL-mapY-1}\nwmargintag{{\nwtagstyle{}\subpageref{NWunbL-mapY-1}}}\moddef{maple subs setup: define functions~{\nwtagstyle{}\subpageref{NWunbL-mapY-1}}}\endmoddef 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); \nwused{\\{NWunbL-unbT-1}}\nwendcode{}\nwbegindocs{32}\nwdocspar % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Next, we generate the C++ code and headers; the listing is followed by the appropriate Makefile entry to invoke maple. % \nwenddocs{}\nwbegincode{33}\sublabel{NWunbL-mapq-1}\nwmargintag{{\nwtagstyle{}\subpageref{NWunbL-mapq-1}}}\moddef{maple subs setup: generate C++ code and header files~{\nwtagstyle{}\subpageref{NWunbL-mapq-1}}}\endmoddef # 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); \nwused{\\{NWunbL-unbT-1}}\nwendcode{}\nwbegindocs{34}\nwdocspar Since the prototypes generated above will be needed by most files, add them to \nwenddocs{}\nwbegincode{35}\sublabel{NWunbL-unbL.2-4}\nwmargintag{{\nwtagstyle{}\subpageref{NWunbL-unbL.2-4}}}\moddef{unbndd-1d-ipl-stinc.h~{\nwtagstyle{}\subpageref{NWunbL-unbL.2-1}}}\plusendmoddef #include "unbndd-1d-ipl-setup-subs.h" \nwendcode{}\nwbegindocs{36}\nwdocspar % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \nwenddocs{}\nwbegincode{37}\sublabel{NWunbL-mapm-1}\nwmargintag{{\nwtagstyle{}\subpageref{NWunbL-mapm-1}}}\moddef{maple subs setup: define some needed derivatives~{\nwtagstyle{}\subpageref{NWunbL-mapm-1}}}\endmoddef # 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: \nwused{\\{NWunbL-unbT-1}}\nwendcode{}\nwbegindocs{38}\nwdocspar % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % With the interpolation/collocation matrix thus set up, it only remains to solve the system and save the answer vector obtained. Thus, \nwenddocs{}\nwbegincode{39}\sublabel{NWunbL-c*sL-1}\nwmargintag{{\nwtagstyle{}\subpageref{NWunbL-c*sL-1}}}\moddef{c setup: solve system~{\nwtagstyle{}\subpageref{NWunbL-c*sL-1}}}\endmoddef PERM *pivot; VEC *u_vec; pivot = px_get(m); LUfactor(A, pivot); u_vec = LUsolve(A, pivot, B, VNULL); \nwused{\\{NWunbL-unbM-1}}\nwendcode{}\nwbegindocs{40}\nwdocspar The vector obtained will be saved in the file unbndd-1d-ipl-answer.mat: % \nwenddocs{}\nwbegincode{41}\sublabel{NWunbL-c*sR-1}\nwmargintag{{\nwtagstyle{}\subpageref{NWunbL-c*sR-1}}}\moddef{c setup: save answer vector~{\nwtagstyle{}\subpageref{NWunbL-c*sR-1}}}\endmoddef \{ 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); \} \nwused{\\{NWunbL-unbM-1}}\nwendcode{}\nwbegindocs{42}\nwdocspar It may be useful to have the matrix and right hand side available in matlab; for this, we use \nwenddocs{}\nwbegincode{43}\sublabel{NWunbL-c*sS-1}\nwmargintag{{\nwtagstyle{}\subpageref{NWunbL-c*sS-1}}}\moddef{c setup: save matrix and rhs~{\nwtagstyle{}\subpageref{NWunbL-c*sS-1}}}\endmoddef \{ 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); \} \nwused{\\{NWunbL-unbM-1}}\nwendcode{}\nwbegindocs{44}\nwdocspar % We thus have the following main program: % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % \nwenddocs{}\nwbegincode{45}\sublabel{NWunbL-unbM-1}\nwmargintag{{\nwtagstyle{}\subpageref{NWunbL-unbM-1}}}\moddef{unbndd-1d-ipl-setup.cc~{\nwtagstyle{}\subpageref{NWunbL-unbM-1}}}\endmoddef #define MAIN_PROG \LA{}c setup: include files~{\nwtagstyle{}\subpageref{NWunbL-c*sM-1}}\RA{} \LA{}c setup: global variable declaration~{\nwtagstyle{}\subpageref{NWunbL-c*sa-1}}\RA{} \LA{}c setup: sinc functions/derivatives~{\nwtagstyle{}\subpageref{NWunbL-derZ-1}}\RA{} int main() \{ \LA{}c setup: global variable assignment~{\nwtagstyle{}\subpageref{NWunbL-c*sZ-1}}\RA{} \LA{}c setup: matrix and rhs~{\nwtagstyle{}\subpageref{NWunbL-c*sN-1}}\RA{} \LA{}c setup: main loop~{\nwtagstyle{}\subpageref{NWunbL-c*sI-1}}\RA{} \LA{}c setup: save matrix and rhs~{\nwtagstyle{}\subpageref{NWunbL-c*sS-1}}\RA{} \LA{}c setup: solve system~{\nwtagstyle{}\subpageref{NWunbL-c*sL-1}}\RA{} \LA{}c setup: save answer vector~{\nwtagstyle{}\subpageref{NWunbL-c*sR-1}}\RA{} \} \nwnotused{unbndd-1d-ipl-setup.cc}\nwendcode{}\nwbegindocs{46}\nwdocspar % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % To properly run the {\sc maple } code and compile the subsequent C++ code, a Makefile is put in order: \nwenddocs{}\nwbegincode{47}\sublabel{NWunbL-Mak8-1}\nwmargintag{{\nwtagstyle{}\subpageref{NWunbL-Mak8-1}}}\moddef{Makefile~{\nwtagstyle{}\subpageref{NWunbL-Mak8-1}}}\endmoddef cc = gcc CFLAGS = -g -I/root/meschach LIB = /root/meschach/meschach.a -lm .SUFFIXES: .map .cc .h \LA{}Makefile: get C++ setup subs from MAPLE~{\nwtagstyle{}\subpageref{NWunbL-Makd-1}}\RA{} \LA{}Makefile: get C++ setup subs 1 from MAPLE~{\nwtagstyle{}\subpageref{NWunbL-Makf-1}}\RA{} # 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) \nwalsodefined{\\{NWunbL-Mak8-2}\\{NWunbL-Mak8-3}}\nwnotused{Makefile}\nwendcode{}\nwbegindocs{48}\nwdocspar % 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 % \nwenddocs{}\nwbegincode{49}\sublabel{NWunbL-UnbM-1}\nwmargintag{{\nwtagstyle{}\subpageref{NWunbL-UnbM-1}}}\moddef{Unbndd1dIplPlotInput.m~{\nwtagstyle{}\subpageref{NWunbL-UnbM-1}}}\endmoddef % 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))'; \nwidentuses{\\{{points}{points}}}\nwindexuse{points}{points}{NWunbL-UnbM-1}\nwnotused{Unbndd1dIplPlotInput.m}\nwendcode{}\nwbegindocs{50}\nwdocspar % The C++ code to do the actual work will have a structure very similar to the setup program: \nwenddocs{}\nwbegincode{51}\sublabel{NWunbL-unbL.3-1}\nwmargintag{{\nwtagstyle{}\subpageref{NWunbL-unbL.3-1}}}\moddef{unbndd-1d-ipl-plot.cc~{\nwtagstyle{}\subpageref{NWunbL-unbL.3-1}}}\endmoddef #define MAIN_PROG \LA{}c setup: include files~{\nwtagstyle{}\subpageref{NWunbL-c*sM-1}}\RA{} \LA{}c setup: global variable declaration~{\nwtagstyle{}\subpageref{NWunbL-c*sa-1}}\RA{} \LA{}c setup: sinc functions/derivatives~{\nwtagstyle{}\subpageref{NWunbL-derZ-1}}\RA{} int main() \{ \LA{}c setup: global variable assignment~{\nwtagstyle{}\subpageref{NWunbL-c*sZ-1}}\RA{} \LA{}c plot: read plot points~{\nwtagstyle{}\subpageref{NWunbL-c*pO-1}}\RA{} \LA{}c plot: main loop~{\nwtagstyle{}\subpageref{NWunbL-c*pH-1}}\RA{} \LA{}c plot: save results ~{\nwtagstyle{}\subpageref{NWunbL-c*pL-1}}\RA{} \} \nwnotused{unbndd-1d-ipl-plot.cc}\nwendcode{}\nwbegindocs{52}\nwdocspar Thus, only three more code segments are required. But first, the Makefile entry: \nwenddocs{}\nwbegincode{53}\sublabel{NWunbL-Mak8-2}\nwmargintag{{\nwtagstyle{}\subpageref{NWunbL-Mak8-2}}}\moddef{Makefile~{\nwtagstyle{}\subpageref{NWunbL-Mak8-1}}}\plusendmoddef # # 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) # \nwendcode{}\nwbegindocs{54}\nwdocspar % Also, the names are somewhat long, so add an ``all'' target: \nwenddocs{}\nwbegincode{55}\sublabel{NWunbL-Mak8-3}\nwmargintag{{\nwtagstyle{}\subpageref{NWunbL-Mak8-3}}}\moddef{Makefile~{\nwtagstyle{}\subpageref{NWunbL-Mak8-1}}}\plusendmoddef # all: unbndd-1d-ipl-plot.exe unbndd-1d-ipl-setup.exe \nwendcode{}\nwbegindocs{56}\nwdocspar % Now the remaining code chunks. \nwenddocs{}\nwbegincode{57}\sublabel{NWunbL-c*pO-1}\nwmargintag{{\nwtagstyle{}\subpageref{NWunbL-c*pO-1}}}\moddef{c plot: read plot points~{\nwtagstyle{}\subpageref{NWunbL-c*pO-1}}}\endmoddef 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); \} \nwindexdefn{points}{points}{NWunbL-c*pO-1}\nwindexdefn{name}{name}{NWunbL-c*pO-1}\eatline \nwidentdefs{\\{{name}{name}}\\{{points}{points}}}\nwused{\\{NWunbL-unbL.3-1}}\nwendcode{}\nwbegindocs{58}% For running the main loop, we need the previously computed answer vector. \nwenddocs{}\nwbegincode{59}\sublabel{NWunbL-c*pH-1}\nwmargintag{{\nwtagstyle{}\subpageref{NWunbL-c*pH-1}}}\moddef{c plot: main loop~{\nwtagstyle{}\subpageref{NWunbL-c*pH-1}}}\endmoddef // 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); \} \nwidentuses{\\{{name}{name}}}\nwindexuse{name}{name}{NWunbL-c*pH-1}\nwalsodefined{\\{NWunbL-c*pH-2}}\nwused{\\{NWunbL-unbL.3-1}}\nwendcode{}\nwbegindocs{60}\nwdocspar Next, create a matrix with 5 columns to hold, in order, $x, f, $ skh, $ f',$ skh\_pri. % \nwenddocs{}\nwbegincode{61}\sublabel{NWunbL-c*pH-2}\nwmargintag{{\nwtagstyle{}\subpageref{NWunbL-c*pH-2}}}\moddef{c plot: main loop~{\nwtagstyle{}\subpageref{NWunbL-c*pH-1}}}\plusendmoddef // // do plotting, put results in matrix: MAT *results; results = m_get(points->m, 5); \{ int i,j, k; for (i=0; im; 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 \} \nwidentuses{\\{{points}{points}}}\nwindexuse{points}{points}{NWunbL-c*pH-2}\nwendcode{}\nwbegindocs{62}\nwdocspar % \nwenddocs{}\nwbegincode{63}\sublabel{NWunbL-c*pL-1}\nwmargintag{{\nwtagstyle{}\subpageref{NWunbL-c*pL-1}}}\moddef{c plot: save results ~{\nwtagstyle{}\subpageref{NWunbL-c*pL-1}}}\endmoddef \{ 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); \} \nwused{\\{NWunbL-unbL.3-1}}\nwendcode{}\nwbegindocs{64}\nwdocspar % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 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. % \nwenddocs{} \nwixlogsorted{c}{{Makefile}{NWunbL-Mak8-1}{\nwixd{NWunbL-Mak8-1}\nwixd{NWunbL-Mak8-2}\nwixd{NWunbL-Mak8-3}}}% \nwixlogsorted{c}{{Makefile: get C++ setup subs 1 from MAPLE}{NWunbL-Makf-1}{\nwixd{NWunbL-Makf-1}\nwixu{NWunbL-Mak8-1}}}% \nwixlogsorted{c}{{Makefile: get C++ setup subs from MAPLE}{NWunbL-Makd-1}{\nwixd{NWunbL-Makd-1}\nwixu{NWunbL-Mak8-1}}}% \nwixlogsorted{c}{{Unbndd1dIplPlotInput.m}{NWunbL-UnbM-1}{\nwixd{NWunbL-UnbM-1}}}% \nwixlogsorted{c}{{c plot: main loop}{NWunbL-c*pH-1}{\nwixu{NWunbL-unbL.3-1}\nwixd{NWunbL-c*pH-1}\nwixd{NWunbL-c*pH-2}}}% \nwixlogsorted{c}{{c plot: read plot points}{NWunbL-c*pO-1}{\nwixu{NWunbL-unbL.3-1}\nwixd{NWunbL-c*pO-1}}}% \nwixlogsorted{c}{{c plot: save results }{NWunbL-c*pL-1}{\nwixu{NWunbL-unbL.3-1}\nwixd{NWunbL-c*pL-1}}}% \nwixlogsorted{c}{{c setup: global variable assignment}{NWunbL-c*sZ-1}{\nwixd{NWunbL-c*sZ-1}\nwixu{NWunbL-unbM-1}\nwixu{NWunbL-unbL.3-1}}}% \nwixlogsorted{c}{{c setup: global variable declaration}{NWunbL-c*sa-1}{\nwixd{NWunbL-c*sa-1}\nwixu{NWunbL-unbM-1}\nwixu{NWunbL-unbL.3-1}}}% \nwixlogsorted{c}{{c setup: include files}{NWunbL-c*sM-1}{\nwixd{NWunbL-c*sM-1}\nwixd{NWunbL-c*sM-2}\nwixu{NWunbL-unbM-1}\nwixu{NWunbL-unbL.3-1}}}% \nwixlogsorted{c}{{c setup: main loop}{NWunbL-c*sI-1}{\nwixd{NWunbL-c*sI-1}\nwixu{NWunbL-unbM-1}}}% \nwixlogsorted{c}{{c setup: matrix and rhs}{NWunbL-c*sN-1}{\nwixd{NWunbL-c*sN-1}\nwixu{NWunbL-unbM-1}}}% \nwixlogsorted{c}{{c setup: save answer vector}{NWunbL-c*sR-1}{\nwixd{NWunbL-c*sR-1}\nwixu{NWunbL-unbM-1}}}% \nwixlogsorted{c}{{c setup: save matrix and rhs}{NWunbL-c*sS-1}{\nwixd{NWunbL-c*sS-1}\nwixu{NWunbL-unbM-1}}}% \nwixlogsorted{c}{{c setup: sinc functions/derivatives}{NWunbL-derZ-1}{\nwixd{NWunbL-derZ-1}\nwixu{NWunbL-unbM-1}\nwixu{NWunbL-unbL.3-1}}}% \nwixlogsorted{c}{{c setup: solve system}{NWunbL-c*sL-1}{\nwixd{NWunbL-c*sL-1}\nwixu{NWunbL-unbM-1}}}% \nwixlogsorted{c}{{maple subs setup: define functions}{NWunbL-mapY-1}{\nwixu{NWunbL-unbT-1}\nwixd{NWunbL-mapY-1}}}% \nwixlogsorted{c}{{maple subs setup: define some needed derivatives}{NWunbL-mapm-1}{\nwixu{NWunbL-unbT-1}\nwixd{NWunbL-mapm-1}}}% \nwixlogsorted{c}{{maple subs setup: generate C++ code and header files}{NWunbL-mapq-1}{\nwixu{NWunbL-unbT-1}\nwixd{NWunbL-mapq-1}}}% \nwixlogsorted{c}{{unbndd-1d-ipl-global-vars.in}{NWunbL-unbS-1}{\nwixd{NWunbL-unbS-1}}}% \nwixlogsorted{c}{{unbndd-1d-ipl-maple-setup-1.map}{NWunbL-unbV-1}{\nwixd{NWunbL-unbV-1}}}% \nwixlogsorted{c}{{unbndd-1d-ipl-maple-setup.map}{NWunbL-unbT-1}{\nwixd{NWunbL-unbT-1}}}% \nwixlogsorted{c}{{unbndd-1d-ipl-plot.cc}{NWunbL-unbL.3-1}{\nwixd{NWunbL-unbL.3-1}}}% \nwixlogsorted{c}{{unbndd-1d-ipl-setup.cc}{NWunbL-unbM-1}{\nwixd{NWunbL-unbM-1}}}% \nwixlogsorted{c}{{unbndd-1d-ipl-stinc.h}{NWunbL-unbL.2-1}{\nwixd{NWunbL-unbL.2-1}\nwixd{NWunbL-unbL.2-2}\nwixd{NWunbL-unbL.2-3}\nwixd{NWunbL-unbL.2-4}}}% \nwixlogsorted{i}{{A{\char95}PUT{\char95}ENTRY}{A:unPUT:unENTRY}}% \nwixlogsorted{i}{{B{\char95}PUT{\char95}ENTRY}{B:unPUT:unENTRY}}% \nwixlogsorted{i}{{name}{name}}% \nwixlogsorted{i}{{points}{points}}% \end{document}