# lagSymHyp.py
# Andrew Davison, ad@coe.psu.ac.th, June 2026
'''
Plots f(x) = 1/x and the Lagrange polynomial for that
function based on the interpolation nodes in xs. The 
polynomial is generated using Sympy's interpolating_poly()
and the result is printed in Algebraic form to stdout in
various forms.

The function and polynomial are plotted between x = 1 to 5
in steps of 1/10 to show how the choice of interpolation
nodes can affect its accuracy by causing it to oscillate.
'''


import math

import sympy as sp
from sympy.polys.specialpolys import interpolating_poly
   # https://docs.sympy.org/latest/modules/polys/reference.html

import matplotlib.pyplot as plt
import faUtils

MAX_X = 5

x = sp.symbols('x')

# Interpolation nodes
xs = [0.3, 0.5, 2, 2.75, 4]   # 0.3, 0.5
ys = [1/v for v in xs]

# Lagrange polynomial
ipoly = interpolating_poly(len(xs), x, X=xs, Y=ys)
print("Poly:")
print(ipoly); print()

p = sp.expand(ipoly)   # multiply out all the terms
print("Expanded:")
print(p); print()
print("Factored:")
print(sp.factor(p)); print()
print("Using Fractions:")
print(sp.nsimplify(p))

# Convert to a Python function
pFunc = sp.lambdify(x, p, "math")
print("\nP(3) approx:", pFunc(3))

# Generate plot data
xVals = [i/10 for i in range(1, MAX_X*10+1)]
hypVals = [1/v for v in xVals]
polyVals = [pFunc(v) for v in xVals]
errVals = [abs(poly - hyp)
                  for poly, hyp in zip(polyVals, hypVals)]

plt.figure(figsize=(8, 5))

# Plot curves
plt.plot(xVals, hypVals, label="1/x")
plt.plot(xVals, polyVals, "--", label="Lagrange")

# Plot interpolation nodes
plt.plot(xs, ys, "o", label="Data")

plt.axhline(0, color="black")
plt.axvline(0, color="black")
plt.xlabel("x")
plt.ylabel("y")
plt.grid(True)
plt.legend()
plt.title("1/x vs Lagrange Approx")
plt.show()