\documentstyle[12pt]{article}
\title{Collocation by Tchebycheffian splines of order 3}
\author{Mladen Rogina}
\date{}
\newcommand{\dj} {\mbox{\=o\hspace{-1mm}l}}
\begin{document}
\maketitle

\mbox{} \\ Math. Subj. Class. (1980) 41A15, 65D07, 65N03 \\[0.25in]
ABSTRACT:\newline
 \begin{em}   A  differential equation  of  the form  $Dp(x)Du(x)  
= f(x,u(x))$ is solved by collocation, using 
 a class of spline functions of order 3. A restriction  of  
each spline function to any subinterval of the partition is  piecewisely 
in the kernel of the operator $DDp(x)D$. A coefficient $p(x)$ may 
be  singular. An  algorithm is described by means of  which  a  B-
spline  basis  could be  efficiently  constructed. Some  numerical 
experiments are considered. \end{em}
\\[0.25in]
KOLOKACIJA \u{C}EBI\u{S}EVLJEVIM SPLINE FUNKCIJAMA REDA 3.~ \newline
\begin{em}    Diferencijalna  jednad\u{z}ba  oblika   $Dp(x)Du(x)  
=f(x,u(x))$  rje\u{s}ava se kolokacijom spline funkcijama reda  3. 
Restrikcija   svake  spline  funkcije  na  podintervale   je   po 
dijelovima element jezgre operatora $DDp(x)D$. Koeficijent $p(x)$ 
mo\u{z}e biti i singularan. Opisan je algoritam kojim se  mo\u{z}e 
efektivno  konstruirati baza B-spline-ova. Promatraju se  i  neki 
numeri\u{c}ki eksperimenti. \end{em}
\vspace{0.5in}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}{Lemma}[section]
\newtheorem{corollary}{Corollary}[section]

\mbox{}

\section{Introduction}
We want to solve the boundary value problem (or maybe an  initial 
value one) of the second order:
\begin{equation}
\frac{d}{dx}p(x)\frac{d}{dx}u(x)=f(x,u(x))
\label{eq:difequ}
\end{equation}
\[x\in(0,1)\]
together    with   suitable   boundary    conditions, such    that 
(~\ref{eq:difequ})  has  a  unique (isolated)  solution.  We  shall 
further  suppose  that  $p(x)$  is such  that  $p(x)\geq  0$  and 
$p(x)\frac{d}{dx}u(x)\in C^{1}[0,1]$. The idea is to project both 
sides  of the equation to some space of polynomial  splines,  the 
most  simple  space being the space of piecewise  constants.  For 
$p(x)\equiv1$   the   problem   was   extensively   studied    by 
\cite[deBoor]{deBoor:1} and many others. An important case arises 
when  an equation can not be easily converted to an  equation  of 
the  form \[\frac{d^{2}}{dx^{2}}u(x)=f(x,u(x))\] and thus we  can 
not  use  polynomial  splines  for  collocation.  This  situation 
frequently  arises  in solving  ordinary  differential  equations 
having  their  origin in elliptic-parabolic  partial  differential 
equations,  such  as in \cite{Parter}.  Usually  $p(x)$  exhibits 
singular  behavior near one or both end points of the  interval, 
such as for instance $p(0)=0$. A typical example is the equation
\begin{eqnarray}
D(x^{0.5}Du)-x^{0.5}u=0 \nonumber\\
0<x<1 \nonumber \\
u(0)=1 \hspace{0.5in} u(1)=0 
\label{eq:model}
\end{eqnarray}
which  we  shall later use for numerical  purposes,  the  example 
itself  being  borrowed  from  \cite[Jamet]{Jamet}, where a 
central   difference  scheme  to  approximate  the  solution   of 
(\ref{eq:difequ}) is constructed in the linear case.
 A variational approach may  be 
used in solving (\ref{eq:difequ}), as in \cite{Varga}, where  in fact  a special kind of spline space was  used  (Tchebycheffian 
"hat"  functions). In order to apply collocation techniques  to  a 
singular  problem, one must first find means of constructing  the 
singular  splines effectively, this task being already  achieved 
in the most satisfactory manner in case of polynomial splines 
in  the form of the well known recurrence \cite{deBoor:2}. In  what 
follows  we  propose a way of calculating  with  the  Tchebycheff 
splines of a very special kind (the polynomial analogue would  be 
the  space  of  parabolic  splines). Singular  splines  have  been 
treated in the most general way in \cite{Jerome} as solutions  of 
certain minimization problems, and in connection with collocation 
they have been used by Reddien \cite{Reddien:1}. 

\mbox{}

\section{Tchebycheff Systems}
Let us choose the knot sequence $t_{1} \ldots t_{m+2}$ in such  a 
way that
\begin{equation}
0=t_{1}=t_{2}=t_{3}<t_{4} \ldots <t_{m}=t_{m+1}=t_{m+2}=1.
\label{knots}
\end{equation}
We shall seek for an approximate solution of (\ref{eq:difequ}) in 
the  space  of  functions  satisfying  \[Dp(x)Du(x)=const.\]  for 
$x\in(t_{i},t_{i+1})$, $i=3 \ldots m-1$. This is equivalent to  a 
condition:
\begin{equation}
DDp(x)Du(x)=0 \hspace{0.5in}  x \in (t_{i},t_{i+1}).
\end {equation}
There are three basis functions in $Ker(DDp(x)D)$, namely:
\begin{eqnarray}
u_{1}=1 \nonumber\\
u_{2}=\int_{0}^{x} \frac{dt_{2}}{p(t_{2})} \label{primal}\\
u_{3}=\int_{0}^{x} \frac{dt_{2}}{p(t_{2})} \int_{0}^{t_{2}}dt_{3} 
     \nonumber   
\end{eqnarray} 
to  be used in construction of the local basis. We shall  further 
need reduced systems obtained by taking  generalized  derivatives 
of  the  function in (\ref{primal}). The  following  notation  is 
useful:
\[\begin{array}{ll}
L_{1}:=p(x)D & L_{1}^{*}:=D \\
L_{2}:=Dp(x)D & L_{2}^{*}:=p(x)DD. \\
\end{array}\]
The first reduced system is therefore:
\begin{eqnarray}
u_{1,1}=1 \label{reduced} \\
u_{1,2}=x \nonumber
\end{eqnarray}
and the second one consists only of unity.
We shall further need dual systems of functions:
\begin{eqnarray}
u_{1}^{*}=1 \nonumber\\
u_{2}^{*}(y)=y \label{dual} \\ 
u_{3}^{*}(y)=\int_{0}^{y}dt_{3} 
\int_{0}^{t_{3}}\frac{dt_{2}}{p(t_{2})} 
\nonumber   
\end{eqnarray} 
and  the  reduced  dual  systems, obtained  by  applying  $L^{*}$ 
operators to system (~\ref{dual}):
\begin{eqnarray}
u_{1,1}^{*}=1  \\
u_{1,2}^{*}(y)=\int_{0}^{y} \frac{dt_{3}}{p(t_{3})} \nonumber
\end{eqnarray}
\[u_{2,1}^{*}=1.\]
Let  $d\sigma:=\frac{1}{p(t)}dt$, and let  us  suppose 
that all integrals in (~\ref{primal}) exist. Then the following 
Lemma holds true:
\begin{lemma}
If $p(x) \geq 0$ then functions in system (~\ref{primal}) form  a 
{\em  Canonical  Complete Tchebycheff system\/}  (CCT-system)  on 
$(0,1).$ 
\label{l1}
\end{lemma}
{\bf  Proof. \/}Since  $\sigma$  is  monotone  by  assumption, the 
conclusion follows form \cite{Schu:1}.$\bullet$\\
\paragraph{}
{\bf  Note. \/}We must be cautious  not to employ  anything  from 
the  Tchebycheff system apparatus involving differentiability  of 
the functions in (~\ref{primal}). \\ 

\mbox{}

\section{Tchebycheff Splines}
We  may  now construct the Tchebycheff splines as it is  done  in 
\cite{Schu:1,Schu:2}.  Let us  define  two  important 
Green' s functions:
\begin{eqnarray*}
g_{3}(x,y)=\int_{y}^{x} \frac{dt_{2}}{p(t_{2})} \int_{y}^{t_{2}}dt_{3} 
            &x>y\\
g_{3}(x,y)=0 & x\leq y
\end{eqnarray*}
and
\begin{eqnarray*}
g_{2}(x,y)=\int_{y}^{x}dt_{3} & x>y\\
g_{2}(x,y)=0 & x\leq y.
\end{eqnarray*}
The first one is a Green's function of the system  (~\ref{primal}) 
and the second one of the system (~\ref{reduced}).\\
{\bf  Note. \/}We have departed here from the standard  definition 
by  which  $g_{2}(x,y)$  is  a  Green's  function  in  CCT-system 
consisting of the first two functions in (~\ref{primal}).\\
\begin{lemma}
\[L_{1}g_{3}(x,y)=g_{2}(x,y)\]  where the derivative is taken  with 
respect to the first variable.
\label{l2}
\end{lemma}
{\bf  Proof. \/}If $x>y$ we have by \cite[Lemma 4.1]{Schu:1}:
\[L_{1}g_{3}(x,y)=L_{1}\sum_{i=1}^{i=3}(-1)^{3-i}u_{i}(x)u_{4-i}^{*}(y)
=\sum_{i=2}^{i=3}(-1)^{3-i}u_{1,i-1}(x)u_{4-i}^{*}(y).\]
The  term on the right-hand side we may recognize as  $g_{2}(x,y)$ 
if    we    apply    the   quoted    Lemma    to    the    system 
(~\ref{reduced}).$\bullet$ \\
\paragraph{}
In  what  follows we shall compute the Tchebycheff  (TB-)  splines  as 
certain weighted integrals of polynomial splines. To distinguish
between  them, we  shall denote the polynomial B-splines  of  order  k 
vanishing   at   $t_{i}$   and   $t_{i+k}$   as   $ {\displaystyle 
B_{i}^{k}}$   and  the  TB-spline  as ${\displaystyle TB_{i}^{k}}$.
The following Theorem is based on Lemma (~\ref{l1}):
\begin{theorem}
\begin{equation}
TB_{i}^{3}(x)=\frac{1}{C(i)} \int_{t_{i}}^{x}B_{i}^{2}d\sigma - 
              \frac{1}{C(i+1)}\int_{t_{i+1}}^{x}B_{i+1}^{2}d\sigma  
\label{fundament}
\end{equation}
where $C(i)$ is defined as: \\
\begin{equation}
C(i):=(t_{i+2}-t_{i})\Delta^2 u_{3}^{*}(t_{i},t_{i+1},t_{i+2})
\label{C(i)}
\end{equation}
${\displaystyle   \Delta^{2}}$   being   the   ordinary   divided 
difference.
\end{theorem}
{\bf Proof. \/}By standard definition a normalized TB-spline is  a 
generalized  divided difference of Green's function  $g_{3}(x,y)$ 
(see   \cite[page  382]{Schu:2}).  By   Sylvester's   determinant 
identity (see \cite{Lyche}) we have:\\
\begin{equation}
TB_{i}^{3}(x)=(-1)^{3}(<<t_{i+1},t_{i+2},t_{i+3}>>g_{3}(x,\cdot)-
              <<t_{i},t_{i+1},t_{i+2}>>g_{3}(x,\cdot)
\label{first}
\end{equation}
where  the differences $<< >>$ are with respect to the dual  system 
(~\ref{dual}). We  can  now  reduce the  order  of  difference  by 
applying M\"{u}lbach's *recursion \cite{Mulbach} to each member of 
(~\ref{first}). The first part of the formula is:
\begin{equation}
\frac{<t_{i+2},t_{i+3}>g_{3}(x,\cdot)-
<t_{i+1},t_{i+2}>g_{3}(x,\cdot)}
{<t_{i+2},t_{i+3}>u_{3}^{*}-
<t_{i+1},t_{i+2}>u_{3}^{*}}
\label{second}
\end{equation}
and the second one is obtained by shifting index $i$. Differences 
are now with respect to the first two functions of (~\ref{dual}). 
But  since  these are $ 1$ and $y$, the differences are up  to  a 
constant  multiplier the standard polynomial ones.  By  *recursion 
for ordinary divided differences (~\ref{second}) is equal to:
\begin{equation}
\frac{(t_{i+3}-t_{i+1)}\Delta^{2}g_{3}(x;t_{i+1},t_{i+2},t_{i+3})}
{(t_{i+3}-t_{i+1)}\Delta^{2}u_{3}^{*}(t_{i+1},t_{i+2},t_{i+3})}.
\label{third}
\end{equation}
We may  apply  the 
$L_{1}$ operator to $TB_{i}^{3}$. By Lemma (~\ref{l2}) we  obtain 
polynomial differences of $g_{2}(x,\cdot)$ on the right-hand side. 
But  $g_{2}$ is a well known function $(x-\cdot)_{+}^{1}$. Thus by 
Definition  of a linear B-spline and (~\ref{C(i)})   we  conclude 
that:
\begin{equation}
L_{1}TB_{i}^{3}=
\frac{1}{C(i)}B_{i}^{2}(x)- 
              \frac{1}{C(i+1)}B_{i+1}^{2}(x).
\label{L1}
\end{equation}
The   conclusion  of  the  Theorem  follows  on  integration   of 
(~\ref{L1}) with respect to $d\sigma$.$\bullet$ \\
\paragraph{}
The  following Lemma may be used to calculate (~\ref{C(i)}) in  a 
stable way:
\begin{lemma}
\begin{equation}
C(i):=\int_{t_{i}}^{t_{i+2}}B_{i}^{2}d\sigma
\label{C1(i)}
\end{equation}
\end{lemma}
{\bf  Proof. \/}We  may explicitly calculate  all  differences  in 
(~\ref{C(i)}). An easier proof follows from Peano representation
Theorem \cite[Theorem 4.23]{Schu:1}.$\bullet$ \\
\paragraph{}
At  the  end,  we summarize  all  the  preceding  facts  in  the 
following\\
\begin{theorem}
\[TB_{i}^{3}(x)=\frac{1}{C(i)} \int_{t_{i}}^{x}B_{i}^{2}d\sigma - 
              \frac{1}{C(i+1)}\int_{t_{i+1}}^{x}B_{i+1}^{2}d\sigma\]
\[L_{1}TB_{i}^{3}(x)=\frac{1}{C(i)}B_{i}^{2}(x)- 
              \frac{1}{C(i+1)}B_{i+1}^{2}(x).\]
\begin{eqnarray*}
L_{2}TB_{i}^{3}(x)=\frac{B_{i}^{1}(x)}{C(i)(t_{i+1}-t_{i})}-
                 \frac{B_{i+1}^{1}(x)}{t_{i+2}-t_{i+1}}
                 (\frac{1}{C(i+1)}+\frac{1}{C(i)})+ \\
                 \frac{B_{i+2}^{1}(x)}{C(i+1)(t_{i+3}-t_{i+2})}
\end{eqnarray*}
\label{final}
\end{theorem}   
{\bf  Proof. \/}We  have only to prove the  last  equation.  Since 
${\displaystyle L_{2}=Dp(x)D}$, we just take a derivative of  the 
second relation and then apply a recurrence for derivatives of B-
splines \cite[page 138]{deBoor:2}$\bullet$ \\

\mbox{}

\section{Collocation}
Equation   (~\ref{eq:difequ})   may  be  written  in   the   form 
${\displaystyle  L_{2}u(x)=f(x,u(x))}$.  We  shall  assume   that 
boundary  conditions  are certain linear functional  of  function 
values   and   their  ${\displaystyle  L_{1}-   }$   derivatives, 
independent   over  $Ker(Dp(x)Du)$.  The   collocation   solution  
$u_{c}$  to  (~\ref{eq:difequ}) is a linear  combination  of  TB-
splines:
\[u_{c}(x)= \sum_{i=1}^{m-1}c_{i}TB_{i}^{3}(x)\]
if we assume the knot sequence (~\ref{knots}).\\
{\bf Note. \/} Interior knots are all of multiplicity 1,  implying 
therefore  ${\displaystyle  L_{1}-}$  differentiability  of   the 
spline space. For a singular coefficient the splines need not  be 
"smooth" in the classical sense.\\
There   is   only  one  collocation  point  per   interval,   say 
$x_{j}$, the  total of $m-3$ points. The collocation matrix is:
\[L_{2}TB_{i}^{3}(x_{j}).\]
Since TB-splines have a support of three intervals, the matrix is 
three  diagonal,  if  the  missing  two  boundary  equations  are 
carefully  placed.  The  nonlinear  problem  may  be  solved   by 
quadratically convergent Newton's method (given sufficiently good 
initial   guess),  in  case  of  an  isolated  solution,  as   in 
\cite{deBoor:1,Bellman}.  The convergence in linear case  may  be 
proved  to  be  ${\displaystyle \Delta  d\sigma  (\Delta)}$.  See 
\cite{Reddien:2} for an idea of proof. 
\paragraph{}
The placement of collocations points is of interest. We have used 
midpoints  of subintervals as collocation points. There seems  to 
be no better choice for collocation with splines of order  three; 
the similar situation arising in the parabolic spline collocation 
for  a second order differential equation. The situation  changes 
with higher order splines.
Let  us  consider the model example (~\ref{eq:model}).  It  is  a 
special  case  of a hypergeometric equation, and  may  be  easily 
solved  via functions of parabolic cylinder. The solution  has  a 
large endslope near zero, so polynomial splines are inapplicable. 
The  collocation matrix is a sum of two matrices: the  matrix  of 
${\displaystyle  L_{2}}$ derivatives and the matrix  of  function 
values  of TB-splines. The former is easily evaluated by  Theorem     
~\ref{final}.  However,   we  must calculate  integrals  for  TB-
splines  itself. This is done in the easiest way by  some  simple 
Gaussian formula, such as:
\begin{equation}
\int_{a}^{b}    \frac{f(x)}{\sqrt{x}}dx   \simeq    2(\sqrt{b}-
\sqrt{a})f(\frac{b+\sqrt{ab}+a}{3})
\label{gauss}
\end{equation}
the  formula being exact on linear functions. By  (~\ref{gauss}) 
we  have constructed an algorithm similar to deBoor's  BSPLVD  in 
\cite[page 287]{deBoor:2}. Since TB-spline values are obtained by 
subtraction  there  is  a question  of  numerical  stability.  By 
(~\ref{C1(i)})  both  terms in (~\ref{fundament}) are  less  then 
unity,  and we have found by numerical experiment that  they  are 
never  comparable, regardless of the knots, as long as $p(x)$  is 
not too oscillatory. 
\paragraph{}
At the end, let us briefly compare some results in  approximating 
the model example (~\ref{eq:model}).
Let  ${\displaystyle  f_{n}(x)}$  denote  various  approximation  as 
follows. Then at $x=0.5$ we have:


\begin{tabbing}
$\Delta_{n}$\hspace{0.5in}                    \=                    Jamet\cite{Jamet}\hspace{0.25in} 
\=Reddien\cite{Reddien:2}\hspace{0.25in}   \=$TB_{3}$\\\\
$1/8$         \>0.29308          \>0.25305                \>0.25206\\
$1/16$        \>0.27826          \>0.25223                \>0.252054\\
$1/8192$      \>0.25310\\
\end{tabbing}

Analytic  solution to (~\ref{eq:model}) at $x=0.5$  is  0.252042. 
The error behaves in the same way throughout the interval , and has two sign changes,  due to the  variation  diminishing 
property of TB-splines and monotonicity of the solution.

\mbox{}

\begin{thebibliography}{7}
\bibitem{Bellman}R. E. Bellman  and R. E. Kalaba:  Quasilinearisation 
and  Nonlinear  Boundary-Value Problems, American  Elsevier  Publ. 
Co., N.Y., (1965)
\bibitem{deBoor:2}  C. deBoor:   A  Practical  Guide  to  Splines, 
Springer-Verlag, N.Y., (1978)
\bibitem{deBoor:1}C. deBoor, B. Swartz: Collocation  at  Gaussian 
points, SIAM J. Numer. Anal. Vol. 10 No 4., 582-606 (1973)
\bibitem{Varga}  P. G. Ciarlet,  F. Natterer,  R. S. Varga:  Numerical 
Methods  of  High-Order  Accuracy  for  Singular  Boundary  Value 
Problem, Numer. Math. 15, 87-99 (1970.)
\bibitem{Jamet}P. Jamet:  On  the  Convergence  of   Finite-Difference 
 approximations to One-Dimensional  Singular  Boundary 
value Problems, Numer. Math. 14, 87-99 (1970)
\bibitem{Jerome}  J.  Jerome: On Spline functions  Determined  by 
Singular  Self-Adjoint  Differential  Operators,  J.  of  Approx. 
Theory 5, 15-40 (1972) 
\bibitem{Lyche}T. Lyche: A Recurrence Relation for Chebyshevian B-splines, Constr. Approx. 155-173, (1985) 
\bibitem{Mulbach}G.   M\"{u}lbach:  A  Recurrence   Formula   for 
Generalized  Divided  Differences  and  Some  Applications,  J. of 
Approx. Theory 9, 165-172 (1973)
\bibitem{Parter} S. V. Parter: Numerical Methods for  Generalized 
Axially Symmetric potentials, J. SIAM ser. B Numer. Anal. 2, 500-516 
(1965)
\bibitem{Reddien:1} G. W. Reddien: Projection Methods and Singular  Two 
Point Boundary Value Problems, Numer. Math. 21, 193-205 (1973)
\bibitem{Reddien:2}  G. W. Reddien, L. L. Schumaker: On a  Collocation 
Method   for   Singular  Two  Point  Boundary   Value   Problems, 
Numer. Math. 25, 427-432 (1976)
\bibitem{Schu:1}   L. L. Schumaker:   On   Tchebycheffian    Spline 
Functions, J. of Approx. Theory 18, 278-303, (1976)
\bibitem{Schu:2}L. L.    Schumaker:   Spline   Functions:    Basic 
Theory, John Wiley and Sons, N.Y., (1981)
\end{thebibliography}
\end{document}
