Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

Chapter Contents
Chapter Introduction
NAG Toolbox

# NAG Toolbox: nag_interp_1d_everett (e01ab)

## Purpose

nag_interp_1d_everett (e01ab) interpolates a function of one variable at a given point $x$ from a table of function values evaluated at equidistant points, using Everett's formula.

## Syntax

[a, g, ifail] = e01ab(n, p, a)
[a, g, ifail] = nag_interp_1d_everett(n, p, a)
Note: the interface to this routine has changed since earlier releases of the toolbox:
 At Mark 23: n1 is no longer an optional input parameter; n2 is no longer an input parameter

## Description

nag_interp_1d_everett (e01ab) interpolates a function of one variable at a given point
 $x=x0+ph,$
where $-1 and $h$ is the interval of differencing, from a table of values ${x}_{m}={x}_{0}+mh$ and ${y}_{m}$ where $m=-\left(n-1\right),-\left(n-2\right),\dots ,-1,0,1,\dots ,n$. The formula used is that of Fröberg (1970), neglecting the remainder term:
 $yp=∑r=0 n-1 1-p+r 2r+1 δ2ry0+∑r=0 n-1 p+r 2r+1 δ2ry1.$
The values of ${\delta }^{2r}{y}_{0}$ and ${\delta }^{2r}{y}_{1}$ are stored on exit from the function in addition to the interpolated function value ${y}_{p}$.

## References

Fröberg C E (1970) Introduction to Numerical Analysis Addison–Wesley

## Parameters

### Compulsory Input Parameters

1:     $\mathrm{n}$int64int32nag_int scalar
$n$, half the number of points to be used in the interpolation.
Constraint: ${\mathbf{n}}>0$.
2:     $\mathrm{p}$ – double scalar
The point $p$ at which the interpolated function value is required, i.e., $p=\left(x-{x}_{0}\right)/h$ with $-1.0.
Constraint: $-1.0<{\mathbf{p}}<1.0$.
3:     $\mathrm{a}\left(\mathit{n1}\right)$ – double array
${\mathbf{a}}\left(\mathit{i}\right)$ must be set to the function value ${y}_{\mathit{i}-n}$, for $\mathit{i}=1,2,\dots ,2n$.

None.

### Output Parameters

1:     $\mathrm{a}\left(\mathit{n1}\right)$ – double array
$\mathit{n1}=2×{\mathbf{n}}$.
The contents of a are unspecified.
2:     $\mathrm{g}\left(\mathit{n2}\right)$ – double array
$\mathit{n2}=2×{\mathbf{n}}+1$.
The array contains
 $\phantom{{\delta }^{2r}}{y}_{0}$ in ${\mathbf{g}}\left(1\right)$ $\phantom{{\delta }^{2r}}{y}_{1}$ in ${\mathbf{g}}\left(2\right)$ ${\delta }^{2r}{y}_{0}$ in ${\mathbf{g}}\left(2r+1\right)$ ${\delta }^{2r}{y}_{1}$ in ${\mathbf{g}}\left(2\mathit{r}+2\right)$, for $\mathit{r}=1,2,\dots ,n-1$.
The interpolated function value ${y}_{p}$ is stored in ${\mathbf{g}}\left(2n+1\right)$.
3:     $\mathrm{ifail}$int64int32nag_int scalar
${\mathbf{ifail}}={\mathbf{0}}$ unless the function detects an error (see Error Indicators and Warnings).

## Error Indicators and Warnings

Errors or warnings detected by the function:
${\mathbf{ifail}}=1$
 On entry, ${\mathbf{p}}\le -1.0$, or ${\mathbf{p}}\ge 1.0$.
${\mathbf{ifail}}=-99$
${\mathbf{ifail}}=-399$
Your licence key may have expired or may not have been installed correctly.
${\mathbf{ifail}}=-999$
Dynamic memory allocation failed.

## Accuracy

In general, increasing $n$ improves the accuracy of the result until full attainable accuracy is reached, after which it might deteriorate. If $x$ lies in the central interval of the data (i.e., $0.0\le p<1.0$), as is desirable, an upper bound on the contribution of the highest order differences (which is usually an upper bound on the error of the result) is given approximately in terms of the elements of the array g by $a×\left(\left|{\mathbf{g}}\left(2n-1\right)\right|+\left|{\mathbf{g}}\left(2n\right)\right|\right)$, where $a=0.1$, $0.02$, $0.005$, $0.001$, $0.0002$ for $n=1,2,3,4,5$ respectively, thereafter decreasing roughly by a factor of $4$ each time.

The computation time increases as the order of $n$ increases.

## Example

This example interpolates at the point $x=0.28$ from the function values
 $xi -1.00 -0.50 0.00 0.50 1.00 1.50 yi 0.00 -0.53 -1.00 -0.46 2.00 11.09 .$
We take $n=3$ and $p=0.56$.
```function e01ab_example

fprintf('e01ab example results\n\n');

a = [-1  -0.50   0    0.50    1    1.50];
b = [ 0  -0.53  -1   -0.46    2   11.09];
n = int64(size(a,2)/2);

x = 0.28;
% We get p from x = a(n) + p*h, where h = 0.5
p = (x-a(n))/0.5;

[bx, c, ifail] = e01ab(n, p, b);

for k = 0:n-1
fprintf('Central differences order %4d of y_0 = %12.5f\n', k, c(2*k+1));
fprintf('%37s = %12.5f\n', 'y_1',c(2*k+2));
end
fprintf('\nFunction value at interpolation point = %12.5f\n', c(end));

```
```e01ab example results

Central differences order    0 of y_0 =     -1.00000
y_1 =     -0.46000
Central differences order    1 of y_0 =      1.01000
y_1 =      1.92000
Central differences order    2 of y_0 =     -0.04000
y_1 =      3.80000

Function value at interpolation point =     -0.83591
```