```
In [1]:
```from IPython.display import display
from dividedDifferences import get_coeff, get_polynomial
from sympy import init_printing
from sympy import symbols, simplify
from sympy import Eq, S, Function
init_printing()

```
In [2]:
```# The values in the points we use for the extrapolation
f0, f1, f2, f3 = symbols('f_0, f_1, f_2, f_3', real=True)
values = [f0, f1, f2, f3]
# Grid spacing
h = symbols('h', real=True)
# The value at the ghost
fg = symbols('f_g', real=True)
# Variables just for illustrative purposes
x = symbols('x')
# Coefficients in Newton polynomial
a0, a1, a2, a3 = symbols('a_0, a_1, a_2, a_3')
coeffs = [a0, a1, a2, a3]
# The points to extrapolate from (when they are unspecified)
x0, x1, x2, x3 = symbols('x_0, x_1, x_2, x_3')
positions = [x0, x1, x2, x3]

We have the following set of points we want to use to build the polynomial

```
In [3]:
```display(positions)

sorted after increasing value of the coordinate $x$. These points takes the following values

```
In [4]:
```display(values)

Our goal is to use these four points to build a Newton polynomial we can use to extrapolate to the ghost point $x_g$.

$x_g$ may be the first point in an array of the $x$ coordinates (that is we can order $x$ as $[x_g, x_0, x_1, x_2, x_3 \ldots]$), or it may be the last point (we can order $x$ as $[\ldots, x_0, x_1, x_2, x_3, x_g]$).

In any case:

- Grid points are equally spaced with a gird space $h$ between them.
- The boundary may be the first or last point in the points we want to create the polynomial from, and is located $\frac{h}{2}$ between the ghost point and first/last inner grid points.
- $x_0$ serves as the reference point

```
In [5]:
```display(Eq(symbols('p_N')(x), get_polynomial(coeffs, positions, symbols('x'))))

```
In [6]:
```solvedCoeffs = get_coeff(values, positions)
for nr, coeff in enumerate(solvedCoeffs):
display(Eq(symbols('a_'+str(nr)), coeff))

Inserted in the polynomial yields

```
In [7]:
```display(Eq(Function('p_N')(x), get_polynomial(solvedCoeffs, positions, x)))

```
In [8]:
```# Specification of the inner points
# Using x_0 as the reference point
x_0 = x0
x_1 = x0 + (S(1)/2)*h
x_2 = x0 + (1+S(1)/2)*h
x_3 = x0 + (2+S(1)/2)*h
specifiedPositions = [x_0, x_1, x_2, x_3]
display(specifiedPositions)

Evaluating the polynomial in $x_g=x_0-\frac{h}{2}$ ,where $x_0$ is the boundary yields

```
In [9]:
```# Specification of the position of the ghost point
x_g = x_0 - h*S(1)/2
# Evaluate the polynomial
p = get_polynomial(get_coeff(values, specifiedPositions), specifiedPositions, x_g)
display(Eq(fg, simplify(p)))

```
In [10]:
```# Specification of the inner points
# Using x_0 as the reference point
x_0 = x0
x_1 = x0 + h
x_2 = x0 + 2*h
x_3 = x0 + 3*h
specifiedPositions = [x_0, x_1, x_2, x_3]
display(specifiedPositions)

Evaluating the polynomial in $x_g=x_0-h$ yields

```
In [11]:
```# Specification of the position of the ghost point
x_g = x_0 - h
# Evaluate the polynomial
p = get_polynomial(get_coeff(values, specifiedPositions), specifiedPositions, x_g)
display(Eq(fg, simplify(p)))

```
In [12]:
```# Specification of the inner points
# Using x_0 as the reference point
x_0 = x0
x_1 = x0 + h
x_2 = x0 + 2*h
x_3 = x0 + (2+S(1)/2)*h
specifiedPositions = [x_0, x_1, x_2, x_3]
display(specifiedPositions)

Evaluating the polynomial in $x_g=x_3+\frac{h}{2}$, where $x_3$ is the boundary yields

```
In [13]:
```# Specification of the position of the ghost point
x_g = x_3 + h*S(1)/2
# Evaluate the polynomial
p = get_polynomial(get_coeff(values, specifiedPositions), specifiedPositions, x_g)
display(Eq(fg, simplify(p)))

```
In [14]:
```# Specification of the points
# Using x_0 as the reference point
x_0 = x0
x_1 = x0 + h
x_2 = x0 + 2*h
x_3 = x0 + 3*h
specifiedPositions = [x_0, x_1, x_2, x_3]
display(specifiedPositions)

Evaluating the polynomial in $x_g=x_3+h$ yields

```
In [15]:
```# Specification of the position of the ghost point
x_g = x_3 + h
# Evaluate the polynomial
p = get_polynomial(get_coeff(values, specifiedPositions), specifiedPositions, x_g)
display(Eq(fg, simplify(p)))