newton_polynomial_interpolation_demo

#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Filename: newton_polynomial_interpolation_demo.py
# Version: 1.0.0
# Author: Jeoi Reqi

"""
Description:
    - This script demonstrates polynomial interpolation using Newton's divided differences method.
    - It calculates coefficients of the interpolating polynomial, evaluates the polynomial at specific points,
      finds the minimum value of a function, and calculates error bounds and errors for given points.

Requirements:
    - Python 3.x
    - NumPy
    - SciPy
    - Matplotlib

Functions:
    - newtdd:
        Computes coefficients of interpolating polynomial using Newton's divided differences method.
    - nest:
        Evaluates the interpolating polynomial at a specific point.
    - f5:
        Defines a function used for finding the minimum value.
    - g:
        Defines a function used for error calculation.
    - h:
        Defines a function used for error calculation.
    - divided_diff:
        Calculates the divided differences table.
    - newton_poly:
        Evaluates the Newton polynomial at specific points.

Usage:
    - Run the script and observe the printed results.

Expected Example Output:

    --------------------------------------------------------
         :: Running The Newton Interpolation Test ::
    --------------------------------------------------------
    [Interpolating Polynomial Coefficients]:
    [1.433329   1.98987    3.2589     3.68066667 4.00041667]

    [Interpolating Polynomial at x = 0.82]:
    1.9589097744

    [Interpolating Polynomial at x = 0.98]:
    2.6128479663999995

    [Minimum value of f5]:
    -848.0804466745346

    --------------------------------------------------------
    Error bounds and errors:
    --------------------------------------------------------
    Ec1: -5.3734377101298326e-05
    Er1: 2.3348514214704963e-05
    Ec2: 0.00021656582286280922
    Er2: 0.00010660542393381434
    --------------------------------------------------------

Additional Notes:
    - Ensure NumPy, SciPy, and Matplotlib are installed in your Python environment.
    - This script replicates functionality from MATLAB code.
"""

import numpy as np
from scipy.optimize import minimize_scalar
import matplotlib.pyplot as plt

def newtdd(x_vals, y_vals, n):
    """
    Calculate coefficients of interpolating polynomial using Newton's divided differences method.

    Parameters:
        x_vals (ndarray): Array of x-coordinates of data points.
        y_vals (ndarray): Array of y-coordinates of data points.
        n (int): Number of data points.

    Returns:
        ndarray: Coefficients of the interpolating polynomial.
    """
    v = np.zeros((n, n))
    for j in range(n):
        v[j][0] = y_vals[j]
    for i in range(1, n):
        for j in range(n - i):
            v[j][i] = (v[j+1][i-1] - v[j][i-1]) / (x_vals[i+j] - x_vals[j])
    return v[0]

def nest(coeff, x_vals, t):
    """
    Evaluate the interpolating polynomial at a specific point.

    Parameters:
        coeff (ndarray): Coefficients of the interpolating polynomial.
        x_vals (ndarray): Array of x-coordinates of data points.
        t (float): Point at which to evaluate the polynomial.

    Returns:
        float: Value of the polynomial at point t.
    """
    n = len(coeff)
    result = coeff[-1]
    for i in range(2, n+1):
        result = result * (t - x_vals[n-i]) + coeff[n-i]
    return result

def f5(x_val):
    """
    Define a function used for finding the minimum value.

    Parameters:
        x_val (float): Input value.

    Returns:
        float: Result of the function.
    """
    return -8 * x_val * (4 * x_val**4 + 20 * x_val**2 + 15) * np.exp(x_val**2)

def g(x_val):
    """
    Define a function used for error calculation.

    Parameters:
        x_val (float): Input value.

    Returns:
        float: Result of the function.
    """
    return np.exp(x_val**2)

def h(x_val):
    """
    Define a function used for error calculation.

    Parameters:
        x_val (float): Input value.

    Returns:
        float: Result of the function.
    """
    return (x_val - 0.6) * (x_val - 0.7) * (x_val - 0.8) * (x_val - 0.9) * (x_val - 1)

def divided_diff(x_vals, y_vals):
    '''
    Function to calculate the divided differences table.

    Parameters:
        x_vals (ndarray): Array of x-coordinates of data points.
        y_vals (ndarray): Array of y-coordinates of data points.

    Returns:
        ndarray: Divided differences table.
    '''
    n = len(y_vals)
    coef = np.zeros([n, n])
    # The first column is y
    coef[:,0] = y_vals

    for j in range(1,n):
        for i in range(n-j):
            coef[i][j] = (coef[i+1][j-1] - coef[i][j-1]) / (x_vals[i+j]-x_vals[i])

    return coef

def newton_poly(coef, x_values, points):
    '''
    Evaluate the Newton polynomial at points.

    Parameters:
        coef (ndarray): Coefficients of the interpolating polynomial.
        x_values (ndarray): Array of x-coordinates of data points.
        points (ndarray): Points at which to evaluate the polynomial.

    Returns:
        ndarray: Value of the polynomial at points.
    '''
    n = len(coef) - 1
    p = coef[n]
    for k in range(1, n + 1):
        p = coef[n-k] + (points - x_values[n-k]) * p
    return p

# Define the data points
x_data = np.array([0.6, 0.7, 0.8, 0.9, 1.0])
y_data = np.array([1.433329, 1.632316, 1.896481, 2.247908, 2.718282])

# Calculate the coefficients of interpolating polynomial
coefficients = newtdd(x_data, y_data, len(x_data))

# Evaluate the interpolating polynomial at specific points
P4_082 = nest(coefficients, x_data, 0.82)
P4_098 = nest(coefficients, x_data, 0.98)

# Find the minimum value of f5
f_5max = minimize_scalar(f5, bounds=(0.6, 1)).fun

# Calculate error bounds and errors for specific points
Ec1 = f_5max * h(0.82) / 120
Er1 = np.abs(P4_082 - g(0.82))
Ec2 = f_5max * h(0.98) / 120
Er2 = np.abs(P4_098 - g(0.98))

# Print the results
print("-" * 56)
print("     :: Running The Newton Interpolation Test ::")
print("-" * 56)
print("[Interpolating Polynomial Coefficients]:")
print(coefficients, "\n")
print("[Interpolating Polynomial at x = 0.82]:")
print(P4_082, "\n")
print("[Interpolating Polynomial at x = 0.98]:")
print(P4_098, "\n")
print("[Minimum value of f5]:")
print(f_5max, "\n")
print("-" * 56)
print("Error bounds and errors:")
print("-" * 56)
print("Ec1:", Ec1)
print("Er1:", Er1)
print("Ec2:", Ec2)
print("Er2:", Er2)
print("-" * 56)

# Test Newton's Polynomial Interpolation
x_test = np.array([-5, -1, 0, 2])
y_test = np.array([-2, 6, 1, 3])
coefficients_test = divided_diff(x_test, y_test)[0, :]
x_new = np.arange(-5, 2.1, 0.1)
y_new = newton_poly(coefficients_test, x_test, x_new)

# Plot the results
plt.figure(figsize=(12, 8))
plt.plot(x_test, y_test, 'bo', label='Data Points')  # Label the data points as 'Data Points'
plt.plot(x_new, y_new, label='Newton Polynomial Interpolation')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Newton Polynomial Interpolation Test')

# Annotate each data point with its corresponding polynomial value
for idx, (x, y) in enumerate(zip(x_test, y_test)):
    plt.text(x, y, f'P{idx+1}: {coefficients_test[idx]:.4f}', fontsize=12, ha='right', va='bottom')

plt.legend()  # Show legend with labels
plt.grid(True)
plt.show()