JakobCHR.com
 
Quick Navigation:
 
Personal:
 Go to Home
 MS Research
 PhD Research
 Curriculum Vitae

General:
 Linux

Soon to come:
 Matlab
 On-line Stores
 Cycling
 Medicine & Health
 LaTeX
 OOP & C++
 Sony PCM-R500 DAT


next up previous contents index
Next: Interpolation of q', and Up: Implementation of the interface Previous: Implementation of the interface   Contents   Index

Interpolation with cubic splines

A cubic spline is a set of 3rd degree polynomials, s(x), one for each function value given by

\begin{displaymath}
s(x) = \sum\limits_{k=0}^{3} d_{jk}(x-x_{j-1})^k \qquad {\mbox{for}} \qquad
x_{j-1} \le x < x_j \;,\; j=\{1,2,\ldots,M\}
\end{displaymath} (13.3)

which satisfies

\begin{displaymath}
\lim\limits_{x\rightarrow x_i^-} s^{(p)}(x) =
\lim\limits_...
... {\mbox{for}} \qquad
p=\{0,1,2\} \;,\; i=\{1,2,\ldots,(M-1)\}
\end{displaymath} (13.4)

where the superscript (p) indicate the pth derivative of s with respect to x and the superscripts + and - indicate right-hand limits and left-hand limits respectively..

In order to fit the cubic spline to the function values we require that

\begin{displaymath}
s_i = y_i \qquad{\mbox{for}}\qquad i=\{0,1,\ldots,M\}
\end{displaymath} (13.5)

where si is defined by

\begin{displaymath}
s_i \;\hbox{$=$\kern-0.68em\raise1.1ex
\hbox{$\scriptscriptstyle\triangle$}}\;s(x_i)
\end{displaymath} (13.6)

If we furthermore assume that


\begin{eqnarray*}
s''(x_0) &=& 0 \\
s''(x_M) &=& 0
\end{eqnarray*}


$\textstyle \parbox{1.50cm}{\begin{eqnarray}
\end{eqnarray}}$

we can transform the problem of finding the d-values in (13.3) which satisfies (13.4), (13.5) and (13.7) to a cubic spline of the form


\begin{eqnarray*}
s(x_{i-1}+h_i t) &=& (1-t)[s_{i-1} + \frac{1}{6} ((1-t)^2 -
...
...\
& & {\mbox{for}}\qquad i=\{1,2,\ldots,M\}\;,\;0\le t \le 1
\end{eqnarray*}


$\textstyle \parbox{1.50cm}{\begin{eqnarray}
\end{eqnarray}}$

where si'' is defined by

\begin{displaymath}
s_i''\;\hbox{$=$\kern-0.68em\raise1.1ex
\hbox{$\scriptscri...
... = \left.\frac{1}{h_i^2}
\frac{d^2 s}{dt^2}\right\vert _{t=1}
\end{displaymath} (13.7)

and hi is defined by

\begin{displaymath}
h_i \;\hbox{$=$\kern-0.68em\raise1.1ex
\hbox{$\scriptscrip...
...$}}\;x_i - x_{i-1} \qquad {\mbox{for}} \qquad i={1,2,\ldots,M}
\end{displaymath} (13.8)

If we apply the continuity conditions (13.4) to (13.8) we find that the constants si'' have to fulfill


\begin{eqnarray*}
h_i s_{i-1}'' + 2(h_i-h_{i+1})s_i'' + h_{i+1}s_{i+1}'' &=& 6\...
...ht] \\
& & \qquad {\mbox{for}} \qquad i=\{1,2,\ldots,(M-1)\}
\end{eqnarray*}


$\textstyle \parbox{1.50cm}{\begin{eqnarray}
\end{eqnarray}}$

Equation (13.11) can be cast into the following matrix equation

\begin{displaymath}
\hspace{0.2ex}\underline{\underline{A}}{}\hspace{0.15ex} \h...
...space{0.15ex}'' = \hspace{0.2ex}\underline{b}{}\hspace{0.15ex}
\end{displaymath} (13.9)

where $\hspace{0.2ex}\underline{\underline{A}}{}\hspace{0.15ex}$ is a symmetric positive definite tri-diagonal matrix of order (M-1).

Since the tri-diagonal coefficient matrix is positive definite we can employ the efficient double-sweep method13.1 (see [6]) which is based on simple Gaussian elimination in order to solve the equations (13.12).

In conclusion we can state that the interpolation of a set of (M+1) points {xi,yi} with a cubic spline is achieved by solving a simple tri-diagonal system of equation which has an order of (M-1). A cubic spline is implemented as a structure of the three vectors $\hspace{0.2ex}\underline{x}{}\hspace{0.15ex}$, $\hspace{0.2ex}\underline{h}{}\hspace{0.15ex}$, $\hspace{0.2ex}\underline{y}{}\hspace{0.15ex}$ and $\hspace{0.2ex}\underline{s}{}\hspace{0.15ex}''$. When we wish to calculate a function value of the spline, s(x*), for $x^*\in
[x_0;x_M]$ we firstly need to search for the j*13.2 which satisfies

\begin{displaymath}
x^*\in [x_{j^*-1};x_{j^*}]
\end{displaymath} (13.10)

With j* known we calculate the function value, s(x*), by (13.8).

Let us for the moment assume that the interpolated function values yi originate from some function f(x), ie yi = f(xi).

We can now define the interpolation error, EI(x), as

\begin{displaymath}
E_I(x) \;\hbox{$=$\kern-0.68em\raise1.1ex
\hbox{$\scriptscriptstyle\triangle$}}\;f(x) - s(x)
\end{displaymath} (13.11)

The maximum interpolation error associated with the cubic spline interpolation described above can be stated as ([62])

\begin{displaymath}
\max\limits_{x_0\le x \le x_M} \vert E_I(x) \vert \le {\cal...
...{\xi={x_0,x_M}} \vert f''(\xi) - s''(\xi) \vert {\cal
O}(H^2)
\end{displaymath} (13.12)

where H is defined by

\begin{displaymath}
H \;\hbox{$=$\kern-0.68em\raise1.1ex
\hbox{$\scriptscriptstyle\triangle$}}\;\max\limits_{i=\{ 1,2,\ldots,M\} } h_i
\end{displaymath} (13.13)

We can conclude that the interpolation error associated with the cubic spline interpolation in the general case where we cannot guarantee f''(x0) = f''(xM) = 0 is of order ${\cal O}(H^2)$. However, in the case where we have a reasonable large number of points (M>15) the interpolation error will decrease as ${\cal O}(H^4)$ in the part of the range [x0;xM] which does not lye close to the two boundary points {x0,xM}.

Note that in our case the points {xi,yi} on which the interpolating function s(x) is based originate from a numerical solution to a given problem13.3which implies that the points from the start is distorted mostly by discretization errors. It is therefore of major importance that the order of the interpolation error is not smaller than the order of the numerical method which produces the points {xi,yi}. If this requirement is not fulfilled the order of the numerical method would be destroyed by interpolation errors.

As we saw in sections 1.5 and 7.2 the order of the numerical methods employed in the coupled model do not exceed 2 and we can therefore state that interpolation with cubic splines is adequate for our case.


next up previous contents index
Next: Interpolation of q', and Up: Implementation of the interface Previous: Implementation of the interface   Contents   Index  
 
 
 
Revision 2.0, Copyright © 1999-2004 Jakob Christensen
http://www.JakobCHR.com
E-Mail: webmaster@JakobCHR.com
Top Quality
Developed with

Danish
Brain Power
Linux Powered!