<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Extrapolation to ghost point
We will here use Newton polynomial in order to extrapolate to the ghost point. This method should be convergent as the grid spacing goes to $0$. This can be done by evaluating the polynomial outside the range of interpolating points.
Step1:  Intro
We have the following set of points we want to use to build the polynomial
Step2:  sorted after increasing value of the coordinate $x$. These points takes the following values
Step3:  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
Step4:  The coefficients
The coefficients in the Newton polynomial reads
Step5:  Inserted in the polynomial yields
Step6:  When the ghost point is the first in the $x$ array
When we are using the boundary value in the polynomial
Our points to interpolate from are
Step7:  Evaluating the polynomial in $x_g=x_0-\frac{h}{2}$ ,where $x_0$ is the boundary yields
Step8:  When we are not using the boundary value in the polynomial
Our points to interpolate from are
Step9:  Evaluating the polynomial in $x_g=x_0-h$ yields
Step10:  When the ghost point is the last in the $x$ array
When we are using the boundary value in the polynomial
Our points to interpolate from are
Step11:  Evaluating the polynomial in $x_g=x_3+\frac{h}{2}$, where $x_3$ is the boundary yields
Step12:  When we are not using the boundary value in the polynomial
Our points to interpolate from are
Step13:  Evaluating the polynomial in $x_g=x_3+h$ yields
<ASSISTANT_TASK:>
Python Code:
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()
# 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]
Explanation: Extrapolation to ghost point
We will here use Newton polynomial in order to extrapolate to the ghost point. This method should be convergent as the grid spacing goes to $0$. This can be done by evaluating the polynomial outside the range of interpolating points.
End of explanation
display(positions)
Explanation: Intro
We have the following set of points we want to use to build the polynomial
End of explanation
display(values)
Explanation: sorted after increasing value of the coordinate $x$. These points takes the following values
End of explanation
display(Eq(symbols('p_N')(x), get_polynomial(coeffs, positions, symbols('x'))))
Explanation: 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
The Netwon polynomial
The newton polynomial reads
$$p_N(x) = a_{0} + \sum_{i=1}^{N} a_{i} \prod_{j=0}^{i-1} (x - x_j)$$
our fourth order polynomial thus reads
End of explanation
solvedCoeffs = get_coeff(values, positions)
for nr, coeff in enumerate(solvedCoeffs):
    display(Eq(symbols('a_'+str(nr)), coeff))
Explanation: The coefficients
The coefficients in the Newton polynomial reads
End of explanation
display(Eq(Function('p_N')(x), get_polynomial(solvedCoeffs, positions, x)))
Explanation: Inserted in the polynomial yields
End of explanation
# 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)
Explanation: When the ghost point is the first in the $x$ array
When we are using the boundary value in the polynomial
Our points to interpolate from are
End of explanation
# 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)))
Explanation: Evaluating the polynomial in $x_g=x_0-\frac{h}{2}$ ,where $x_0$ is the boundary yields
End of explanation
# 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)
Explanation: When we are not using the boundary value in the polynomial
Our points to interpolate from are
End of explanation
# 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)))
Explanation: Evaluating the polynomial in $x_g=x_0-h$ yields
End of explanation
# 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)
Explanation: When the ghost point is the last in the $x$ array
When we are using the boundary value in the polynomial
Our points to interpolate from are
End of explanation
# 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)))
Explanation: Evaluating the polynomial in $x_g=x_3+\frac{h}{2}$, where $x_3$ is the boundary yields
End of explanation
# 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)
Explanation: When we are not using the boundary value in the polynomial
Our points to interpolate from are
End of explanation
# 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)))
Explanation: Evaluating the polynomial in $x_g=x_3+h$ yields
End of explanation
<END_TASK>