Homec4science
Diffusion SPClibs (feature/refactoring) Browse

spclibs/docs/manual/bsplines.tex
03540012fedefeature/refactoring

Last Change
Raw File

bsplines.tex
View Options

%
% @file bsplines.tex
%
% @brief
%
% @copyright
% Copyright (©) 2021 EPFL (Ecole Polytechnique Fédérale de Lausanne)
% SPC (Swiss Plasma Center)
%
% spclibs is free software: you can redistribute it and/or modify it under
% the terms of the GNU Lesser General Public License as published by the Free
% Software Foundation, either version 3 of the License, or (at your option)
% any later version.
%
% spclibs is distributed in the hope that it will be useful, but WITHOUT ANY
% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
%
% You should have received a copy of the GNU Lesser General Public License
% along with this program. If not, see <https://www.gnu.org/licenses/>.
%
% @authors
% (in alphabetical order)
% @author Trach-Minh Tran <trach-minh.tran@epfl.ch>
% @author Stephan Brunner <stephan.brunner@epfl.ch>
%
\documentclass[a4paper]{article}
\usepackage{linuxdoc-sgml}
%\usepackage{a4wide}
\usepackage{graphicx}
\usepackage{hyperref}
\usepackage{amsmath}
%\usepackage{fancybox}
%\usepackage[notref]{showkeys}
\title{\tt BSPLINES Reference Guide}
\author{Trach-Minh Tran, Stephan Brunner, Kurt Appert}
\date{v0.3, February 2012}
\abstract{Generalized splines of any order on irregular grids for
interpolation and solving PDEs with FEM.}
\begin{document}
\maketitle
\tableofcontents
\section{Properties of Splines}
In this section, several properties of splines will be shown in
a more or less rigorous way. The aim is mainly to provide a
minimum mathematical background for using the module \texttt{BSPLINES} for the
\emph{interpolation} problem as well as the \emph{Finite Element Method}
to solve PDEs. More rigorous mathematical proofs can be found
in the book by de Boor~\cite{deBoor}.
\subsection{Recurrence Relation}
We start by defining a finite interval $[a,b]$ subdivided into $N_x$
intervals:
\begin{equation}
a=t_0 \le t_1 \le \ldots \le t_{N_x}=b.
\end{equation}
The sequence $t_i, i=0,\ldots,N_x$ can be irregularly spaced. The $j^{th}$
spline of degree $p$ defined on this sequence of grid points (also called
\textbf{knots}), is denoted by $\Lambda_{j}^{p}$ and can be constructed using
the following recurrence relation. Starting with the \emph{constant} spline
\begin{equation}
\Lambda_i^0(x) =
\begin{cases}
1& \text{if $t_i \le x < t_{i+1}$}, \\
0& \text{otherwise}.
\end{cases}
\end{equation}
the splines of degree $p>0$ for $t_i \le x < t_{i+1}$ can be constructed from
\begin{eqnarray}
\Lambda_i^p &=& w_i^p\Lambda_i^{p-1} + (1-
w_{i+1}^p)\Lambda_{i+1}^{p-1}, \label{eq:recRel}\\
w_{i}^{p} &=& \frac{x-t_i}{t_{i+p}-t_i}.
\end{eqnarray}
Thus the values of all \emph{non-zero}
splines up to degree $p$ in the interval $[t_i,t_{i+1}]$ fit into
the triangular array as shown in Fig.~\ref{fig:allSpl}. Starting from the
first column with $\Lambda_i^0=1$, one can compute each of the $p+1$ entries
in a subsequent column with Eq.~(\ref{eq:recRel}). Applying this procedure
to generate splines on every intervals $[t_i,t_{i+1}], i=0,\ldots,N_{x}-1$
would produce the sequence of $N_{x}+p$ splines of degree $p$:
\( \Lambda_{-p}^{p},\ldots, \Lambda_{N_x-1}^{p}\).
\subsection{Support and positivity}
The linear spline
\begin{equation*}
\Lambda_i^1 = w_i^1\Lambda_i^0 + (1-w_{i+1}^1)\Lambda_{i+1}^{0}
=\frac{x-t_i}{t_{i+1}-t_i}\Lambda_i^0 +
\frac{t_{i+2}-x}{t_{i+2}-t_{i+1}}\Lambda_{i+1}^0
\end{equation*}
consists of 2 \emph{linear pieces} on $[t_i,t_{i+2}]$, forming a $C^0$
function which breaks at $t_{i+1}$ and vanishes outside of this
interval. Likewise, the quadratic spline
\begin{eqnarray*}
\Lambda_i^2 &=& w_i^2\Lambda_i^1 + (1-w_{i+1}^2)\Lambda_{i+1}^{1} \\
&=& w^2_i w^1_i \Lambda_i^0 + [w^2_i(1-w^1_{i+1})+w^1_{i+1}(1-w^2_{i+1})]
\Lambda_{i+1}^0 + (1-w^2_{i+1})(1-w^1_{i+2})\Lambda_{i+2}^0
\end{eqnarray*}
consists of 3 \emph{parabolic pieces} on $[t_i,t_{i+3}]$ that
join to form a $C^1$ function which breaks at $t_{i+1}$ and $t_{i+2}$
and vanishes outside of this interval. In general the spline of
degree $p$ can be expressed as:
\begin{equation}
\Lambda^p_i = \sum_{r=0}^{p}\, b_{i+r}^p\Lambda_{i+r}^0
\end{equation}
where $b_{i+r}^p$ is a sum of products of $p$ linear functions, resulting
in $p+1$ polynomials of degree $p$, joining to form
a $C^{p-1}$ function which breaks at $t_i,\ldots,t_{i+p+1}$
and vanishes outside of the \emph{support} $[t_i,t_{i+p+1}]$. From its construction,
$\Lambda^p_i$ is clearly \emph{strictly positive} on the interior of $[t_i,t_{i+p+1}]$.
\begin{equation}
\Lambda^p_i (x) > 0, \qquad t_i<x<t_{i+p+1}.
\end{equation}
\begin{figure}
\centering
\( \begin{array}{cccccc}
& & & & & \\
& & & & & 0 \\
& & & & 0 & \\
& & & \cdot& & \Lambda^p_{i-p} \\
& & 0 & & \Lambda^{p-1}_{i-p+1} & \\
& 0 & & \cdot& & \Lambda^p_{i-p+1} \\
0 & &\Lambda^{2}_{i-2}& &\Lambda^{p-1}_{i-p+2} & \\
& \Lambda^{1}_{i-1}& & \cdot& & \Lambda^p_{i-p+2} \\
\Lambda^{0}_{i} & &\Lambda^{2}_{i-1}& & \cdot & \\
& \Lambda^{1}_{i} & & \cdot& & \cdot \\
0 & &\Lambda^{2}_{i} & & \Lambda^{p-1}_{i-1} & \\
& 0 & & \cdot& & \Lambda^p_{i-1} \\
& & 0 & & \Lambda^{p-1}_{i} & \\
& & & \cdot& & \Lambda^p_{i} \\
& & & & & \\
& & & & 0 & \\
& & & & & 0
\end{array} \)
\caption{The array of all the splines of degree up to $p$ that are non-zero in
$[t_i,t_{i+1}]$.}
\label{fig:allSpl}
\end{figure}
\subsection{Sum of Splines}
For $t_i\leq x < t_{i+1}$
\begin{eqnarray*}
\sum_j\,\Lambda_j^0 &=& \Lambda_i^0 = 1, \\
\sum_j\,\Lambda_j^1 &=& \Lambda_{i-1}^1 + \Lambda_{i}^1= (1-w^1_i)\Lambda_i^0
+ w^1_i\Lambda_i^0 = 1.
\end{eqnarray*}
Thus assuming that for $p>1$:
\[ \sum_{j=i-p+1}^{i}\,\Lambda_j^{p-1} = 1, \]
or that the sum of the next to last column in Fig.~\ref{fig:allSpl} is $1$,
we have, using the recurrence relation (\ref{eq:recRel})
\begin{eqnarray*}
\sum_{j=i-p}^{i}\,\Lambda_j^{p} &=& \sum_{j=i-p}^{i}\,
\left( w^p_j\Lambda_j^{p-1} +(1-w^p_{j+1})\Lambda_{j+1}^{p-1} \right) \\
&=& \sum_{j=i-p+1}^{i}\,w^p_j\Lambda_j^{p-1} +
\sum_{j=i-p+1}^{i}\,(1-w^p_j)\Lambda_j^{p-1} \\
&=& \sum_{j=i-p+1}^{i}\,\Lambda_j^{p-1} = 1.
\end{eqnarray*}
\subsection{Derivative of Splines}
The derivative of the splines of degree $p$ can be expressed in terms of the
splines of degree $p-1$ by the following relation:
\begin{equation}
\label{derivative of splines}
\frac{d}{dx}\Lambda_i^p =
p\left(
\frac{\Lambda_i^{p-1}}{t_{i+p}-t_i}
- \frac{\Lambda_{i+1}^{p-1}}{t_{i+p+1}-t_{i+1}}
\right).
\end{equation}
A straightforward consequence of this relation is that the splines of order
$p$ are $C^{p-1}$ continuous. The demonstration of Eq.(\ref{derivative of
splines}) is done by induction. One starts with the case $p=1$:
\begin{eqnarray*}
\frac{d}{dx}\Lambda_i^1
& = &
\frac{d}{dx}\left[
w_i^1\Lambda_i^0 + (1-w^1_{i+1})\Lambda_{i+1}^0
\right]
\\
& = &
\frac{d\,w_i^1}{dx}\Lambda_i^0 + \frac{d\,(1-w^1_{i+1})}{dx}\Lambda_{i+1}^0
+ w_i^1\frac{d\,\Lambda_i^0}{dx} + (1-w^1_{i+1})\frac{d\,\Lambda_{i+1}^0}{dx}
\\
& = &
\frac{1}{t_{i+1}-t_i}\Lambda_i^0 - \frac{1}{t_{i+2}-t_{i+1}}\Lambda_{i+1}^0,
\end{eqnarray*}
having used Eq.(\ref{eq:recRel}) and $d\,\Lambda_i^0/dx = 0$. One then assumes
Eq.(\ref{derivative of splines}) true for $p-1$ and demonstrates that it
remains true for $p$. This is done as follows:
\begin{eqnarray}
\label{demo deriv. 1}
\frac{d}{dx}\Lambda_i^p
& = &
\frac{d}{dx}\left[
w_i^p\Lambda_i^{p-1} + (1-w^p_{i+1})\Lambda_{i+1}^{p-1}
\right]
\\
\nonumber
& = &
\frac{d\,w_i^p}{dx}\Lambda_i^{p-1}
+ \frac{d\,(1-w^p_{i+1})}{dx}\Lambda_{i+1}^{p-1}
+ w_i^p\frac{d\,\Lambda_i^{p-1}}{dx}
+ (1-w^p_{i+1})\frac{d\,\Lambda_{i+1}^{p-1}}{dx}
\\
\nonumber
& = &
\frac{\Lambda_i^{p-1}}{t_{i+p}-t_i}
-\frac{\Lambda_{i+1}^{p-1}}{t_{i+p+1}-t_{i+1}}
\\
\label{demo deriv. 2}
&&
+ w_i^p (p-1)\left(
\frac{\Lambda_i^{p-2}}{t_{i+p-1}-t_i}
- \frac{\Lambda_{i+1}^{p-2}}{t_{i+p}-t_{i+1}}
\right)
+ (1-w_{i+1}^p) (p-1)\left(
\frac{\Lambda_{i+1}^{p-2}}{t_{i+p}-t_{i+1}}
- \frac{\Lambda_{i+2}^{p-2}}{t_{i+p+1}-t_{i+2}}
\right)
\end{eqnarray}
having used Eq.(\ref{eq:recRel}) to obtain (\ref{demo deriv. 1}), and the
induction hypothesis to obtain Eq.(\ref{demo deriv. 2}). Now, rearranging the
last two terms of Eq.(\ref{demo deriv. 2}), one easily obtains:
\begin{eqnarray}
\nonumber
\frac{d}{dx}\Lambda_i^p
& = &
\frac{\Lambda_i^{p-1}}{t_{i+p}-t_i}
-\frac{\Lambda_{i+1}^{p-1}}{t_{i+p+1}-t_{i+1}}
\\
\nonumber
&&
+(p-1)\left[
\frac{1}{t_{i+p}-t_i}
\left(
\frac{x-t_i}{t_{i+p-1}-t_i}\Lambda_i^{p-2}
+\frac{t_{i+p}-x}{t_{i+p}-t_{i+1}}\Lambda_{i+1}^{p-2}
\right)
\right.
\\
\nonumber
&&
\hspace{2.cm}
\left.
-\frac{1}{t_{i+p+1}-t_{i+1}}
\left(
\frac{x-t_{i+1}}{t_{i+p}-t_{i+1}}\Lambda_{i+1}^{p-2}
+\frac{t_{i+p+1}-x}{t_{i+p+1}-t_{i+2}}\Lambda_{i+2}^{p-2}
\right)
\right]
\\
\label{demo deriv. 3}
& = &
\frac{\Lambda_i^{p-1}}{t_{i+p}-t_i}
-\frac{\Lambda_{i+1}^{p-1}}{t_{i+p+1}-t_{i+1}}
+ (p-1)\left(
\frac{\Lambda_i^{p-1}}{t_{i+p}-t_i}
-\frac{\Lambda_{i+1}^{p-1}}{t_{i+p+1}-t_{i+1}}
\right)
\\
\nonumber
& = &
p\left(
\frac{\Lambda_i^{p-1} }{t_{i+p}-t_i}
-\frac{\Lambda_{i+1}^{p-1}}{t_{i+p+1}-t_{i+1}}
\right)
\end{eqnarray}
having again used Eq.(\ref{eq:recRel}) to obtain (\ref{demo deriv. 3}). This
completes the demonstration of relation (\ref{derivative of splines}).
\subsection{Integrals of Splines}
With the proper normalization all splines of all degrees have unitary surface:
\begin{equation}
\label{integrals of splines}
\frac{p+1}{t_{i+p+1}-t_i}\int \Lambda_i^p(x)dx = 1.
\end{equation}
This relation holds trivially for $p=0$ and $p=1$. A recursive proof of the
general statement (\ref{integrals of splines}) starts assuming
\begin{equation}
\label{previousInt}
\frac{p}{t_{i+p}-t_i}\int \Lambda_i^{p-1}(x)dx = 1
\end{equation}
to be true. Then using Eq.(\ref{derivative of splines}) multiplied by $x$ and integrating one obtains:
\begin{equation}
\nonumber
\int x \frac{d}{dx}\Lambda_i^p dx = -\int \Lambda_i^p dx =
p\int\left( \frac{x\Lambda_i^{p-1}}{t_{i+p}-t_i}
-\frac{x\Lambda_{i+1}^{p-1}}{t_{i+p+1}-t_{i+1}}
\right)dx.
\end{equation}
Completing the fractions in the big parentheses in view of using
Eq.(\ref{eq:recRel}) one has
\begin{eqnarray}
\nonumber
\int \Lambda_i^p dx
& = &
-\, p\int \frac{x-t_i}{t_{i+p}-t_i}\Lambda_i^{p-1} dx
- p\int \frac{t_i}{t_{i+p}-t_i}\Lambda_i^{p-1} dx \\
\nonumber
& &
-\, p\int \frac{t_{i+p+1}-x}{t_{i+p+1}-t_{i+1}}\Lambda_{i+1}^{p-1} dx
+ p\int \frac{t_{i+p+1}}{t_{i+p+1}-t_{i+1}}\Lambda_{i+1}^{p-1} dx,
\end{eqnarray}
where the first and the third terms on the right correspond to $-p\int\Lambda_i^p dx $,
Eq.(\ref{eq:recRel}), and can be combined with the left side to yield
\begin{equation}
\label{proofIntegrals}
(1+p)\int\Lambda_i^p dx = t_{i+p+1}-t_i,
\end{equation}
where relation (\ref{previousInt}) has been used for the rest on the right hand side.
This concludes the proof of Eq.(\ref{integrals of splines}).
\subsection{Boundary Conditions}
Applying the recurrence relation to generate \emph{all} the splines on
the finite domain $[t_0,t_{N_x}]$ yields the $N_x+p$ splines
of degree $p$:
\begin{equation}
\Lambda^p_{-p},\Lambda^p_{-p+1},\ldots, \Lambda^p_{N_x-1}. \label{eq:splSeq}
\end{equation}
Note that \emph{additional} knots beyond both ends of $[t_0,t_{N_x}]$ have to
be defined to generate all these splines.
\subsubsection{Periodic splines}
The extra knots are simply defined through periodicity.:
\begin{eqnarray}
t_{-\nu} &=& t_{N_x-\nu}-(b-a), \\
t_{N_{x}+\nu} &=& t_{\nu}+(b-a), \qquad \nu=0,\ldots,p.
\end{eqnarray}
The $p+1$ leftmost splines in (\ref{eq:splSeq}) are thus identical to the rightmost splines:
\begin{equation}
\Lambda^p_{-\nu} = \Lambda^p_{N_x-\nu}, \qquad \nu=0,\ldots,p.
\end{equation}
\subsubsection{Non-periodic splines}
The choice made in \texttt{BSPLINES} is simply:
\begin{equation}
t_{-p} = \cdots = t_{0} = a, \qquad b=t_{N_x}=\cdots=t_{N_x+p}.
\end{equation}
Thus in the first interval $[t_0,t_1]$, the first spline $\Lambda^p_{-p}$
is constructed (refer to the first entry on each of the column of
Fig.~\ref{fig:allSpl}, with $i=0$) as follow:
\begin{eqnarray*}
\Lambda^1_{-1} &=& (1-w^1_{0})\Lambda^0_{0} = \frac{t_{1}-x}{t_{1}-t_{0}}\Lambda^0_{0}\\
\Lambda^2_{-2} &=& (1-w^2_{-1})\Lambda^1_{-1} =\frac{t_{1}-x}{t_{1}-t_{-1}}\Lambda^1_{-1}
=\left(\frac{t_{1}-x}{t_{1}-t_{0}}\right)^{2}\Lambda^0_{0}\\
\cdot & & \qquad\cdot\qquad\qquad\qquad\qquad \cdot \\
\Lambda^p_{-p} &=& (1-w^p_{-p+1})\Lambda^p_{-p+1} =\frac{t_{1}-x}{t_{1}-t_{-p+1}}\Lambda^{p-1}_{-p+1}
=\left(\frac{t_{1}-x}{t_{1}-t_{0}}\right)^{p}\Lambda^0_{0}
\end{eqnarray*}
In the same manner, the generation of the \emph{last} spline $\Lambda^p_{N_x-1}$ (last entry
on each of the column of Fig.~\ref{fig:allSpl}, with $i=N_{x}-1$) yields:
\begin{eqnarray*}
\Lambda^1_{N_x-1} &=& w^1_{N_x-1}\Lambda^0_{N_x-1} =
\frac{x-t_{N_x-1}}{t_{N_x}-t_{N_x-1}}\Lambda^0_{N_x-1} \\
\Lambda^2_{N_x-1} &=& w^2_{N_x-1}\Lambda^1_{N_x-1} =
\frac{x-t_{N_x-1}}{t_{N_x+1}-t_{N_x-1}}\Lambda^1_{N_x-1}
= \left(\frac{x-t_{N_x-1}}{t_{N_x}-t_{N_x-1}}\right)^{2}\Lambda^0_{N_x-1}\\
\cdot & & \qquad\cdot\qquad\qquad\qquad\qquad \cdot \\
\Lambda^p_{N_x-1} &=& w^p_{N_x-1}\Lambda^p_{N_x-1} =
\frac{x-t_{N_x-1}}{t_{N_x+p-1}-t_{N_x-1}}\Lambda^1_{N_x-1}
= \left(\frac{x-t_{N_x-1}}{t_{N_x}-t_{N_x-1}}\right)^{p}\Lambda^0_{N_x-1}
\end{eqnarray*}
Since the sum of all splines is 1 and using the positivity of splines, all the
non-periodic splines, except the first (last) spline should vanish at $x=a$ ($x=b$):
\begin{equation}
\Lambda^p_{r}(a) = \delta_{r,-p}, \qquad\Lambda^p_{r}(b) = \delta_{r,N_x-1}
\end{equation}
The spline derivatives at the boundaries $x=a$ and $x=b$ can be
derived using Eq.(\ref{derivative of splines}) as follow. At $x=a$
(interval $[t_0,t_1]$), by noting that only the spline
$\Lambda^{p-1}_{-p+1}$ is non-zero at $x=a$ (see next to last column
of Fig.{\ref{fig:allSpl}, with $i=0$), it is easy to see that
there are only 2 non-zero derivatives
given by
\begin{equation}
\begin{split}
\frac{d}{dx}\Lambda^p_{-p}(a) =&
-\frac{p\,\Lambda^{p-1}_{-p+1}(a)}{t_1-t_{-p}} = -\frac{p}{t_1-t_0}, \\
\frac{d}{dx}\Lambda^p_{-p+1}(a) =&
\frac{p\,\Lambda^{p-1}_{-p+1}(a)}{t_1-t_{-p+1}} = \frac{p}{t_1-t_0},
\end{split}
\end{equation}
where we have used $t_0=t_{-1}=\ldots=t_{-p}=a$. Likewise, the 2
non-zero derivatives of spline at the other boundary $x=b$ are
\begin{equation}
\begin{split}
\frac{d}{dx}\Lambda^p_{N_x-2}(b) =&
-\frac{p\,\Lambda^{p-1}_{N_x-1}(b)}{t_{N_x+p-1}-t_{N_x-1}} = -\frac{p}{t_{N_x}-t_{N_x-1}}, \\
\frac{d}{dx}\Lambda^p_{N_x-1}(b) =&
\frac{p\,\Lambda^{p-1}_{N_x-1}(b)}{t_{N_x+p-1}-t_{N_x-1}} = \frac{p}{t_{N_x}-t_{N_x-1}},
\end{split}
\end{equation}
where we have used $t_{N_x}=t_{N_x+1}=\ldots=t_{N_x+p}=b$.
\subsubsection{Spline expansion}
In summary, the approximation of a function $f$ defined
in the interval $[a,b]$ using a basis (\textbf{Is this obvious?) }of
splines of degree $p$ associated with the sequence of knots
$t_i, i=-p,\ldots,N_{x}+p$
can be written as
\begin{equation}
f(x) = \sum_{j=-p}^{N_x-1}\, c_j\Lambda^p_j(x), \qquad
\begin{array}{l}
\mbox{support of $\Lambda^p_j$:}\quad [t_{j},t_{j+p+1}],\\
t_i \leq x < t_{i+1} \Longrightarrow \Lambda^p_{i-p}(x),\ldots,
\Lambda^p_{i}(x) \ge 0.
\end{array}
\end{equation}
Note that the \emph{last} spline in the interval $[t_i,t_{i+1}]$, which can be written as
\[ \Lambda^p_{i}(x)=w^p_i(x) \Lambda^{p-1}_{i}(x)=\ldots=w^p_i(x) w^{p-1}_i(x) \ldots
w^{1}_i(x)\Lambda^{0}_{i}(x) \]
\emph{vanishes at the knot} $x=t_i$. Thus at any position $x$, the sum
involves $p+1$ terms except at the knots $t_i$ where there are only $p$ terms.
It is sometimes more convenient to renumber the spline index $j$ so that it
starts from $0$. With this new numbering, the spline expansion becomes
\begin{equation}
f(x) = \sum_{j=0}^{N_x+p-1}\, c_j\Lambda^p_j(x), \qquad
\begin{array}{l}
\mbox{support of $\Lambda^p_j$:}\quad [t_{j-p},t_{j+1}], \\
t_i \leq x < t_{i+1} \Longrightarrow \Lambda^p_{i}(x),\ldots,
\Lambda^p_{i+p}(x) \ge 0.
\end{array}
\label{eq:splExp}
\end{equation}
In the \emph{periodic} case, there are $N_{x}$ \emph{independent}
spline coefficients since
\begin{equation}
c_{N_{x}+\nu} = c_{\nu}, \qquad \nu=0,\ldots,p-1. \label{eq:perSp}
\end{equation}
In the \emph{non-periodic} case,
the first and the last spline coefficients $c_{0},\,c_{N_x+p-1}$ are
respectively the values of $f$ at the end points $a$ and $b$.
The basis functions for both non-periodic and periodic cubic splines
($p=3$) are shown in Fig~.\ref{fig:cubic_splines} where this new numbering is
used.
\begin{figure}[htbp]
\centering
\includegraphics[angle=0,width=\hsize]{driv1}
\caption{The basis of non-periodic and periodic cubic splines. The periodic
splines $\Lambda_{10}$, $\Lambda_{11}$, $\Lambda_{12}$ denote the same
splines as $\Lambda_{0}$, $\Lambda_{1}$, $\Lambda_{2}$ respectively. }
\label{fig:cubic_splines}
\end{figure}
\subsection{Spline Initialization with \texttt{SET\_SPLINE}}
The initialization of a spline is performed by calling the routine
\texttt{SET\_SPLINE}, passing the desired degree $p$ and the
sequence of grid points (or knots) $t_j, j=0,\ldots,N_x$. If
Gauss points on each of the intervals $[t_j,t_{j+1}]$ are needed, a non-zero
value of \texttt{NGAUSS} should be specified. The other input
argument is the \emph{optional} \texttt{LOGICAL} argument \texttt{PERIOD}
to define the periodicity of the splines. By default it is \texttt{.FALSE.}.
The routine returns the 1d spline \texttt{SP} which is of type
\texttt{TYPE(spline1d)}:
\par
\addvspace{\medskipamount}
\nopagebreak\hrule
\begin{verbatim}
SUBROUTINE set_spline(p, ngauss, grid, sp, period)
INTEGER, INTENT(in) :: p, ngauss
DOUBLE PRECISION, INTENT(in) :: grid(:)
LOGICAL, OPTIONAL, INTENT(in) :: period
TYPE(spline1d), INTENT(out) :: sp
LOGICAL, OPTIONAL, INTENT(in) :: period
\end{verbatim}
\nopagebreak\hrule
\addvspace{\medskipamount}
Besides the main characteristics of the spline (degree $p$ of splines,
number of grid intervals, dimension of splines $N_x+p$, etc.)
the following quantities will be determined and stored in \texttt{SP}:
\begin{itemize}
\item values and all the $p$ derivatives of the $p+1$ non-vanishing splines
on each knots $t_j$. These quantities will be used to speed up computation
of the spline expansion (\ref{eq:splExp}).
\item integrals of splines \(I_i=\int\Lambda_i(x)\,dx\).
\end{itemize}
For a 2d spline
\begin{equation}
\Lambda^{p+q}_{ij}(x,y) = \Lambda^p_i(x)\Lambda^q_j(y),
\end{equation}
on a 2d structured mesh defined by the grid points
\texttt{grid1(0:N1), grid2(0:N2)}, the same call as in the 1d case can be
used, except that the scalars \texttt{p, ngauss, period} become 2 element arrays and
the output \texttt{SP} is now of type \texttt{TYPE(spline2d)}:
\par
\addvspace{\medskipamount}
\nopagebreak\hrule
\begin{verbatim}
INTEGER :: p(2), ngauss(2)
LOGICAL, OPTIONAL :: period(2)
DOUBLE PRECISION, dimension(:) :: grid1, grid2
TYPE(spline2d) :: sp2d
...
CALL set_spline(p, ngauss, grid1, grid2, sp2d, period)
\end{verbatim}
\nopagebreak\hrule
\addvspace{\medskipamount}
The derived type \texttt{spline2d} is a \emph{wrapper} of 2 \texttt{spline1d}
objects which can be accessed through \texttt{sp2d\%sp1} and \texttt{sp2d\%sp2}.
Once \texttt{SET\_SPLINE} is called, the routine \texttt{GET\_DIM} can be
called to inquire the spline's essential characteristics such as dimension,
number of intervals and degree, for both 1d and 2d splines: \par
\addvspace{\medskipamount}
\nopagebreak\hrule
\begin{verbatim}
SUBROUTINE get_dim(sp, dim, nx, nidbas)
TYPE(spline1d), INTENT(in) :: sp
INTEGER, INTENT(out) :: dim
INTEGER, OPTIONAL, INTENT(out) :: nx, nidbas
\end{verbatim}
\nopagebreak\hrule
\addvspace{\medskipamount}
Integral of function $\int^b_a\,f(x)dx$ is computed from its spline
\texttt{sp} and splines coefficients in:
\par
\addvspace{\medskipamount}
\nopagebreak\hrule
\begin{verbatim}
DOUBLE PRECISION FUNCTION fintg(sp, c)
TYPE(spline1d), INTENT(in) :: sp
DOUBLE PRECISION, INTENT(in) :: c(:)
\end{verbatim}
\nopagebreak\hrule
\addvspace{\medskipamount}
For a 2d functions, the same function should be called with a 2d spline
\texttt{sp} and 2d array $c$.
Finally \texttt{DESTROY\_SP(sp)} should be called when a spline
\texttt{sp} is not needed anymore to clean up memory space.
\subsection{Generating Splines with \texttt{DEF\_BASFUN}}
\par
\addvspace{\medskipamount}
\nopagebreak\hrule
\begin{verbatim}
SUBROUTINE def_basfun(xp, sp, fun, left)
DOUBLE PRECISION, INTENT(in) :: xp
TYPE(spline1d), INTENT(in) :: sp
DOUBLE PRECISION, INTENT(out) :: fun(:,:)
INTEGER, OPTIONAL, INTENT(out) :: left
\end{verbatim}
\nopagebreak\hrule
\addvspace{\medskipamount}
This routine computes, for a given point $\mbox{xp}\in [t_0,t_{N_x}]$, the value
and optionally the $m$ derivatives of the $p+1$ splines \texttt{sp} which were previously
defined and returns them in \texttt{fun(1:p+1,1:m+1)} with $m\leq p$. The maximum
number of computed derivatives $m$ is determined by the size of the second dimension
of the array \texttt{fun}. The subroutine will return the optional integer
\texttt{left} defined such that:
\[
t_{\mbox{left}} \leq xp < t_{\mbox{left+1}}, \qquad 0\leq \mbox{left} \leq N_{x.-1}.
\]
\subsection{Example 1: Values and derivatives of all splines}
In this example, we first initialize a cubic spline with
the knot sequence $t_0,\ldots,t_{N_x}$ with \texttt{SET\_SPLINE}
and then call \texttt{DEF\_BASFUN} to compute its values, first and second
derivatives on the mesh points \texttt{xp(1:npts)}.
\par
\addvspace{\medskipamount}
\nopagebreak\hrule
\begin{verbatim}
USE BSPLINES
INTEGER, PARAMETER :: nx=10, npts=100
DOUBLE PRECISION :: t(0:nx), xp(npts)
DOUBLE PRECISION, ALLOCATABLE :: fxp0(:,:), fxp1(:,:), fxp2(:,:)
DOUBLE PRECISION :: fun(4,3) ! 4 cubic splines at a given xp
! plus first and second derivatives.
INTEGER :: i, dim, left
TYPE(spline1d) :: sp
!
! Define t(0:nx), xp(npts)
!
CALL set_spline(3, 0, t, sp, period=.FALSE.)
CALL get_dim(sp, dim)
ALLOCATE(fxp0(npts,0:dim-1), fxp1(npts,0:dim-1), fxp2(npts,0:dim-1)
fxp0 = 0.0
fxp1 = 0.0
fxp2 = 0.0
DO i=1,npts
CALL def_basfun(xp(i), sp, fun, left=left)
fxp0(i, left:left+3) = fun(1:4, 1) ! Value
fxp1(i, left:left+3) = fun(1:4, 2) ! 1st derivative
fxp2(i, left:left+3) = fun(1:4, 3) ! 2nd derivative
END DO
DEALLOCATE(fxp0, fxp1, fxp2)
CALL destroy_sp(sp)
\end{verbatim}
\nopagebreak\hrule
\addvspace{\medskipamount}
This code fragment will store \texttt{dim=nx+3=13} splines
and theirs first 2 derivatives in \texttt{fxp0}, \texttt{fxp1} and
\texttt{fxp2}. Change the \texttt{period} to \texttt{.TRUE.} to obtain
\emph{periodic} splines.
\section{Spline Interpolation}
Given the interval $[a,b]$ discretized into $\{x_k,\,k=0,\ldots,N_g\}$ with
$x_0=a$ and $x_{N_g}=b$, the problem of interpolating $f(x), x\in [a,b]$ with
splines of degree $p$ is to solve the following equations
for the spline coefficients $c_i$:
\begin{equation}
\sum_{i=0}^{N_x+p-1}\,c_i\Lambda^p_{i}(x_k) = f(x_k),\quad k=0,\ldots,N_g.
\label{eq:intEq}
\end{equation}
The sequence of knots $t_0,\ldots,t_{N_x}$ defines completely the splines
$\Lambda^p_i$ and its choice will be described in the following section.
\subsection{Choice of knots}
If Eqs.~(\ref{eq:intEq}) are the only conditions for our interpolation
problem, the number of equations should match the number
of unknowns $c_i$. The number of knot intervals $N_x$ hence has to verify
\begin{equation}
N_x = N_g-p+1.\label{eq:knotNum}
\end{equation}
For the \emph{periodic} case, taking into account the $p$ periodic spline conditions
(\ref{eq:perSp}) on
$c_i$ and $f(a)=f(b)$, this condition reduces to:
\begin{equation}
N_x = N_g.
\end{equation}
For \emph{odd} values of the spline degree $p$, the knots $t_i$ could be placed on
the \emph{interpolation sites} $x_k$ while when $p$ is even, $t_i$ should not
be on $x_k$ to avoid a \emph{badly conditioned} linear system when solving
Eq.~(\ref{eq:intEq}). This leads to the following choice for $t_i$ in
\texttt{BSPLINES}:
\subsubsection{Periodic splines}
The number of knots $N_x+1$ is \emph{equal} to the number of interpolation points
$N_g+1$ with
\begin{equation}
t_i =
\begin{cases}
x_i & \text{$p$ odd} \\
(x_{i-1}+x_i)/2 & \text{$p$ even}
\end{cases}
,\qquad i=0,\ldots,N_x
\end{equation}
\subsubsection{Non-Periodic splines}
In order to satisfy the equality (\ref{eq:knotNum}), first, the 2 end points
are retained as knots:
\begin{equation}
t_0=x_0, \qquad t_{N_x} = x_{N_g}.
\end{equation}
For even $p$, the first $p/2$ interpolation intervals are \emph{skipped}:
\begin{equation}
t_i = (x_{i+p/2-1} + x_{i+p/2})/2, \qquad i=1,\ldots,N_x-1, \label{eq:evenKnots}
\end{equation}
while for odd $p$, $(p-1)/2$ interpolation points are \emph{skipped}:
\begin{equation}
t_i = x_{i+(p-1)/2 }, \qquad i=1,\ldots,N_x-1. \label{eq:oddKnots}
\end{equation}
Instead of skipping grid points, an alternative would be to supplement the
system of equations (\ref{eq:intEq}) with conditions on derivatives of $f(x)$
at one or both ends of $[a,b]$. This type of boundary conditions is not
implemented in the present version of the \texttt{BSPLINES} module.
\subsection{The collocation matrix}
The \emph{collocation matrix} $\Lambda^p_i(x_k)$ of the interpolation problem
(\ref{eq:intEq}) is a square matrix. Each row has at most $p+1$ non-zero
terms. Let us consider separately the non-periodic and the periodic cases.
\subsubsection{The non-periodic case}
\paragraph{Even spline degree}
From (\ref{eq:evenKnots}), there are $p/2+1$ interpolation points \(x_{0}, \ldots, x_{p/2}\)
in the first knot interval $[t_0,t_1)$. Since there are at most $p+1$ non-zero
splines for any points in each interval (except for $x_0$ where $\Lambda_i(x_0)=\delta_{i,0}$,
the collocation matrix starts as:
\begin{equation}
\left(\begin{array}{llllll}
\Lambda_0(x_{0}) & 0 &\cdots & \cdots & \cdots &\cdots \\
\Lambda_0(x_{1}) &\Lambda_1(x_{1}) &\cdots &\Lambda_p(x_{1}) & 0 &\cdots \\
\vdots &\vdots &\cdots &\vdots & 0 &\cdots \\
\Lambda_0(x_{p/2})&\Lambda_1(x_{p/2}) &\cdots &\Lambda_p(x_{p/2}) & 0 &\cdots \\
0 &\Lambda_1(x_{p/2+1})&\cdots &\Lambda_p(x_{p/2+1})& \Lambda_{p+1}(x_{p/2+1})&0 \\
0 &\ddots &\ddots &\ddots & \ddots &\ddots
\end{array}\right)
\end{equation}
The number of \emph{upper-diagonals} (non including the diagonal) is
obviously determined by the second row of the matrix above, which yields $p-1$.
Since the knot placement is identical for both ends of the interpolation mesh,
the matrix
$\Lambda_i(x_k)$ is \emph{banded} with half-bandwidths
\begin{equation}
kl=ku=p-1
\end{equation}
\paragraph{Odd spline degree}
Applying the same procedure, it is straightforward to show for $p$ odd and
from (\ref{eq:oddKnots}), that $x_0,\ldots.x_{(p-1)/2}$ are located in the
first knot interval $[t_0,t_1)$ and that the matrix has again the same
half-bandwidths as in the even $p$ case.
The resulting interpolation problem can then be solved with
the usual \emph{banded matrix factorization} followed by a \emph{back-solve}
phase.
\subsubsection{The periodic case}
Let consider the matrix for $p=3$ and $N_x=10$ (see lower figure of
Fig.~(\ref{fig:cubic_splines}):
\begin{equation}
\left(\begin{array}{llllll}
\Lambda_0(x_{0}) & \Lambda_1(x_{0}) & \Lambda_2(x_{0}) & 0 & \cdots & \\
0 & \Lambda_1(x_{1}) & \Lambda_2(x_{1}) & \Lambda_3(x_{1}) & 0 &\cdots \\
\vdots & \ddots & \ddots & \ddots & \ddots &\ddots \\
0 & \cdots & 0 & \Lambda_7(x_{7}) & \Lambda_8(x_{7}) & \Lambda_9(x_{7}) \\
\Lambda_0(x_{8}) & 0 & 0 & \cdots &\Lambda_8(x_{8}) & \Lambda_9(x_{8}) \\
\Lambda_0(x_{9}) &\Lambda_1(x_{9}) & 0 & 0 & \cdots & \Lambda_9(x_{9})
\end{array}\right) \label{eq:perMat}
\end{equation}
The matrix is ``almost triangular'' (except for the last 2 rows) and is not
\emph{diagonally dominant}! A more satisfactory (and symmetric in shape) matrix is
however obtained by simply renumbering the splines such that the sequence
starts with $-\lfloor p/2 \rfloor$ instead of $0$. This renumbered splines
are shown in Fig.~\ref{fig:fitSpl} for the cubic
and quadratic periodic splines. With this renumbering, the matrix
(\ref{eq:perMat}) has a more symmetric shape and is diagonally dominant:
\begin{equation}
\left(\begin{array}{lllll}
\Lambda_0(x_{0}) & \Lambda_1(x_{0}) & 0 & \cdots & \Lambda_9(x_{0}) \\
\Lambda_0(x_{1}) & \Lambda_1(x_{1}) & \Lambda_2(x_{1}) & 0 &\cdots \\
\vdots & \ddots & \ddots & \ddots &\ddots \\
0 & \cdots & \Lambda_7(x_{8}) & \Lambda_8(x_{8}) & \Lambda_9(x_{8}) \\
\Lambda_0(x_{9}) & 0 & 0 & \Lambda_8(x_{9}) & \Lambda_9(x_{9})
\end{array}\right) \label{eq:perMatnew}
\end{equation}
In general, for arbitrary $p$ (even and odd values), the collocation matrix
$A=\Lambda_j(x_i)$ can be written as
\begin{equation}
A = B + UV^T
\end{equation}
where $B$ is a banded matrix with half-bandwidths $kl=ku=b=\lfloor p/2\rfloor$ and
rank $N_x$. $U$ and $V$ are $N_x\times 2b$ sparse matrices:
\begin{equation}
U = \left(
\begin{matrix}
I & 0 \\
0 & 0 \\
0 & I
\end{matrix}\right), \qquad
V = \left(
\begin{matrix}
0 & D^T \\
0 & 0 \\
C^T & 0
\end{matrix}\right), \qquad
V^T = \left(
\begin{matrix}
0 & 0 & C \\
D & 0 & 0
\end{matrix}\right), \qquad
\end{equation}
where $C$, $D$ are the $b\times b$ \emph{off-band} sub-matrices and $I$, the
identity matrix. In the cubic spline example considered above, the
\emph{off-band} matrices are simply $1\times 1$ matrices with
$C=\Lambda_9(x_0)$ and $D=\Lambda_0(x_9)$.
The inverse of $A$ can be deduced from the \emph{Sherman-Morrison-Woodbury formula} \cite{Golub}:
\begin{eqnarray*}
A^{-1} &=& B^{-1} - B^{-1}U(1+V^{T}B^{-1}U)^{-1}V^{T}B^{-1} \\
&=& B^{-1} - ZW^{T}B^{-1},
\end{eqnarray*}
where
\begin{eqnarray*}
Z &=& B^{-1}U, \\
H &=& 1+V^{T}B^{-1}U \\
W^T &=& H^{-1}V^{T}.
\end{eqnarray*}
The solution of the interpolation problem $Ax=f$ can then be reduced to a
\emph{factorization} and a \emph{back-solve} phase:
\begin{enumerate}
\item Factorization
\begin{enumerate}
\item Factor: \( B \longleftarrow L_BU_B \)
\item Solve: \( (L_BU_B)Z = U, \quad U\longleftarrow Z \)
\item Compute: \( H = 1+V^{T}Z \)
\item Factor: \( H=L_HU_H \)
\item Solve: \( (L_HU_H)W^{T} = V^{T}, \quad V^{T}\longleftarrow W^{T} \)
\end{enumerate}
\item Back-solve
\begin{enumerate}
\item Solve: \( (L_BU_B)y = f \)
\item Compute: \( t = W^{T}y \)
\item Compute: \( x = y - Zt \)
\end{enumerate}
\end{enumerate}
At the end of the factorization, only the (updated) matrices $B$, $U$ and
$V^{T}$, required in the back-solve phase, need to be saved.
Note that we avoid to store the product
$ZW^T$ because it is a \emph{big} $N_x\times N_x$ matrix.
After the \emph{back-solve} step, the solution $x$ is \emph{shifted back} (by
$\lfloor p/2\rfloor$) and the appropriate periodicity condition is applied to
obtain the spline coefficients $c_j,\, j=0,\ldots,N_x+p-1$, as defined in
(\ref{eq:splExp}).
\begin{figure}[htbp]
\centering
\includegraphics[angle=0,width=\hsize]{fit}
\caption{The periodic cubic and quadratic splines used
for interpolation. The spline knots are indicated by \emph{blue full circles} and
the interpolation points, by \emph{dashed vertical lines} }
\label{fig:fitSpl}
\end{figure}
\subsection{\texttt{PP} representation}
The computation of $f(x)$ using directly the spline expansion Eq.~(\ref{eq:splExp}) can
be costly, because of the evaluation of the splines $\Lambda^p_j(x)$,
especially when interpolating on large number of points. Expanding $f(x)$,
using truncated Taylor series in each interval $[t_\mu,t_{\mu+1}]$, we obtain
the following \emph{Piecewise Polynomial Function} representation (or
\texttt{ppform}) of $f(x)$:
\begin{equation}
f(x) = \sum^p_{k=0}\, \Pi_{k\mu}(x-t_\mu)^k, \quad t_\mu\leq x<t_{\mu+1},
\label{eq:ppRep}
\end{equation}
where
\begin{equation}
\Pi_{k\mu} = \frac{1}{k!} \frac{d^k}{dx^k} f(t_\mu)
= \frac{1}{k!} \sum_j\,c_j\frac{d^k}{dx^k} \Lambda_j(t_\mu)
= \sum_j\,c_j\alpha_{j\mu k}.
\label{eq:ppCoef}
\end{equation}
Note that
\begin{equation}
\alpha_{j\mu k} = \frac{1}{k!} \frac{d^k}{dx^k} \Lambda_j(t_\mu)
\label{eq:derSpl}
\end{equation}
depend only on the spline specifications. They are \emph{pre-calculated}
in the spline setup routine \texttt{SET\_SPLINE} and stored in the
3d arrays \texttt{sp\%val0}. The \texttt{PP} coefficients
$\Pi_{k\mu}$ can be computed once the spline coefficients $c_j$ are available,
using (\ref{eq:ppCoef}).
Then the interpolated function values together with the $p$ derivatives
can be calculated \emph{efficiently} using the power series.
\begin{eqnarray*}
f(x) &=& \sum^p_{k=0}\, \Pi_{k\mu}(x-t_\mu)^k \\
f'(x) &=& \sum^p_{k=1}\, k\,\Pi_{k\mu}(x-t_\mu)^{k-1} \\
f''(x) &=& \sum^p_{k=2}\, k(k-1)\,\Pi_{k\mu}(x-t_\mu)^{k-2} \\
\vdots
\end{eqnarray*}
These 2 steps are performed in \texttt{GRIDVAL}. Note that the first step
(computation of $\Pi_{k\mu}$ from $c_j$) can be skipped for subsequent calls
to \texttt{GRIDVAL} with the same function $f$, for example to compute $f$ or
its derivatives at any others points $x$.
\subsection{Example 2: Cubic spline interpolation}
Given a function $f$ with its grid values \texttt{f(1:nx)}, the
following example determines the spline coefficients \texttt{c(1:dim)}. Using these
coefficients, the interpolated $f$ and derivative $f'$ are then computed on the mesh points
\texttt{xp(1:npts)}. Note that the second call to \texttt{GRIDVAL} does not
include the spline coefficients \texttt{c} to signal that the previously calculated
\texttt{PP} coefficients will be reused.
\par
\addvspace{\medskipamount}
\nopagebreak\hrule
\begin{verbatim}
USE BSPLINES
INTEGER, PARAMETER :: nx=10, npts=100
DOUBLE PRECISION :: x(0:nx), f(0:nx), xp(npts), fp0(npts), fp1(npts)
DOUBLE PRECISION, ALLOCATABLE :: c(:)
INTEGER :: dim
TYPE(spline1d) :: sp
!
! Define x and f
!
CALL set_splcoef(3, x, sp, period=.FALSE.)
CALL get_dim(sp, dim)
ALLOCATE(c(dim))
CALL get_splcoef(sp, f, c) ! compute spline coefs c(1:dim)
!
! Compute interpolated f and f' on xp(npts)
!
CALL gridval(sp, xp, fp0, 0, c)
CALL gridval(sp, xp, fp1, 1)
!
DEALLOCATE(c)
CALL destroy_sp(sp)
\end{verbatim}
\nopagebreak\hrule
\addvspace{\medskipamount}
The description of each of the routines called in the example above is briefly
given below:
\begin{description}
\item[\texttt{SET\_SPLCOEF}] Determines the spline knots, sets up the splines
with \texttt{SET\_SPLINE}, assembles the collocation matrix
and performs its \emph{factorization}.
\item[\texttt{GET\_SPLCOEF}] Computes the spline \emph{coefficients} $c$ from the input
grid values of function $f$ (\emph{back-solve phase}), using the factorized matrix.
\item[\texttt{GRIDVAL}] Compute the \texttt{PP} coefficients using
(\ref{eq:ppCoef}) \emph{if $c$ is provided}, locates the interval containing
the given point $x$ and computes interpolated function values or derivatives
using the \texttt{PP} representation (\ref{eq:ppCoef}).
\end{description}
\subsection{2d interpolation}
Consider the spline interpolation on the plane $(x,y)$, using a tensor
product of splines defined as follow
\begin{equation}
\Lambda^{p, q}_{ij}(x,y)=\Lambda^p_i(x)\Lambda^q_j(y), \qquad
\begin{aligned}[c]
i&=1,\ldots,d_1=N_1+p,\\
j&=1,\ldots,d_2=N_2+q,
\end{aligned}
\end{equation}
where $(p,q)$ are the spline degrees , $(N_1,N_2)$, the number of knot
intervals in each direction:
\[ t_0 \leq x < t_{N_1}, \quad s_0 \leq y < s_{N_2}.\]
\subsubsection{Spline coefficients (\texttt{GET\_SPLCOEF})}
The 2d version of (\ref{eq:intEq}) can be written as:
\begin{equation}
\begin{split}
\sum_{ij}\,c_{ij}\Lambda^p_{i}(x_\mu)\Lambda^q_{j}(y_\nu) &= f(x_\mu,y_\nu)\\
&=f_{\mu\nu}\\
\end{split}
\qquad
\begin{aligned}[c]
\mu&=0,\ldots,N_{g1}\\
\nu&=0,\ldots,N_{g2}
\end{aligned}
\end{equation}
where $(x_\mu,y_\nu)$ are the \emph{interpolation sites} on the $(x,y)$ plane.
These equations can be rearranged into
\begin{eqnarray}
\sum_i\,\bar{c}_{i\nu}\Lambda_i(x_\mu) &=& f_{\mu\nu},\\
\sum_j\,c_{ij}\Lambda_j(y_\nu) &=& \bar{c}_{i\nu}.
\end{eqnarray}
Such a 2 step procedure is implemented, using the 1d version
of \texttt{GET\_SPLCOEF} in the 2d version of
\texttt{GET\_SPLCOEF} by the following code fragment:
\par
\addvspace{\medskipamount}
\nopagebreak\hrule
\begin{verbatim}
TYPE(spline2d) :: sp
DOUBLE PRECISION :: ctr(SIZE(c,2), SIZE(c,1))
CALL get_splcoefn(sp%sp1, f, c)
CALL get_splcoefn(sp%sp2, TRANSPOSE(c), ctr)
c = TRANSPOSE(ctr)
\end{verbatim}
\nopagebreak\hrule
\addvspace{\medskipamount}
\subsubsection{\texttt{PP} representation (\texttt{GRIDVAL})}
Let us start with the \emph{spline representation}, for $f(x,y)$ with
$t_\mu\leq x < t_{\mu+1}$ and $s_\nu \leq y < s_{\nu+1}$:
\begin{equation}
f(x,y) = \sum_{j=1}^{d_2}\left( \sum_{i=1}^{d_1}\,c_{ij}\Lambda^p_i(x)\right)\Lambda^q_j(y).
\end{equation}
Applying successively (\ref{eq:ppRep}) to the $x$ and $y$ functional dependency
yields the following \texttt{PP} representation for $f(x,y)$:
\begin{equation}
f(x,y) = \sum^q_{k'=0}\left(\sum^p_{k=0}\,\Pi_{k\mu k'\nu}(x-t_\mu)^k\right)(y-s_\nu)^{k'}
\label{eq:pp2Rep}
\end{equation}
The 2d \texttt{PP} coefficient is the \emph{tensor product} of 2 1d
\texttt{PP} coefficients:
\begin{equation}
\Pi_{k\mu k'\nu} = \sum_{l'=1}^{d_2}\left(\sum_{l=1}^{d_1}\,c_{ll'}\,\alpha_{l\mu k}\right)\alpha'_{l'\nu k'}.
\label{eq:pp2Coef}
\end{equation}
where $\alpha$ and $\alpha'$ are respectively the derivatives (\ref{eq:derSpl}) of all splines
in $x$ and $y$ direction. All the derivatives of $f$ can be deduced straightforwardly
from the \texttt{PP} representation (\ref{eq:pp2Rep}).
In the present version, the 2d \texttt{GRIDVAL} can be called, either for (1) points
on a 2d structured mesh: \texttt{XP(1:NPX), YP(1:NPY)} and returns the array
\texttt{FP(1:NPX,1:NPY)}, or (2) with a 1d sequence of points
\texttt{XP(1:NP), YP(1:NP)} and returns the 1d array \texttt{FP(1:NP)}.
\section{Finite Elements using Splines}
\subsection{The weak form}
\subsection{The matrix assembly}
\subsection{The boundary conditions}
\subsubsection{Dirichlet condition}
Dirichlet BC can be simply applied by imposing the conditions on the
spline coefficients and the boundary point.
For the BC $u(x=0)=u_{1} = c$ for example,
the discrete equations can be expressed as:
\begin{equation}
\left(\begin{matrix}
1 & 0 & \cdots \\
A_{21} & A_{22} & \cdots \\
\vdots & \vdots & \ddots \\
\end{matrix}\right)
\left(\begin{matrix}
u_{1} \\
u_{2} \\
\vdots\\
u_{N} \\
\end{matrix}\right) =
\left(\begin{matrix}
c \\
f_{2} \\
\vdots\\
f_{N} \\
\end{matrix}\right)
\end{equation}
A more appropriate transformed system which
preserves any symmetry of the original system is:
\begin{equation}
\left(\begin{matrix}
1 & 0 & \cdots \\
0 & A_{22} & \cdots \\
\vdots & \vdots & \ddots \\
\end{matrix}\right)
\left(\begin{matrix}
u_{1} \\
u_{2} \\
\vdots\\
u_{N} \\
\end{matrix}\right) =
\left(\begin{matrix}
c \\
f_{2} -cA_{21}\\
\vdots\\
f_{N} -cA_{N1}\\
\end{matrix}\right)
\end{equation}
\subsubsection{Unicity on the axis}
Denoting the $N$ solutions at the axis by $(u_1, \ldots, u_N)$ , and
their transforms by $(\hat u_1, \ldots, \hat u_N)$ defined by
\begin{equation} \begin{array}{ccc}
u_1-u_N = \hat u_1 & & u_1 = \hat u_1 + \hat u_N \\
u_2-u_N = \hat u_2 & & u_2 = \hat u_2 + \hat u_N \\
\vdots & \Longrightarrow & \vdots \\
u_{N-1}-u_N = \hat u_{N-1} & & u_{N-1} = \hat u_{1-1} + \hat u_N \\
u_N = \hat u_N & & u_N = \hat u_N,
\end{array} \label{eq:unicity1} \end{equation}
the unicity condition can be specified by simply imposing
\begin{equation}
\hat u_1=\hat u_2=\ldots=\hat u_{N-1}=0. \label{eq:unicity2}
\end{equation}
From (\ref{eq:unicity1}), the \emph{transformation matrix} \(\mathbf U\) is defined
as
\begin{equation}
\mathbf{u} = \mathbf{ U \cdot\hat u}, \qquad \mathbf{U} =
\left(\begin{matrix}
1 & 0 & \dots & 0 & 1 \\
0 & 1 & \dots & 0 & 1 \\
& & \ddots& & \vdots \\
0 & 0 & \dots & 1 & 1 \\
0 & 0 & \dots & 0 & 1
\end{matrix}\right), \quad \mathbf{U^{T}} =
\left(\begin{matrix}
1 & 0 & \dots & 0 & 0 \\
0 & 1 & \dots & 0 & 0 \\
& & \ddots& & \vdots \\
0 & 0 & \dots & 1 & 0 \\
1 & 1 & \dots & 1 & 1
\end{matrix}\right).
\end{equation}
\paragraph{Matrix product \( \mathbf{A\cdot U}\)}
\begin{equation}
\mathbf{ A\cdot U} =
\left(\begin{array}{lllll}
A_{1,1} & A_{1,2} & \dots & A_{1,N-1} & \sum_{j} A_{1,j} \\
A_{2,1} & A_{2,2} & \dots & A_{2,N-1} & \sum_{j} A_{2,j} \\
& & \ddots& & \vdots \\
A_{N-1,1} & A_{N-1,2} & \dots & A_{N-1,N-1} & \sum_{j}A_{N-1,j} \\
A_{N,1} & A_{N,2} & \dots & A_{N,N-1} & \sum_{j}A_{N,j}
\end{array}\right).
\end{equation}
Thus \emph{right multiply by \(\mathbf{U}\)} is equivalent to put the
\emph{the sum of each row on the last column}.
\paragraph{Matrix product \( \mathbf{ U^T \cdot A}\)}
\begin{equation}
\mathbf{ U^T \cdot A} =
\left(\begin{array}{lllll}
A_{1,1} & A_{1,2} & \dots & A_{1,N-1} & A_{1,N} \\
A_{2,1} & A_{2,2} & \dots & A_{2,N-1} & A_{2,N} \\
& & \ddots& & \vdots \\
A_{N-1,1} & A_{N-1,2} & \dots & A_{N-1,N-1} & A_{N-1,N} \\
\sum_{i}A_{i,1} & \sum_{i}A_{i,2} & \dots & \sum_{i}A_{i,N-1} &
\sum_{i}A_{i,N}
\end{array}\right).
\end{equation}
Thus \emph{left multiply by \(\mathbf{\hat U}\)} is equivalent to put the
\emph{the sum of each column on the last row}.
\paragraph{Product \( \mathbf{\hat U \cdot b}\)}
\begin{equation}
\mathbf{\hat b} = \mathbf{U^T\cdot b} =
\left(\begin{array}{l}
b_1 \\
b_2 \\
\vdots \\
b_{N-1} \\
\sum_{i} b_{i}
\end{array}\right),
\end{equation}
\paragraph{Transformation of the original matrix equation}
The full original linear system, obtained from the discretization of the
2D \(r,\theta\) polar coordinates can be written as:
\begin{equation}
\left(\begin{array}{ll}
\mathbf{A} & \mathbf{B} \\
\mathbf{C} & \mathbf{D}
\end{array}\right)
\left(\begin{array}{l}
\mathbf{u} \\
\mathbf{v}
\end{array}\right) =
\left(\begin{array}{l}
\mathbf{b} \\
\mathbf{c}
\end{array}\right), \label{eq:orig_matrix_eq}
\end{equation}
where the solution array is split into the solutions \(\mathbf{u}\) at \(r=0\) and
the solutions \(\mathbf{v}\) on the remaining domain. The transformed system can
thus be written as
\begin{equation*}
\left(\begin{array}{ll}
\mathbf{U^T} & 0 \\
0 & \mathbf{I}
\end{array}\right)
\left(\begin{array}{ll}
\mathbf{A} & \mathbf{B} \\
\mathbf{C} & \mathbf{D}
\end{array}\right)
\left(\begin{array}{ll}
\mathbf{U} & 0 \\
0 & \mathbf{I}
\end{array}\right)
\left(\begin{array}{l}
\mathbf{\hat u} \\
\mathbf{v}
\end{array}\right) =
\left(\begin{array}{ll}
\mathbf{U^T} &0 \\
0 & \mathbf{I}
\end{array}\right)
\left(\begin{array}{l}
\mathbf{b} \\
\mathbf{c}
\end{array}\right),
\end{equation*}
\begin{equation}
\Longrightarrow
\left(\begin{array}{cc}
\mathbf{U^TAU} & \mathbf{U^TB} \\
\mathbf{CU} & \mathbf{D}
\end{array}\right)
\left(\begin{array}{l}
\mathbf{\hat u} \\
\mathbf{v}
\end{array}\right) =
\left(\begin{array}{c}
\mathbf{U^Tb} \\
\mathbf{c}
\end{array}\right),
\end{equation}
Notice that the transformation preserves any symmetry existing in the original system
(\ref{eq:orig_matrix_eq}). The transformed matrix is finally given in the following where
only the modified elements are shown and the sum is only over the first \(N\)
points in \(\theta\) direction. The \(\times\) symbol denotes unmodified elements.
\begin{equation}
\left(\begin{array}{lllllll}
\times & \times & \times & \times & \sum_{j} A_{1,j} & \times & \times \\
\times & \times & \times & \times & \sum_{j} A_{2,j} & \times & \times \\
\times & \times & \times & \times & \vdots & \times & \times \\
\times & \times & \times & \times & \sum_{j} A_{N-1,j} & \times & \times \\
\sum_{i}A_{i,1} & \sum_{i}A_{i,2} & \dots & \sum_{i}A_{i,N-1} &
\sum_{i,j}A_{i,j} & \sum_{i}A_{i,N+1} & \dots \\
\times & \times & \times & \times & \sum_{j} A_{N+1,j} & \times & \times \\
\times & \times & \times & \times & \vdots & \times & \times
\end{array}\right)
\end{equation}
Only the \(N^{th}\) column and the \(N^{th}\) row are affected by the transformation.
Applying now the unicity condition (\ref{eq:unicity2}) the final transformed system
reads:
\begin{equation}
\left(\begin{array}{lllllll}
1 & 0 & \dots & 0 & 0 & 0 & 0 \\
0 & 1 & \dots & 0 & 0 & 0 & 0 \\
0 & 0 & \ddots & 0 & \vdots & 0 & 0 \\
0 & 0 & \dots & 1 & 0 & 0 & 0 \\
0 & 0 & \dots & 0 & \sum_{i,j}A_{i,j} & \sum_{i}A_{i,N+1} & \dots \\
0 & 0 & \dots & 0 & \sum_{j} A_{N+1,j} & \times & \times \\
0 & 0 & \dots & 0 & \vdots & \times & \times
\end{array}\right)
\left(\begin{array}{l}
\hat u_1 \\
\hat u_2 \\
\vdots\\
\hat u_{N-1}\\
\hat u_{N} \\
u_{N+1} \\
\vdots
\end{array}\right) =
\left(\begin{array}{l}
0 \\
0 \\
\vdots\\
0 \\
\sum_{i} b_{i} \\
b_{N+1} \\
\vdots
\end{array}\right).
\end{equation}
\begin{thebibliography}{99}
\bibitem{deBoor} C. de Boor, \emph{A Practical Guide to Splines}, Applied
Mathematical Sciences, vol.~27 (Springer, NY, 2001).
\bibitem{Golub} G.H. Golub, C.F. Van Loan, \emph{Matrix Computation, 3rd
Edition}, p.5 (The John Hopkins University Press, 1996).
\end{thebibliography}
\end{document}
c4science · Help