> Very interesting... but, I think you sent me a copy of my worksheet...
> Yours, very sincerely
While the name is the same, it's a different worksheet. Anyways,
here's the faster code:
def NewtonInterpolation(x,y,f):
poly=f[0]
q=1
s=f[0:len(f)]
stride=1
for k in range(len(f)-1,0,-1):
for i in range(k):
s[i]=(s[i+1]-s[i])/(y[i+stride]-y[i])
q*=(x-y[stride-1])
poly+=s[0]*q
stride+=1
return poly
x = var('x')
fonc = 1/(1+25*x^2)
ffonc = fonc._fast_float_('x')
x = RDF['x'].gen()
pf = plot(fonc,-1,1, thickness=1,color='red')
@interact
def Unif(points=(3..40)):
t = cputime()
pi = RDF.pi()
y=[RDF(-1+2*i/(points-1)) for i in range(points)]
col=point([(y[i],0) for i in range(points)],rgbcolor='red',pointsize=30)
p = NewtonInterpolation(x, y, [ffonc(z) for z in y])
html('$%s$'%latex('Le phenomene de Runge'))
html('Interpolation polynomiale, points de collocation equidistants, et')
html('interpolation polynomiale aux points de Tchebychev.')
html('$f(x)\;=\;%s$'%latex(fonc))
print points,' points.'
pp= plot(p,-1,1, thickness=1)
show(col+pp+pf,ymin=-1,ymax=1.5)
ytche=[cos(i*pi/(points-1)) for i in range(points)]
coltche=point([(ytche[i],0) for i in
range(points)],rgbcolor='red',pointsize=30)
ptche=NewtonInterpolation(x,ytche,[ffonc(z) for z in ytche])
pptche= plot(ptche,-1,1, thickness=1)
show(coltche+pf+pptche,ymin=-1,ymax=1.5)
print cputime()-t
--Mike
--~--~---------~--~----~------------~-------~--~----~
To post to this group, send email to [email protected]
To unsubscribe from this group, send email to [EMAIL PROTECTED]
For more options, visit this group at
http://groups.google.com/group/sage-support
URLs: http://www.sagemath.org
-~----------~----~----~----~------~----~------~--~---
