1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
|
from pylab import *
def eval(a,xpoints,x):
n = len(xpoints) - 1
p = a[n]
for k in range(1,n+1):
p = a[n-k] + (x -xpoints[n-k]) * p
return p
def coef(x,y):
a = copy(y)
m = len(x)
for k in range(1,m):
a[k:m] = (a[k:m] - a[k-1])/(x[k:m]-x[k-1])
return a
x = array([0,1,2,3])
y = array([0,3,14,39])
coef = coef(x,y)
NP = 20
newx = linspace(0,3, NP) # New x-values
newy = zeros(NP)
for i in range(NP):
newy[i] = eval(coef, x, newx[i])
plot(newx, newy,'-x')
plot(x, y,'ro')
show()
|
