ajout newton

python/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()