newton.py

import numpy as np
import matplotlib.pyplot as plt
import math

def echant(x0, x1, n):
	dx = (x1-x0)/n
	x = np.arange(x0,x1,dx)
	return (x,np.sin(x))

def delta(x, y):
	if len(x) == 1:
		return y[0]
	else:
		return (delta(x[1:],y[1:]) - delta(x[0:len(x)-1], y[0:len(y)-1])) / (x[len(x)-1] - x[0])

def prod(x0, x):
	if len(x) == 1:
		return x0-x[0]
	else:
		return (x0-x[len(x)-1]) * prod(x0, x[0:len(x)-1])

def newt(deltas, x0, x):
	tot = deltas[0]
	for i in range(1,len(x)):
		tot += deltas[i] * prod(x0, x[0:i])
	return tot

def delta_proc(x, y, d):
	delta = np.array([])
	for i in range(0,len(x)-1-d):
		delta = np.append(delta, (y[i+1] - y[i]) / (x[i+1+d] - x[i]))

	return delta


# x = np.array([ 0, 2, 4, 5, 8, 10])
# y = np.array([-1, 1, 6, 0, 2,  5])
n = 67
(x,y) = echant(0, math.pi, n)
print(x,y)

# deltas = np.array([])
# for i in range(0,len(x)):
# 	deltas = np.append(deltas, delta(x[0:i+1],y[0:i+1]))

# print(deltas)

test = y
deltas_proc = np.array([y[0]])
for d in range(0,len(x)-1):
	print(d)
	test = delta_proc(x, test, d)
	deltas_proc = np.append(deltas_proc, test[0])

print(deltas_proc)

# test_y = np.array([])
# for i in x:
# 	test_y = np.append(test_y, newt(deltas, i, x))

# print(np.sum((y-test_y)**2))
# print(np.sum((deltas-deltas_proc)**2))

plt.plot(x,y,'ko')

dx = math.pi / (n)
xx = np.arange(0, math.pi, dx)
yy = np.array([])
for i in xx:
	yy = np.append(yy, newt(deltas_proc, i, x))

plt.plot(xx,yy,'r-')
plt.grid(True)
plt.show()