\appendix

\section{Summary of abbreviations and notations}\label{sec:summ}
\begin{description}
\item[CZ:] convective zone
\item[EOS:] equation of state
\item[HR:] Hertzsprung-Russel
\item[IRK:] implicit Runge-Kutta scheme
\item[LMR:] limit between a mixed zone and a radiative zone
\item[MZ:] mixed zone, \ie, a convective zone contingently extended by 
overshooting
\item[PMS:] pre main-sequence
\item[ZAMS:] zero age main-sequence
\item[POST:] post main-sequence
\end{description}
\medskip
\begin{description}
\item[$n,\ n_{\rm a};\ n_{\rm X};\ n_\star$:] number of shells for internal 
structure, 
in the atmosphere; number of chemical species; index of the shell, in the 
atmosphere where the radius of the star is defined

\item[$X,\ Y,\ Z$:] mass fraction in mass unit of, hydrogen, helium and heavy 
elements

\item[$X_i,\ \bar X_i;\ {\cal X},\ \bar{\cal X};\ {\cal R}_{Z0}$:] 
abundance for the chemical species \# $i$; 
mean abundance; vector of chemical species; vector of mean chemical species; 
initial ratio, in mass, of heavy elements to hydrogen contents 

\item[$P,\ T,\ R,\ L,\ M,\ t,\ \tau$:] Pressure, temperature, radius, 
luminosity, mass, time, Rosseland optical depth

\item[$P_{\rm b},\ T_{\rm b},\ R_{\rm b},\ L_{\rm b},\ M_{\rm b},\ \tau_{\rm 
b},\ \tau_i$:] Pressure, temperature, radius, luminosity, mass, Rosseland 
optical depth 
at the bottom of the atmosphere and at the shell \#$i$

\item[$R_\star,\ L_\star,\ \tau_\star$:] radius, luminosity of the star, 
Rosseland optical depth at the level taken as the radius of the star

\item[$T_{\rm eff},\ P_{\rm ext} ,\ \tau_{\rm ext},\ {\cal T}$:] effective 
temperature, pressure, optical depth at the external limit of the atmosphere, 
symbol for $T(\tau)$ law

\item[$\ M_\odot,\ R_\odot,\ L_\odot,\ t_\odot,{\cal R}_{Z\odot}$:] at present 
time, mass, radius, luminosity, age and ratio, in mass, of hydrogen to heavy 
elements for the Sun

\item[$M_\star,\ M_{\rm ini};\ \dot M,\ \dot{\cal M}$:] total and initial mass
of the star; external 
mass loss rate, change of mass due to mass defect

\item[$\xi,\ \eta,\ \zeta,\ \zeta_\star,\ \lambda,\ \mu,\ \nu$:] variables used 
in the calculations for pressure, temperature, radius, luminosity, mass or mean 
molecular weight

\item[$\xi_{\rm b},\ \eta_{\rm b},\ \zeta_{\rm b},\ \lambda_{\rm b},\ \mu_{\rm 
b}$:] variables used in the calculations for pressure, temperature, radius, 
luminosity, mass at the bottom of the atmosphere

\item[$\Delta_\xi,\ \Delta_\mu;\ \Delta_i,\ \Delta_g$:] 
distribution factors for pressure and mass; maximum changes of the chemical 
species labeled by $i$, of gravothermal energy

\item[$T_6,\ A;\ \sigma;\ {\rm d}M_{\rm d},\ {\rm d}\mu_{\rm d}$:] variables of 
OPAL opacity tabulations; 
Stefan-Boltzman constant; directions, in mass, taking into account the mass 
defect

\item[$\nabla_M,\ \nabla,\ \nabla_{\rm ad},\ \nabla_{\rm rad}$:] gradient 
operator with respect to the mass, gradient, adiabatic gradient, 
radiative gradient

\item[$\rho,\ S,\ U,\ \Omega,\ {\cal N},\ \Gamma_1$:] density, specific entropy, 
specific internal energy, angular velocity, Brunt-V\"ais\"ala frequency, adiabatic exponent

\item[$\kappa,\ \epsilon,\ \epsilon_G,\ \epsilon_\nu,\ {\cal C},\ c$:] 
Rosseland mean opacity, nuclear 
energy rate, gravothermal energy, energy of neutrinos, velocity of light, 
contraction factor


\item[$\varphi,\ \varphi_\star;\ \Delta\varphi^\pm$:] variable index for the 
reconstitution of the atmosphere; derivatives in the upper and lower parts of 
the atmosphere

\item[$Q,\ C(t),\ q,\ \psi,\ \theta,\ \omega$:] repartition function,
repartition constant, index real, 
derivatives of the repartition function, weight function 

\item[$Q_i,\ \psi_i,\ \theta_i,\ \omega_i$:] repartition function, derivatives
of the repartition function, weight function for the shell \#$i$

\item[$F_i,\ F_i^\pm$:] flow, incoming/outcoming flow in a MZ of the chemical 
species \#$i$

\item[$F_\Omega^\pm,\ {\cal M}_\Omega,\ \dot{\cal M}_\Omega$,\ $R_{\rm A}$:] 
incoming/outcoming flow of angular momentum in a MZ, angular momentum, lost rate 
of 
angular momentum, Alf\`en radius

\item[$\Psi_i,\ J,\ \lambda_{\rm min},\ \lambda_{\rm max}$:] creation rate of 
the chemical species \#$i$, jacobian matrix, minimum and maximum eigenvalues of
the jacobian matrix

\item[${\cal D}_i,\ d_{i,j},
\ d^\star_{i,j},\ d_{\rm T},\ d_{\rm M},\ v_i $:] vector of diffusion 
coefficients, total diffusion 
coefficient, microscopic diffusion coefficient of the chemical species \#$i$ with 
respect to \#$j$, turbulent and mixing diffusion coefficients,
velocity of diffusion for the chemical species \#$i$

\item[$l,\ G,\ g,\ g_\star$:] mixing-length parameter or, also, index, 
gravitational 
constant, local gravity, total gravity of the star

\item[$h,\ O_j^h,\ h_q,\ h_c$:] Order of B-splines, $j$-ith B-spline of order 
$h$, order of B-splines for the integration of the quasi-static equilibrium, for 
the interpolation of the chemical species 

\item[${\cal A},\ a_{i,j}\ c_i,\ b_i,\ d_i,\ p,\ s$:] matrix, coefficients, 
order and number of steps of the IRK scheme

\item[$y,\ y_0,\ \bar y,\ \bar y_j,\ y_{j,k+1/2},\ z_i,\ \bar z_i$:] abundances 
at several levels of the IRK scheme

\item[$\nu_{k+1/2},\ \Delta\nu_{k+1/2},\ T_{i,k+1/2},\ \rho_{i,k+1/2}$:] 
intermediate times, interval between intermediate times, temperature, density 
at intermediate times in the IRK scheme

\item[$\Delta t,\ \Delta t_{\rm opt},\ \Delta t_{\rm new}$:] time step, optimal 
and new time steps

\item[$\varepsilon_{\rm r},\ \varepsilon_{\rm e}$:] required and estimated 
accuracy 
for the integration of chemical composition 

\end{description}


\section{}
\subsection{Sets of variables without central singularities}\label{sec:ap1}
With the lagrangian independent variable $\mu\equiv M^\alpha$, and for the 
variables
$\zeta\equiv R^\beta$ and $\lambda\equiv L^\delta$ in place of $R$ and $L$,
one seeks exponents $\alpha$, $\beta$ and $\delta$ in such a way
that the gradients of pressure, radius and luminosity
are $O(1)$ in a neighborhood of the center, $M=0$.
The gradients are written:
\begin{eqnarray*}
 \dsdx P\mu &=& -\frac{G}{4\pi\alpha}\frac{\mu^{2/\alpha-1}}
 {\zeta^{4/\beta}}, \\
 \dsdx\zeta\mu &=& \frac \beta{4\pi\rho\alpha}\frac{\mu^{1/\alpha-1}}
 {\zeta^{3/\beta-1}}, \\
 \dsdx\lambda\mu &=& \frac\delta\alpha \left[\epsilon -T\dsdx S t \right]
 \frac{ \mu^{1/\alpha-1}} {\lambda^{1/\delta-1}}.
\end{eqnarray*}
\noi In the vicinity of the center the density, $\rho$, is $\rho=O(1)$,
therefore $R=O(M^{1/3})$ and $\zeta=O(\mu^{\beta/3\alpha})$; similarly,
$L=O(M)$ and $\lambda=O(\mu^{\delta/\alpha})$ and, at the limit $\mu=0$:
\begin{eqnarray*}
 \dsdx P\mu      &=& O(\mu^{2/3\alpha-1}) \equiv O(\mu^p), \\
 \dsdx\zeta\mu   &=& O(\mu^{\beta/3\alpha-1}) \equiv O(\mu^r), \\
 \dsdx\lambda\mu &=& O(\mu^{\delta/\alpha-1}) \equiv O(\mu^l).
\end{eqnarray*}
The singularities at $M=0$ are avoided if the exponents $p$, $r$ and $l$ 
are non negative integers; in \tab{tab:expl} the exponents $\alpha$, $\beta$ 
and 
$\delta$, deduced from the previous set of equations, are given for the values 
0, 1 and 2, for $p$, $r$ and $l$.

\begin{table}
\caption[]{With a lagrangian independent variable, the exponents
$\alpha$, $\beta$ and $\delta$ are given for the values 0, 1 and 2, for $p$,
$r$ and $l$. As seen
the optimal accuracy ($p=r=l=0$), is reached only with non integer
exponents for the mass and the luminosity, \ie, with $\alpha=2/3$, $\beta=2$ and 
$\lambda=2/3$.}
\label{tab:expl}
\begin{tabular}{lllllllllll} \hline
&& $\alpha$ & \multicolumn{4}{c}{$\beta$} & \multicolumn{4}{c}{$\delta$} \\
$p$ &&&& $r=0$ & $r=1$ & $r=2$ && $l=0$ & $l=1$ & $l=2$ \\ \hline

0 && 2/3  && 2   & 4   & 6 && 2/3 & 4/3 & 2 \\
1 && 1/3  && 1   & 2   & 3 && 1/3 & 2/3 & 1 \\
2 && 2/9  && 2/3 & 4/3 & 2 && 2/9 & 4/9 & 2/3 \\ \hline

\end{tabular}
\end{table}

\noi From \tab{tab:expl}, with the exponents 
$\alpha=2/3$, $\beta=2$ and $\delta=2/3$, \ie, with the variables $M^{2/3}$, 
$R^2$ and
$L^{2/3}$, one ensures $p=0$, $r=0$ and $l=0$, therefore $\nabla P=O(1)$, 
$\nabla\zeta=O(1)$ and $\nabla\lambda=O(1)$ as expected. 

Obviously other sets of exponents can be employed, the best of them is 
chosen according to the following numerical criterion:
working with a scheme of, let say, third order, between two grid
points the dependent variables are quadratic polynomials
in $\mu$; any of them, $\xi$ for example, is written:
$\xi=\xi_0+\xi_1\mu+\xi_2\mu^2+O(\mu^3)$. In a neighbourhood of $\mu=0$, 
the gradient are $O(\mu^p)$; with $p=0$ one has: 
\begin{displaymath}
 \dsdx\xi\mu=O(\mu^0)=O(1)\sim\xi_1+2\xi_2\mu+O(\mu^2),
\end{displaymath}
\noi while, if $p=1$:
\begin{displaymath}
 \dsdx\xi\mu=O(\mu^1)\sim 2\xi_2\mu+O(\mu^2).
\end{displaymath}
\noi The last approximation only uses the term in $\mu^2$, therefore
information is lost. Hence the set of exponents
$\alpha=2/3$, $\beta=2$ and $\delta=2/3$, \ie, the variables $M^{2/3}$, $R^2$ 
and $L^{2/3}$, are optimal.

\begin{table}
\caption[]{The same as Table~\ref{tab:expl} but for eulerian
variables. Note that the optimal accuracy ($p=r=l=0$) is not reached
with $\alpha=\beta=\delta=1$, \ie, with the variables used in \eq{eq:eq5}.}
\label{tab:expe}
\begin{tabular}{lllllllllll} \hline
&& $\beta$ & \multicolumn{4}{c}{$\alpha$} & \multicolumn{4}{c}{$\delta$} \\
$p$&&&& $r=0$ & $r=1$ & $r=2$ && $l=0$ & $l=1$ & $l=2$ \\ \hline

0 && 2   && 2/3 & 4/3 & 2   && 2/3 & 4/3 & 2 \\
1 && 1   && 1/3 & 2/3 & 1   && 1/3 & 2/3 & 1 \\
2 && 2/3 && 2/9 & 4/9 & 2/3 && 2/9 & 4/9 & 2/3 \\ \hline

\end{tabular}
\end{table}

\noi From table~\ref{tab:expe}, where the exponents for
eulerian variables are given, one sees that
the highest accuracy cannot be obtained with the set of variables of
\eq{eq:eq5}.

\subsection{Setting a grid point at a given location.}\label{sec:set}
At each LMR the derivative of the temperature gradient is a discontinuous 
function. The numerical technique used for the integration of the equations of
internal structure is based on an approximation of the unknown functions by 
piecewise polynomials (see \App{sec:ide}).
Between two mesh points, each polynomial piece is a differentiable function, 
therefore their derivatives are allowed to be discontinuous only at grid points
where the polynomial pieces are connected; as a consequence, each LMR must be 
localized precisely on a knot; this allows, with an {\it ad hoc} choice of the 
basis, to have different values on left and right for the derivatives of the
unknown functions.
For sake of clarity in the discussion and the notation, one assumes that there
is only one LMR where a discontinuity arises (labeled hereafter with the
subscript ``d''); the extension to several discontinuities is straightforward. 
The
location 
of each discontinuity on a knot is obtained by a slight change of the 
algorithm used for the automatic mesh refinement; it is, basically,
an inverse linear interpolation: The grid points are localized in such a way
that, from a mesh to the next, the change of the
distribution function is constant; inside each mesh, indexed by $i$,
the relative weights of the variables are defined by the distribution 
factors $\Delta_\xi$ and $\Delta_\mu$ known within a factor 
$\omega_i(q)$, 
$i=1,2,\ldots,n-1$, which was, up to now,
implicitly assumed to be the unity for each mesh, but it could be
adjusted (iteratively)
in order to set each LMR on the closest grid point; in other words, locally, for 
both sides of each LMR, instead of using the constant value $C(t)$
the repartition constant is adjusted to the amounts just needed
for fitting a grid point right on the LMR.
Doing that, instead of $Q(q)$, for $q\in [q_i,q_{i+1}[$,
one uses a "weighted distribution function", defined by pieces, as:

\begin{equation}\label{eq:qi}
Q_i(q)=Q(q)/\omega_i(q),\ \omega_i(q)\neq 0. 
\end{equation}
Analogously one also defines:

\begin{displaymath}
\psi_i(q)=\dsdx{Q(q)}q/\omega_i(q),\ \theta_i(q)=\dsdx{Q(q(\mu))}{\mu} 
/\omega_i(q).
\end{displaymath}
\noi The pieces of the function $\psi_i$ being connected by continuity:
 $\psi_i(q_i)=\psi_{i-1}(q_i)$, $i=2,\ldots,n-1$.
Therefore, instead of \eq{eq:muq}, one uses:

\begin{equation}\label{eq:mui}
\left( \dsdx{\mu}q \right)_t=\frac{\psi_i}{\theta_i} ,\qquad
\left( \dsdx{\psi_i}q \right)_t=0.
\end{equation}
Note that $\omega_i(q)$ and $\theta_i(q)$ are step functions. 
Due the facts (see Appendix~\ref{sec:ide}), that at grid points first, the 
differential equations are not written and, second, discontinuous derivatives
are allowed for, the discontinuities of $\partial\psi_i(q)/\partial q$
are implicitly managed.
Let, at time $t$, a solution of \eq{eq:syq} be known and $Q_i^{(0)}(q)$ be
the weighted distribution function defined by \eq{eq:qi} with 
$\omega_i^{(0)}(q)\equiv 1$. Here each quantity with a `$(j)$' superscript is 
related to the $j$-th iterated model, $j=0,1,\ldots$ At the location 
$q^{(j)}_\LoD$, where the discontinuity occurs, the distribution function, 
calculated with \eq{eq:qrp}, amounts to  $Q^{(j)}_\LoD$; let $l$ be the index
of the nearest mesh point from the discontinuity ($l\not=1,n$) (in this
discussion 
one simplifies by assuming that the index $l$ does not change as it does, 
sometimes, in the calculations); since 
$Q^{(j)}(q)$ is a strictly monotonous function, there is only a
unique value $q^{(j)}_\LoD$ such that $Q^{(j)}(q^{(j)}_\LoD)=Q^{(j)}_\LoD$,
therefore, for $q\in [q_i,q_{i+1}[$, $i=1,2,\ldots,n-1$, the weights 
$\omega^{(j)}_i(q)$ can be defined by: $\omega^{(0)}_i=1$ and, for $j\geq 1$:
\begin{displaymath}
\begin{array}{l}
{\omega^{(j)}_i(q)=} \\ \\
{\left\{
\begin{array}{ll}
\left[Q^{(j-1)}(q_\LoD)-Q^{(j-1)}(q_{l-1})\right]\ /\ \psi^{(j-1)}_i,&
 \mbox{if $i=l-1$}, \\
\left[Q^{(j-1)}(q_{l+1})-Q^{(j-1)}(q_\LoD)\right]\ /\ \psi^{(j-1)}_i,&
 \mbox{if $i=l$}, \\
1,& \mbox{otherwise.}
\end{array}
\right.} 
\end{array}
\end{displaymath}
\noi Hence, at the iteration $j\geq 1$, the LMR will be
closer, than previously, to the grid point with the index $l$.
After convergence ($j=\infty$), $q^{(\infty)}_D=q_l$,
$\psi^{(\infty)}_i(q)={\sl Cte}=Q^{(\infty)}_{l-1}(q_l)
-Q^{(\infty)}_{l-1}(q_{l-1})$.
The jump of the distribution function, between $q_{l-1}$ and 
$q_l=q^{(\infty)}_D$ will satisfy: 
\begin{displaymath}
\begin{array}{l}
Q^{(\infty)}(q_l)-Q^{(\infty)}(q_{l-1})
=Q^{(\infty)}(q^{(\infty)}_D)-Q^{(\infty)}(q_{l-1})=\\
\left[Q^{(\infty)}_{l-1}(q^{(\infty)}_D)-Q^{(\infty)}_{l-1}(q_{l-1})\right]
\omega^{(\infty)}_{l-1}(q_{l-1}), \\
\end{array}
\end{displaymath}
just that is needed.
For $j=0$ the weights $\omega^{(j)}_i(q)$ are initialized to unity and
updated between each Newton-Raphson iteration of the quasi-static equilibrium 
problem with $Q^{(j)}(q)$ and $q^{(j)}_D$ derived from the locations of the
LMR at the $(j-1)$-th iteration, for $j>0$.
\Fig{fig:LMR} shows a typical displacement of a grid point
toward a LMR (see also \fig{fig:vais}). In this simplified case only
three iterations have been necessary to set a grid point on the LMR.
\begin{figure}
\defbox{\psfig{figure=LMR.ps,%
height=7cm,bbllx=0.5cm,bblly=0.5cm,bburx=27.5cm,bbury=20.5cm,clip=}}
\caption{Automatic adjustment of a grid point on the limit between a radiative
and a convective zone.
Two iterations have been needed to push a grid point towards the LMR; during 
the process the radiative gradient $\nabla_{\rm rad}$ slightly decreases and
the LMR moves 
right. The adiabatic gradient $\nabla_{\rm ad}$ , full line, remains unchanged. 
The initial 
locations of the grid points (empty squares), are moved to the right at the
first iteration (empty circles) and slightly right again at the second on
(full circles).}\label{fig:LMR}
\end{figure}
\noi This iterative process is only of first order. It will converge
only if the 
functions considered are sufficiently smooth and, moreover, if the LMR is 
sufficiently well defined, as in the case illustrated in \fig{fig:LMR}; it is 
the case
for the external LMR of the Sun and for the convective core of stars of 
intermediate mass. On the other hand, for semi-convective zones or convective 
cores with diffuse LMR, as for the convective core of the young Sun, the 
closest 
index $l$, from the LMR, is not well fixed and the algorithm may not converge; 
fortunately, in such cases, the jump of the slope of the gradient is almost 
zero, therefore a mislabelling of a small jump in the slope of
the gradient has almost no consequence.

\subsection{Precise external boundary conditions.}\label{sec:pbc}
As described \Sec{sec:ebc}, the precise reconstitution of an atmosphere consists
in a differential problem, with boundary conditions at three different levels:
$\tau_b\goa 10$, $\tau_{\rm min}\loa 10^{-4}$ and $0.3\loa\tau_\star\loa 0.6$; 
recall that, at $\tau=\tau_\star(t)$, defined by \eq{eq:star1}, there is a {\em 
open} 
inner limit. In the numerical solution of \eq{eq:sysa}, one assigns, among the 
$n_{\rm a}$ grid points allowed for the atmosphere, a {\em fixed} index 
$n_\star$ 
to the {\em moving} limit $\tau_\star$; that is done by the piecewise 
logarithmic linear mapping: $\tau\mapsto\omega$ of $[\tau_{\rm b},\tau_{\rm 
min}]$ in $[1,n_{\rm a}]$ defined by:
\begin{equation}\label{eq:map}\displaystyle
\begin{array}{l}
{\ln\tau\mapsto\varphi(\omega)=} \\
{\cases{
\ln\tau_{\rm b}+(\omega-1)\Delta\varphi^+,&if $\omega\in[1,n_\star]$;\cr
\ln\tau_{\rm min}+(\omega-n_{\rm a})\Delta\varphi^-,&otherwise.\cr}}
\end{array}
\end{equation}
\noi Here the slopes $\Delta\varphi^+$ and $\Delta\varphi^-$ are defined by:
\begin{displaymath}
\Delta\varphi^+\equiv\frac{\ln\tau_\star-\ln\tau_{\rm min}}{n_\star-n_{\rm a}},
\ \Delta\varphi^-\equiv\frac{\ln\tau_\star-\ln\tau_{\rm b}}{n_\star-1};
\end{displaymath}
\noi hence one has:
\begin{description}
 \item $\omega=1$, at the bottom of the atmosphere where $\tau=\tau_{\rm b}$,
 \item $\omega=n_\star$, at the inner limit where $\tau=\tau_\star$,
 \item $\omega=n_{\rm a}$, at the upper
  limit of the atmosphere where $\tau=\tau_{\rm min}$;
\end{description}
\noi the index $n_\star$ is a open parameter; with $n_\star\simeq n_{\rm a}/2$, 
$\ln\tau_\star$ will lie in the vicinity of the middle of the interval
$[\ln\tau_{\rm b},\ln\tau_{\rm min}]$. From
this mapping results a constant mesh size, it is an advantage of the algorithm.
For the numerical solution it is convenient to add one more unknown: the 
temperature $T$, which fulfills the {\em scalar} equation $T  =  {\cal 
T}(\tau,T_{\rm eff},g)$, and the boundary condition
${\cal T}(\tau_\star,T_{\rm eff},g)=T_{\rm eff}$. With  the following 
variables, employed for the reconstitution of the atmosphere:
$\xi=\ln P$, $\eta=\ln T$, $\zeta_{\rm a}=R/R_\odot $,
$\mu_{\rm a}=M/M_\odot$, 
$\zeta_\star=R_\star/R_\odot$ and $\varphi_\star=\ln\tau_\star$, \eq{eq:sysa} 
are written:

\begin{equation}\displaystyle
\label{eq:sysa2}
\left\{
\begin{array}{l}\vspace{0.05 in}
\ddx \xi\omega=\left[\frac{GM_\odot }{R_\odot^2}\frac{\mu_{\rm a}}{\zeta_{\rm 
a}^2}
-\frac{2R_\odot\Omega^2\zeta_{\rm 
a}}3\right]\frac{\Delta\varphi^\pm}{\kappa}\exp(\varphi-\xi)
\\ \vspace{0.05 in}

\eta=\ln\left({\cal T}(\exp(\varphi),T_{\rm eff},g_\star)\right) \\
\vspace{0.05 in}

\ddx {\zeta_{\rm a}}\omega  = -\frac 1{R_\odot 
}\frac{\Delta\varphi^\pm\exp(\varphi)}
{\kappa\rho} \\ \vspace{0.05 in}

\ddx {\zeta_\star}\omega  = 0 \\ \vspace{0.05 in}

\ddx {\mu_{\rm a}}\omega  = -\frac{4\pi R_\odot ^2}{M_\odot }
\frac{\Delta\varphi^\pm}\kappa\exp(\varphi)\zeta_{\rm a}^2 \\ \vspace{0.05 in}

\ddx {\varphi_\star}\omega =  0 \\ \vspace{0.05 in}
\ddx \varphi\omega  =  \Delta\varphi^\pm \\
\end{array}
\right.
\end{equation}
\noi where $\Delta\varphi^\pm=\Delta\varphi^+$ \respt{$\Delta\varphi^-$} if 
$\omega\in[1,n_\star]$ \respt{$\omega\in[n_\star,n_{\rm a}]$}. The boundary 
conditions are:

\begin{displaymath}\displaystyle
\left\{
\begin{array}{l}
\zeta_{\rm a}(1)  =  \zeta_{\rm b} \\ \vspace{0.05 in}

\mu_{\rm a}(n_\star)  =  {M_\star}/{M_\odot } \\ \vspace{0.05 in}

\zeta_\star(n_\star)  =  {R_\star}/{R_\odot } \\ \vspace{0.05 in}

\eta(n_\star)=\ln{\cal T}(\exp(\varphi_\star),T_{\rm eff},g_\star) \\ 
\vspace{0.05 in}

\eta(n_\star)=\ln T_{\rm eff} \\ \vspace{0.05 in}

\varphi(n_{\rm a})=\ln \tau_{\rm min} \\ \vspace{0.05 in}

\rho(n_{\rm a})  = \rho_{\rm ext} \\
\end{array}
\right.
\end{displaymath}
\noi here:
\[
T_{\rm eff}=\sqrt[4]{\frac{L_\odot\sqrt{\lambda^3_{\rm b}}}{4\pi R_\odot 
^2\sigma\zeta_\star^2}},
\quad g=\frac{GM_\odot\mu_{\rm a}}{R_\odot ^2\zeta_{\rm a}^2}.
\]
\noi The numerical solution of the non-linear differential pro\-blem
Eq.~(\ref{eq:sysa2}), with its boundary conditions at three different 
levels, is solved using the B-spline collocation method
(see Appendix~\ref{sec:ide}). The iterative process is initialized either, 
with an approximate model of the solar atmosphere for PMS or ZAMS initial
models or, along the evolution, with the atmosphere of the previous model .
The solution gives $\xi_{\rm b}(\zeta_{\rm b},\lambda_{\rm b})$,
$\eta_{\rm b}(\zeta_{\rm b},\lambda_{\rm b})$,
$\mu_{\rm b}(\zeta_{\rm b},\lambda_{\rm b})$ and their derivatives with respect
to $\eta_{\rm b}$ and $\lambda_{\rm b}$.
Ty\-pi\-cal\-ly, one uses from $n_{\rm a}=30$ to $n_{\rm a}=50$ grid points for
the interval $[\tau_{\rm b},\tau_{\rm min}]$.
With B-splines of order 3, there are from 59 to 99 points where 
Eqs.~(\ref{eq:sysa2}) are fulfilled; due to the 
superconvergence (see \sec{sec:ide}), the order of the scheme is 4.

 

\begin{table}
 \caption[]{Coefficients in \eq{eq:IRK1} for the IRK formulas Lobatto IIIC
 (see text) of order $p=1,2,4$  with $s=1,2,3$ steps. 
 As usual, the $c_i,\ i=1,\ldots ,s$ are in the first column, the $b_i$
 in the last row, the remainder entries are the coefficient $a_{ij}$
 of the IRK matrix ${\cal A}$.}
 \label{tab:lob}
 \hbox{
 \begin{tabular}{c|c}
  1 & 1 \\ \\
  \hline \\
  & 1 \\
 \end{tabular}
 \qquad\qquad

 \begin{tabular}{c|cc}
  0 & $\frac12$ & $-\frac12$ \\ \\
  1 & $\frac12$ & $\frac12$ \\ \\
  \hline \\
  & $\frac12$ & $\frac12$
 \end{tabular}

 \qquad\qquad
 \begin{tabular}{c|ccc}
  0 & $\frac16$ & $-\frac13$ & $\frac16$\\ \\
  $\frac12$ & $\frac16$ & $\frac5{12}$ & $-\frac1{12}$\\ \\
  1 & $\frac16$ & $\frac23$ & $\frac16$\\ \\
  \hline \\
   & $\frac16$ & $\frac23$ & $\frac16$
  \end{tabular}
}
\end{table}
\subsection{Stable integration with IRK formulas.}\label{sec:sirk}
With stiff problems, care must be taken for saving the stability;
the procedures are briefly recalled here in the case of only one chemical 
element; let $y_0$ \respt{$y$} be the abundance of that 
element at time $t$ \respt{$t+\Delta t$}.
With the coefficients reproduced in \tab{tab:lob}, the IRK scheme is written:
 \begin{equation}\label{eq:IRK1}
 \begin{array}{l}
  y_i=y_0+\Delta t\sum_{j=1}^sa_{ij}\Psi(T_j,\rho_j,y_j),\ i=1,\ldots, s \\
  y=y_0+\Delta t\sum_{j=1}^sb_j\Psi(T_j,\rho_j,y_j)
 \end{array}
\end{equation}
where $T_j$, $\rho_j$ and $y_j$ are, respectively,
the temperature, density and abundance at intermediate times
$t+c_j\Delta t$, $j=1,\ldots, s$; $s$ is the number of steps of the IRK 
formulae.
This standard formulation for IRK schemes does not give satisfactory
results for the evolution of the chemical composition: negative values for
the abundances are observed; that is
due to the amplification of the errors by large Lipschitz constants.
The round-off errors are reduced
if the first \eq{eq:IRK1} is rewritten
for the smaller quantities $z_i=y_i-y_0$ (Hairer \& Wanner, \loc, \S IV.8) as:
\begin{equation}\label{eq:IRK21}
 \left(
 \begin{array}{c} 
  z_1 \\ z_2 \\ \vdots \\ z_s
 \end{array}
 \right)\ =\ \Delta t\ {\cal A}\ \left(
 \begin{array}{c}
  \Psi(T_1,\rho_1,y_0+z_1) \\
  \Psi(T_2,\rho_2,y_0+z_2) \\
  \vdots \\ 
  \Psi(T_s,\rho_s,y_0+z_s)
 \end{array}
 \right).
\end{equation}
The matrix ${\cal A}\equiv(a_{i,j})_{i,j=1}^s$ of any IRK Lobatto IIIC 
formulas is invertible therefore, with
$(d_1,d_2,\ldots,d_s)\equiv (b_1,b_2,\ldots,b_s){\cal A}^{-1}$,
one has\footnote{This property only results from the fact that the last row of 
$A$ and the vector $b$ are identical, \eg all the Lobatto IIIC formulae are 
established in such a way.} $d_j=0$ if $j\not= s$ and $d_s=1$;
then the second \eq{eq:IRK1} is simply written:
\begin{equation}\label{eq:IRK22}
 y=y_0+z_s.
\end{equation}

\paragraph{Instantaneous mixing with IRK formulas.}
In MZ the abundances have constant values as described \sec{sec:mixc};
for the chemical species $y$, let $\bar y$ be this mean. At time $t+c_j\Delta t,\ 
j=1,\ldots, s$, the mean, $\bar y_j$, is approximated using the mid-point rule:
\[
 \bar y_j=\sum_{k\in{\rm MZ}}y_{j,k+1/2}\Delta\nu_{k+1/2},
\]
\noi here the sum involves all the meshes inside the MZ; $y_{0,k+1/2}$
is the abundance at the mid-point, $\nu_{k+1/2}\equiv(\nu_{k+1}+\nu_k)/2$
and: 
\begin{equation}\label{eq:dnu}
\Delta\nu_{k+1/2}\equiv\frac{ (\nu_{k+1}-\nu_k) \sqrt{\nu_{k+1/2}} }
{ \sum_{ j\in{\rm MZ} } (\nu_{j+1}-\nu_j) \sqrt{ \nu_{j+1/2} }  }.
\end{equation}
Therefore Eq.~(\ref{eq:IRK21}) is written:
\begin{equation}\label{eq:IRKm1}\displaystyle
\begin{array}{l}
{
 \left(
 \begin{array}{c} 
  \bar z_1 \\ \bar z_2 \\ \vdots \\ \bar z_s
 \end{array}
 \right)=\Delta t\ {\cal A} 
\ \times} \\ \\
{\times 
\left(
 \begin{array}{c}
  \sum_k\Psi(T_{1,k+\frac 12},\rho_{1,k+\frac 12},\bar y_0+\bar z_1)
\Delta\nu_{k+\frac 12} \\
   \sum_k\Psi(T_{2,k+\frac 12},\rho_{2,k+\frac 12},\bar y_0+\bar z_2)
\Delta\nu_{k+\frac 12} \\
  \vdots \\ 
   \sum_k\Psi(T_{s,k+\frac 12},\rho_{s,k+\frac 12},\bar y_0+\bar z_s)
\Delta\nu_{k+\frac 12}
 \end{array}
 \right) 
} \end{array}
\end{equation}
where $T_{j,k+\frac 12}$ \respt{$\rho_{j,k+\frac 12}$} is the temperature 
\respt{the density} at the intermediate point at time $t+c_j\Delta t$;
$\bar y$ then is written:
\begin{equation}\label{eq:IRKm2}
\bar y=\bar y_0+\bar z_s.
\end{equation}
The extension to several chemical elements of \eq{eq:IRK21} to 
\eq{eq:IRKm2} is straightforward. Note that, with the use of the 
mid-point rule, \eq{eq:dnu} for 
the numerical approximation the integral of \eq{eq:dx2}
over the MZ, is only of second order; a two-points Gaussian formula
(\ie third order accuracy) will be more appropriate. The error is not
important in the solar case since the thermonuclear 
reactions are not really efficient in MZ.

In CESAM, \eq{eq:IRK21} to \eq{eq:IRKm2} are solved iteratively using
Newton-Raphson scheme\footnote{due to the stiffness, the Lipschitz constants
are large and predictor-corrector methods
are useless.} with an analytical jacobian updated at each iteration;
the initial values of the $z_i$ are fixed to zero; the quadratic iterative
process is halted as soon as the corrections on the $z_j$ are less
than $10^{-8}$ in relative values; most of the time, the
convergence is reached in less than five iterations.

\paragraph{Conservation of nucleus and charge numbers.}\label{sec:cirk}
For a number of chemicals $n_{\rm X}>1$, let at time $t$:
\[
 \sum_{i=1}^{n_{\rm X}}\alpha_i X_i^t=cte,
\]
\noi be any linear relationship between the abundances
with {\em fixed} coefficient $\alpha_i$, \eg conservation of baryon and charge 
numbers are relationships of that form, then:

\[
 \dsdx{\ }t\sum_{i=1}^{n_{\rm X}}\alpha_i X_i=0\Longrightarrow
 \sum_{i=1}^{n_{\rm X}}\alpha_i\Psi_i=0.
\]
\noi Owing to the linearity of \eq{eq:IRK21} and \eq{eq:IRKm1}, for $n_{\rm X}$
elements, with obvious notations:
\[
 \sum_{i=1}^{n_{\rm X}}\alpha_i z_i=0,
\]
\noi and:
\[
 \sum_{i=1}^{n_{\rm X}}\alpha_i X_i^{t+\Delta t}=
 \sum_{i=1}^{n_{\rm X}}\alpha_i X_i^t+\sum_{i=1}^{n_{\rm X}}\alpha_i z_i=cte;
\]
\noi therefore, with IRK formulae, the number of baryons and the number of 
charges
are exactly conserved up to the accuracy of the closure of the Newton-Raphson 
scheme.

\subsection{Integration of the diffusion equation}\label{sec:ided}
One seeks the $X_i$ as piecewise polynomials
of order $h_c\ge 2$ {\em with discontinuous first derivative} at each LMR, 
\ie, $X_i\in{\cal C}^0[0,\nu_b]$. For the calculations it is
convenient to project them on their normalized B-spline basis,
${\cal N}\equiv\{N^{h_c}_k\}_{k=1,K_{\rm X}}$
(see \app{sec:ap2}), ${\cal N}\subset{\cal C}^0[0,\nu_b]$, $K_{\rm X}$ 
is the 
dimension of the basis; therefore $X_i(\nu,t)$ is written:
\begin{equation}\label{eq:spli}
X_i(\nu,t)=\sum_{k=1}^{K_{\rm X}} x_{i,k}(t)N^{h_c}_j(\nu).
\end{equation}
One looks for the {\em weak solutions} of \eq{eq:diff} using
the Petrov-Galerkin's formalism\footnote{In 
fact it is a {\em "Tau"} method since the $O^{h_c}_j$ do not
fit the boundary conditions.} (Schatzman \etal{sm}; Marchouk \& Agochkov, 
\cite{mk}; Quarteroni \& Valli, \cite{qv}):
Let ${\cal O}\equiv\{O^{h_c}_k\}_{k=1,K_{\rm X}}$ be
the B-spline basis {\em without discontinuity}, \ie, ${\cal O}\subset{\cal 
C}^{h_c-2}$, constructed on the grid of
$\cal N$ -- schematic representations of nodal
vectors of $\cal N$ and $\cal O$ can be found \app{sec:ga};
for $j=1,\ldots K_X$, multiply both sides of Eq.~(\ref{eq:diff}) with
$O^{h_c}_j$, and integrate over the whole star:
\begin{equation}\label{eq:sca1}
<\frac32M_\odot 
\sqrt\nu\left(\dsdx{X_i}t-\Psi_i\right),O_j^{h_c}>+<\dsdx{F_i}\nu,O_j^{h_c}>=0,
\end{equation}
with the standard notation $<f,g>\equiv\int_0^{\nu_{\rm b}}f(\nu)g(\nu) 
{\rm d}\nu$ for the
inner product of two functions $f$ and $g$.
Due to the continuity of the fluxes $F_i$, the standard integration by 
parts gives:
\begin{equation}\label{eq:inn}
 <\dsdx{F_i}\nu,O_j^{h_c}>=\Big[ F_iO_j^{h_c}\Big]_0^{\nu_{\rm b}}-
 <F_i,\dsdx{O_j^{h_c}}\nu >.
\end{equation}
Owing to the boundary conditions, $F_i(\nu=0)=0$
at the center $\nu\equiv 0$, $F_i(\nu=\nu_{\rm b})=0$ at the outer limit 
$\nu\equiv\nu_{\rm b}$ of the envelope and to the continuity of $X_i$ and $F_i$, 
the set of Eq.~(\ref{eq:sca1}) is written:
\begin{eqnarray}\label{eq:sca2}
 &&<\frac32M_\odot \sqrt\nu\bigg(X_i-X_i^{(t)}-\Psi_i\Delta t\bigg),O_j^{h_c}>+
\nonumber \\
 &&+<\frac{32\pi^2 R^4\rho^2\,\Delta t}{3 \sqrt\nu}
 D_i \bullet \nabla\!_\nu{\cal X},\dsdx{O_j^h}\nu>-\nonumber \\
 &&-<4\pi R^2\rho\Delta t\, V_iX_i,\dsdx{O_j^{h_c}}\nu> = 0,\ j=1,\ldots,K_X,
\end{eqnarray}
here all the quantities are taken at the time $t+\Delta t$ except 
$X_i^{(t)}$ taken at time $t$.
This fully implicit scheme is of order $h_c$ in space and of first order
in time.

One emphasizes the following three facts: (i) in \eq{eq:sca1} the 
discontinuities of the first
derivatives of the $X_i$ are implicitly taken into account by the
integral formulation,
(ii) the use of the B-spline basis is convenient because
the algorithms remain simple even if the locations of the LMR change with 
time,
(iii) the Petrov-Galerkin formulation only involves first order derivatives,
(iv) the fully implicit scheme has no restrictive time step condition
(Richmyer \& Morton, \cite{rm}).

\section{Calculations with B-splines}\label{sec:ap2}
Computation with B-spline is not familiar in stellar modelling,
so I briefly describe the basic of the method directed towards 
interpolation and integration of differential equations. 
For more advanced description and proofs refers to de\,Boor~(\cite{db}), 
Shumaker~(\cite{sk}) and to the literature cited herein.

With this technique the functions to be interpolated \respt{integrated}
\eg pressure, temperature, \etc are approximated by piecewise polynomials
which coincide with the function at interpolation points
\respt{fulfill the differential equation}.
With a lattice with points at $q_i$, for algebra, the simplest way for writing 
piecewise polynomials at the location $q$, is the use of the canonical basis,
\ie, the powers of $(q-q_i)$. For the calculations that
basis leads to algorithms having poor stability conditions; among
the bases available for the linear set of piecewise polynomials, a local basis,
the normalized B-spline basis, leads to stable algorithms.

A B-spline is itself, obviously, a piecewise polynomial.
It can be defined either from divided differences (de\,Boor, \loc;
Gaches, 1988),
or from the algorithm used for their calculation; they are
convolution products of the $\Pi$ (door) function (Marchouk \& Agochkov, 1985).

\subsection{The nodal vector: the key of the computations with B-splines}
\label{sec:nod}
Before going ahead in the definition of the B-splines, one must 
define more precisely what is a piecewise polynomial ${\cal P}(q)$
of order $s$, on a given
partition $\Delta$ of $[1,n]$, $\Delta\equiv \{ 1=q_1 < q_2< \ldots<
q_{n-1} < q_n=n\}$\footnote{of course, the abscissa are not necessarily
equidistant integers.}: ${\cal P}(q)$ is defined as a function 
which coincides, on each sub-interval $[q_i,q_{i+1}],\ 1\leq i<n$,
with a polynomial of degree $s-1$. 
The existence and uniqueness of ${\cal P}(q)$ on $\Delta$
results from the knowledge, at each breakpoint $q_i$, of the kind of connection
which is required between the right and the left pieces of polynomials.
As an example, the so-called "natural spline"
is a piecewise polynomial ${\cal P}_N(q)$ of order 4
(degree 3), constructed in such a way that, at the inner breakpoints
$q_i,\ 1 < i < n$, the second derivatives ${\cal P}^{\prime\prime}_N(q_i)$
are continuous. 

Formally, the rules of connection can differ from one breakpoint
to the next: at some of them, ${\cal P}(q_i)$
can be discontinuous, at some others the left and right pieces
can be tied fulfilling
the continuity of their first derivatives or only of ${\cal P}(q_i)$, \etc
For each $q_i$, the connection is characterized by the so-called
multiplicity $m_i$ obtained as follows:
let $d_i<s,\ 1<i<n$ be the order of the oscularity of the connection
between the left and right polynomial pieces at the breakpoint $q_i$;
$d_i=0$ means ${\cal C}^0$ continuity of
${\cal P}(q_i)$, $d_i=1$ means ${\cal C}^1$ continuity of the first derivative,
$d_i=2$ (${\cal C}^2$) of the second ones \etc and one sets $d_i=-1$ for a
${\cal C}^{-1}$ discontinuity of the function itself.
At $i=1$ and $i=n$ (inner and outer boundaries) one fixes the value
of ${\cal P}(q_i)$. This is equivalent to assume that ${\cal P}$ is
discontinuous at these points. Hence, one has $d_1=d_n=-1$.
The multiplicity of the point $q_i$ is then defined as:
\[
m_i\equiv s-(d_i+1)
\] \noi
From the vector of multiplicities: ${\cal M}=(m_1,\ldots,m_n)^T$
one constructs the heart of the calculation with B-spline:
the so-called "nodal vector": ${\cal N}(\Delta,\ {\cal M},\ s)$
associated to a B-spline basis; it is the vector 
which has the $q_i$ for coordinates, each of them appearing
$m_i$ times. ${\cal N}$ is written:
\begin{eqnarray*}
{\cal N}&=&(\overbrace{q_1,q_1,\ldots,q_1}^{m_1},
              \overbrace{q_2,q_2,\ldots,q_2}^{m_2},\ldots,
              \overbrace{q_n,q_n,\ldots,q_n}^{m_n})^T \\
& \equiv &
({t}_1,{t}_2,{t}_3,\ldots,{t}_{S+m_n})^T.
\end{eqnarray*}
The dimension of $\cal N$ is then $S+m_n$, with $S=\sum_{i=1}^{n-1}
m_i$. 
A convenient schematic representation of a nodal vector consists in
writing one column of $m_i$ elements for each breakpoint $q_i$.

As an example, with $s=4$, $n=7$, 
continuity of ${\cal P}(q)$ at $q_3$, 
of ${\cal P}^\prime(q)$ at $q_6$, 
of ${\cal P}^{\prime\prime}(q)$ at $q_2$ and $q_5$, 
and discontinuity of ${\cal P}(q)$ at $q_4$, 
the nodal vector is represented as:
\[\begin{array}{lllllll}
\times\quad & \times\quad & \times\quad & \times\quad & \times\quad &
\times\quad & \times\quad \\
\times &        & \times & \times &        & \times & \times \\
\times &        & \times & \times &        &        & \times \\
\times & 	&        & \times &	   &	    & \times
\end{array}
\] \noi

Once given this nodal vector, 
the value of the $k$-th normalized B-spline
on ${\cal N}$ ($k=1,2,\ldots,S$) is obtained by induction:
\begin{equation}
 N_k^1(q)=\left\{
 \begin{array}{l}\vspace{0.05 in}
  1,\ {\rm if}\ {t}_k\leq q<{t}_{k+1},\\ \vspace{0.05 in}
  0,\ {\rm otherwise},\\
 \end{array}
 \right.
\end{equation}
\noi and, for $1<\sigma\leq s$:
\begin{equation}
 N_k^\sigma(q)=\frac{q-{t}_k}{{t}_{k+\sigma-1}
 -{t}_k}N_k^{\sigma-1}(q)+
 \frac{{t}_{k+\sigma}-q}{{t}_{k+\sigma}-{t}_{k+1}}N_{k+1}^{\sigma-1}(q).
\end{equation}
Note that its first derivative can easily be calculated:
\[
 \ddx{N_k^s(q)}q=(s-1)\left[\frac{-N_{j+1}^{s-1}(q)}
  {{t}_{k+s}-{t}_{k+1}}
  +\frac{N_j^{s-1}(q)}{{t}_{k+s-1}-{t}_k}\right].
\] \noi
From these relationships one can see that $N_k^s(q)\equiv 0$ if
$q\not\in [{t}_{k-s+1},{t}_{k+s}]$. Hence the B-spline basis
is a local basis. Moreover, $\forall q\in [1,n]$, only $s$ B-splines
of order $s$ are not identically zero and their sum is equal to 1.
For their calculation, algorithms which avoid zero denominator have
been designed. The present work uses the Schumaker's (\cite{sk}) algorithms
5-5, 5-11 and 5-13.

It can be shown that the set ${\cal B}(\Delta,\ {\cal N},\ s)$ of all
the B-splines $N_k^s(q),\ 1\leq k\leq S$, forms a
basis for the linear space of the set of the piecewise polynomials having
${\cal N}$ for nodal vector. Thus, the dimension of ${\cal B}$ is
precisely $S$.

For discontinuous piecewise polynomials the use of B-spline basis is 
particularly 
convenient since only the algorithm which constructs the nodal vector is 
concerned with the discontinuity.

The extension to several dimensions are found in Schumaker (\loc) and in
de\,Boor (\loc)

\subsection{Interpolation}\label{sec:int}
Let $f(q)$ be a real function to be interpolated on the grid
$q_i,\ i=1,\ldots,n$, ($q_i\not= q_j$ if $j\not= i$), at each $q_i$,
$f(q_i)$ is known.
For interpolation of order $s$ the so-called Schoenberg's basis allows to reach
an accuracy almost optimal (de\,Boor, \loc, Chapt. XIII, p. 219). The $n$ 
coordinates of this nodal vector are written:
\begin{eqnarray*}
 && t_i=q_1, i=1,\ldots,s, \\
 && t_{s+i} = \frac{q_{i+1}+\ldots+q_{i+s-1}}{s-1},\ i=1,\ldots,n-s,\\
 &&\ t_i=q_n, i=n,\ldots,n+s.
\end{eqnarray*}
\noi The dimension of this basis is exactly $n$.
As an example, the standard natural spline of interpolation employed this basis, 
with $s=4$, the interpolating function is a cubic piecewise polynomial with 
pieces continuously connected up to second derivative.   

\subsection{Computation of the B-spline coefficients.}\label{sec:bsc}
For $q\in[q_1,q_n]$ the interpolating piecewise polynomial is written:
\[
 {\cal P}(q)=\sum_{i=1}^n p_iN_i^s(q)
\]
\noi where the coefficients $p_i$ fulfill the linear system:
\begin{equation}\label{eq:lin}
 f(q_j)=\sum_{i=1}^n p_iN_i^s(q_j),\ j=1,\ldots,n.
\end{equation}
\noi It is a banded, totally positive linear system, thus it can be accurately 
solved by Gauss elimination without pivoting (de\,Boor, \loc,
Chapt. XIII, p. 203).

\subsection{Piecewise polynomial interpolation is conservative.}\label{sec:cpw}
First briefly recall that standard polynomial interpolation is conservative.
Let $\{f_i\}_{i=1,\ldots,m}$ be a set of real functions
to be interpolated on the grid $\{x_k\}_{k=1,n}$, $(x_k\not=x_j$ if $j\not= k)$.
The polynomial $P(x)$ for interpolation is written (Stoer \& Bulirsch, 
\cite{sb}):
\[
 P(x)=\sum_{j=1}^n f(x_j) L_j(x),
\]
\noi where $L_j(x)$ is the $j$-th. Lagrange's polynomial.
Owing to the property:
\[
\sum_{j=1}^n L_j(x)=1,
\]
\noi if any linear combination holds for the functions, \eg conservation of 
baryon and charge 
numbers are relationships of that form, then:
\[
 \sum_i \alpha_i f_i(x)=Cte,
\]
\noi it also holds for the polynomials of interpolation:
\begin{eqnarray*}
 &&\sum_i \alpha_i P(x)=\sum_i \alpha_i\sum_j f(x_j) L_j(x)=\\
 &&=\sum_j L_j(x)\sum_i\alpha_i f(x_j)=\sum_j L_j(x)\,Cte=Cte.
\end{eqnarray*}
Piecewise polynomials have the same property, since between to break points
(i) they are polynomials and (ii) according 
to the Schoenberg-Whitney theorem (de\,Boor, \loc, p. 200, Theorem XIII.1)
where is at least a location where the linear combination is fulfilled. As seen 
\sec{sec:chim}, the chemicals are interpolated using piecewise polynomials, the 
former property ensures the conservation of baryon and charge numbers.

\subsection{Integration of differential equations}\label{sec:ide}
One seeks, on $[a,b]$, the solution of a differential equation in the form of 
piecewise 
polynomials.
Then, with respect to their B-spline basis, $\{N^s_k(q)\}_{k=1}^S$ of order $s$,
at $q$, for the variables, the piecewise polynomials are written:
\[
 {\cal P}_j(q)=\sum_{k=1}^S p_{k,j}N_k^s(q),
\]
\noi where the $p_{k,j}$ are the coordinates of the $j$-th
unknown with respect to
the B-spline basis. The quantities to be calculated consist in the B-spline
coefficients $p_{k,j}$, rather than in the functions
themselves. Formally, each of the differential equations
can be expressed, for the $p_{k,j}$, in the form:
\begin{equation}\label{eq:dp}
\dsdx{{\cal P}_j(q)}q=\sum_{k=1}^S
 p_{k,j}\dsdx{N_k^s(q)}q={\cal F}_j(q;{\cal P})
\end{equation}
\noi where the function ${\cal F}_j(q;{\cal P})$ is the $j$th. coordinate
in the right hand side of the differential equation.
Hence, derivatives of the functions can be readily calculated by differentiation
of the basis function $N_k^s(q)$.

There are two standard methods for the calculation of the coordinates $p_{k,j}$:
i) in the {\em collocation method} one fulfills Eq.~(\ref{eq:dp}) at a
number of grid points, $q_1,\ldots,q_S$, equal to the dimension $S$ of the 
B-spline basis, ii) in the 
{\em Galerkin method} one defines an inner product and calculates the {\em 
weak solutions} which ensure that the {\em residual}:
\[
\dsdx{{\cal P}_j(q)}q-{\cal F}_j(q;{\cal P})
\]
\noi is orthogonal to the basis:
\begin{equation}\label{eq:gal}
 <\sum_{k=1}^S  p_{k,j}\dsdx{N_k^s(q)}q-{\cal F}_j(q;{\cal P}),N_j^s>=0,
 \ j=1,\ldots, S.
\end{equation}

\noi With non-linear differential problems,
\eq{eq:dp} and \eq{eq:gal} are also non-linear.
They are linearized for the $p_{k,j}$ and solved by iterations
(Newton-Raphson method). For the sake of brevity details are omitted.

\subsubsection{Collocation with the de\,Boor's basis.}\label{sec:dbb}
With the de\,Boor's basis, in the case of first order
differential equations, the continuity of the piecewise polynomials
is ensured at the
mesh points but not their subsequent derivatives. Hence, the nodal
vector is represented by:
\[ \begin{array}{ccccccc}
\times & \times & \times & \dots & \times & \times & \times \\
\times & 	&	 &	 &	  &	   & \times
\end{array}
\mbox{\qquad for $s$=2,}
\] \noi
and
\[ \begin{array}{ccccccc}
\times & \times & \times & \dots & \times & \times & \times \\
\times & \times & \times & \dots & \times & \times & \times \\
\times & 	&	 &	 &	  &	   & \times
\end{array}
\mbox{\qquad for $s$=3.}
\] \noi
This allows to obtain the B-spline functions shown on
\fig{fig:splines}; with $s=2$ the scheme is similar to the centered finite 
difference mid-point scheme. It is clear 
that the continuity of the derivatives of the functions is not
ensured at the mesh points. This does not mean that the functions
derivatives will be discontinuous but rather that they {\em can} be
discontinuous at these points. (This property has been used in the
\sec{sec:set}.) 

\begin{figure}
%\boite{15cm}{2cm}
\defbox{\psfig{figure=splines.ps,%
height=10cm,bbllx=0.5cm,bblly=2cm,bburx=20cm,bbury=24.5cm,clip=}}
\caption{Normalized B-splines of order $s$=2 (left panel) and $s$=3
(right panel) calculated with de\,Boor basis for $n$=3
equidistant mesh points. The nodal vectors are indicated by crosses
below both graphs.}
\label{fig:splines}
\end{figure}
\noi For each mesh and for each piecewise polynomial, $s$ parameters have
to be known (recall that ${\cal P}(q)$ is, on $[q_i,q_{i+1}]$, a
polynomial of degree $s-1$). The continuity of the function
or the boundary conditions of the problem provide one constraint,
the equations have to be written at $s-1$ points within each mesh.
These are called the {\em collocation points}, noted
$(\tau_k)_{k=1}^S$, where $S=(n-1)(s-1)+1$ and are defined as follow:
\begin{equation}
\left\{\begin{array}{l}
-1\le\alpha_1<\alpha_2<\dots<\alpha_{s-1}\le 1, \vspace{.5ex} \\
\tau_{(i-1)s+\sigma}\equiv\dsty\frac 12 \left[q_{i+1}+q_i
+\alpha_\sigma(q_{i+1}-q_i)\right],\\
\qquad\qquad\qquad\qquad i=1,\dots,n;\ \sigma=1,\dots,s.
\end{array}\right.
\end{equation}
The quantities $\alpha_\sigma$ can be almost arbitrarily chosen. However,
for the de\,Boor's nodal vector, they are chosen as the zeros of
the $(s-1)$-th Legendre polynomial.
Note that the collocation points are therefore {\em different from the mesh 
points} $q_i$; this disposition is particularly convenient at the center.
For this particular
collocation pattern, it can be demonstrated that, for first order
differential equations, the overall approximation is of order
O$(h)^s$, but that it is of order up to O$(h)^{2(s-1)}$ at the mesh
points (de\,Boor \& Swartz, \cite{bs}); here $h$ is the size of the mesh. This 
is the so-called {\em superconvergence}. 
 
Once written the differential equation on this particular basis, including
boundary conditions at both ends, is formally written:
\begin{equation}
{\cal F}_j(q; {\bf P})=0,
\label{eq:nl}
\end{equation}
\noi where for the $j$-th unknown ${\cal F}_j$ is a function of $q$ and of
the vector ${\bf P}$ of the coordinates $p_{k,j}$,
$k=1,\ldots,S$.
The resolution of the non-linear system Eq.~(\ref{eq:nl})
is obtained by 
linearization (Newton-Raphson method) starting from the initial model.
Since the B-splines basis is a {\em local}
basis, the jacobian matrix is block-diagonal.
This allows storage ease and faster calculations. 
Here one usually uses either
$s=2$, $s=3$ or $s=4$ for the order of the piecewise polynomials.
As an example for system \eq{eq:syq}, the linear
system to be solved has then one of the following structures for $s=2$ (left)
and $s=3$ (right):
\[
 \left(
 \begin{array}{lllllll}
  \circ &\circ\\
  \Box&\Box\\
  &\Box&\Box\\
  &&\cdot\\
  &&&\cdot\\
  &&&&\cdot\\
  &&&&\Box&\Box\\
  &&&&&\Box&\Box\\
  &&&&&\circ&\circ\\
 \end{array}
 \right)\qquad
\left(
 \begin{array}{lllllllllll}
  \circ &\circ &\circ\\
  \Box&\Box&\Box\\
  \Box&\Box&\Box\\
  &&\Box&\Box&\Box\\
  &&\Box&\Box&\Box\\
  &&&&\cdot\\
  &&&&&\cdot\\
  &&&&&&\cdot\\
  &&&&&&\Box&\Box&\Box\\
  &&&&&&\Box&\Box&\Box\\
  &&&&&&&&\Box&\Box&\Box\\
  &&&&&&&&\Box&\Box&\Box\\
  &&&&&&&&\circ&\circ&\circ\\
 \end{array}
 \right),
\]
\noi here the "$\Box$" are $6\times 6$ matrix for the collocation 
points and the $\circ$ are $3\times6$ matrix for the boundaries.

\subsubsection{Basis for Petrov-Galerkin's method.}\label{sec:ga}
With the Galerkin type methods, the choice of the basis is more flexible.
Here a basis with continuity up to ${\cal C}^{s-2}$ for the derivatives at the
grid points and ${\cal C}^1$ (continuity) or ${\cal C}^0$ (discontinuity) at 
LMR. As examples
the nodal vectors with $s=4$ with continuity at LMR ($\downarrow$) is written: 
\begin{displaymath}
\begin{array}{cccccccccc}
&&&&&\downarrow \\
\times&\times&\times&\dots&\times&\times&\times&\ldots&\times&\times \\
\times&      &      &     &      &\times&      &      &      &\times \\
\times&      &      &     &      &\times&      &      &      &\times \\
\times&      &      &     &      &      &      &      &      &\times
\end{array}
\end{displaymath}
\noi and with discontinuity ($\Downarrow$):

\begin{displaymath}
\begin{array}{cccccccccc}
&&&&&\Downarrow \\
\times&\times&\times&\dots&\times&\times&\times&\ldots&\times&\times \\
\times&      &      &     &      &\times&      &      &      &\times \\
\times&      &      &     &      &\times&      &      &      &\times \\
\times&      &      &     &      &\times&      &      &      &\times
\end{array}
\end{displaymath}
%\noi Note that for parabolic equations the Galerkin method is suitable since 
%the 
%difficulty due to discontinuities is avoided as seen \app{sec:ided}.

 
As an example from \eq{eq:syq} results a jacobian matrix $2s+1$ bloc-diagonal
of $6\times 6$ blocks.
For s=2, the linear
system to be solved has the following structure:
\[
 \left(
 \begin{array}{cccccccc}
  \circ\\
  \Box&\Box\\
  \Box&\Box&\Box\\
  &\Box&\Box&\Box\\
  &&&\cdot\\
  &&&&\Box&\Box&\Box\\
  &&&&&\Box&\Box&\Box\\
  &&&&&&\Box&\Box\\
  &&&&&&&\circ\\
 \end{array}
 \right).
\]
\noi Note that the inner products between two B-splines of
order less or equal to
$s$ are obtained without error with a Gaussian type integration formula with
$s$ points; the same integration formula is also used for all the inner 
products; the calculation of these inner products is conveniently made with 
the algorithm 5-22 of Schumaker (\loc).
Emphasize is made again on the fact that the values of the functions 
are not needed at the grid points, particularly at the center. 

\subsubsection{Solving the linearized systems}\label{sec:sol}
The resolution of the linear system is performed with partial
pivoting. The robustness of the solution is improved by using the modified
Newton-Raphson method which is based on the fact that, for unconstrained 
optimization problems, the "Newton direction is 
a descent direction" (Conte \& de\, Boor, \cite{cdb}, p. 219), giving the 
values of
the unknown functions at iteration $l+1$ from their values at
iteration $l$:
\begin{equation}
p_{k,j}^{(l+1)} = p_{k,j}^{(l)} - \alpha\, dp_{k,j}^{(l)}
\qquad k=1,\dots,S.
\label{eq:NR}
\end{equation}
where $dp_{k,j}^{(l)}$ are the corrections and $\alpha\leq 1$ is a parameter
which ensures that the variations of 
the unknowns $p_{k,j}$ remain small even when $dp_{k,j}^{(l)}$ is large; here 
$\alpha$ is adjusted in order that $|p_{k,j}^{(l)}|\times 0.1 < \alpha\,
\ |dp_{k,j}^{(l)}|$
The iteration then stops when a sufficient accuracy is reached.



