Newton Forward Divided Difference Formula

H
{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Newton polynomial can be expressed in a simplified form when $ x_0,x_1,...,x_k $ are arranged consecutively with equal spacing. Introducing the notation $h=x_{i+1} - x_i$ for each $ i = 0, 1,...,k-1$ and $ x = x_0+sh$, the difference $x-x_i$ can be written as $(s-i)h$. So the Newton polynomial becomes $$\\sum_{i=0}^k = {s \\choose i}i!h^i[y_0,...,y_i]. $$This is called the Newton forward divided difference formula."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's check this method for the next function: $$f(x*3^x)$$ with $\\varepsilon = 0.0001$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "f(x 0 ) = -0.3333333333333333\n",
      "f(x 1 ) = -0.28867513459481287\n",
      "f(x 2 ) = 0.0\n",
      "f(x 3 ) = 0.8660254037844386\n",
      "f(x 4 ) = 3.0\n",
      "[-0.33333333  0.0446582   0.24401694  0.33333333  0.35726559]\n",
      "Count of parts in the interpolation formula: 5\n",
      "Value of functions f(x) in -0.25 = -0.06957804087824197\n"
     ]
    }
   ],
   "source": [
    "import numpy as np \n",
    "import math\n",
    "x = np.array([-1.0, -0.5, 0, 0.5, 1.0], float)\n",
    "y = np.array([-0.3333333333333333,-0.28867513459481287, 0.0, 0.8660254037844386, 3.0], float)\n",
    "eps = 0.0001\n",
    "n = len(x)\n",
    "\n",
    "for i in range(n):\n",
    "    a = float(x[i]*(3**x[i]))\n",
    "    print(\"f(x\",i,\") =\",a)\n",
    "\n",
    "def coef(x, y): \n",
    "    x.astype(float) \n",
    "    y.astype(float) \n",
    "    n = len(x) \n",
    "    a = [] \n",
    "    for i in range(n): \n",
    "        a.append(y[i]) \n",
    "    for j in range(1, n): \n",
    "        for i in range(n-1, j-1, -1): \n",
    "            a[i] = float(a[i]-a[i-1])\n",
    "    return np.array(a)  \n",
    "print(coef(x, y))\n",
    "r = float(-0.25)\n",
    "t = float((r - x[0])/(x[1]-x[0]))\n",
    "\n",
    "def tsum(t, n):\n",
    "    sum = float(1.0)\n",
    "    for i in range(0, n, 1):\n",
    "        sum = float(sum*(t-i))\n",
    "    return float(sum)\n",
    "    \n",
    "def Eval(a, x, t): \n",
    "    x.astype(float) \n",
    "    n = len(a) \n",
    "    temp = a[0]\n",
    "    count = 1\n",
    "    for i in range(1, n, 1):\n",
    "        temp += (tsum(t, i-1)*a[i])/math.factorial(i)\n",
    "        count = count + 1\n",
    "    print(\"Count of parts in the interpolation formula:\",count)\n",
    "    return temp\n",
    "result = Eval(coef(x,y),x,t)\n",
    "\n",
    "print(\"Value of functions f(x) in\",r,\"=\",result)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
About this Algorithm

The Newton polynomial can be expressed in a simplified form when $ x_0,x_1,...,x_k $ are arranged consecutively with equal spacing. Introducing the notation $h=x_{i+1} - x_i$ for each $ i = 0, 1,...,k-1$ and $ x = x_0+sh$, the difference $x-x_i$ can be written as $(s-i)h$. So the Newton polynomial becomes $$\sum_{i=0}^k = {s \choose i}i!h^i[y_0,...,y_i]. $$This is called the Newton forward divided difference formula.

Let's check this method for the next function: $$f(x*3^x)$$ with $\varepsilon = 0.0001$

import numpy as np 
import math
x = np.array([-1.0, -0.5, 0, 0.5, 1.0], float)
y = np.array([-0.3333333333333333,-0.28867513459481287, 0.0, 0.8660254037844386, 3.0], float)
eps = 0.0001
n = len(x)

for i in range(n):
    a = float(x[i]*(3**x[i]))
    print("f(x",i,") =",a)

def coef(x, y): 
    x.astype(float) 
    y.astype(float) 
    n = len(x) 
    a = [] 
    for i in range(n): 
        a.append(y[i]) 
    for j in range(1, n): 
        for i in range(n-1, j-1, -1): 
            a[i] = float(a[i]-a[i-1])
    return np.array(a)  
print(coef(x, y))
r = float(-0.25)
t = float((r - x[0])/(x[1]-x[0]))

def tsum(t, n):
    sum = float(1.0)
    for i in range(0, n, 1):
        sum = float(sum*(t-i))
    return float(sum)
    
def Eval(a, x, t): 
    x.astype(float) 
    n = len(a) 
    temp = a[0]
    count = 1
    for i in range(1, n, 1):
        temp += (tsum(t, i-1)*a[i])/math.factorial(i)
        count = count + 1
    print("Count of parts in the interpolation formula:",count)
    return temp
result = Eval(coef(x,y),x,t)

print("Value of functions f(x) in",r,"=",result)
f(x 0 ) = -0.3333333333333333
f(x 1 ) = -0.28867513459481287
f(x 2 ) = 0.0
f(x 3 ) = 0.8660254037844386
f(x 4 ) = 3.0
[-0.33333333  0.0446582   0.24401694  0.33333333  0.35726559]
Count of parts in the interpolation formula: 5
Value of functions f(x) in -0.25 = -0.06957804087824197