1209.sagews

f = lambda x: x^2+2*x
x0=var('x0')
x0=2

P=[x0,f(x0)]
Qh=[[x0+1/n,f(x0+1/n),1/n] for n in range(1,10)]
print 'P(',P[0],',',P[1],')'
for k in range(9):
    print "h=",Qh[k][2],": Qh(",Qh[k][0],",",Qh[k][1],")"

P( 2 , 8 ) h= 1 : Qh( 3 , 15 ) h= 1/2 : Qh( 5/2 , 45/4 ) h= 1/3 : Qh( 7/3 , 91/9 ) h= 1/4 : Qh( 9/4 , 153/16 ) h= 1/5 : Qh( 11/5 , 231/25 ) h= 1/6 : Qh( 13/6 , 325/36 ) h= 1/7 : Qh( 15/7 , 435/49 ) h= 1/8 : Qh( 17/8 , 561/64 ) h= 1/9 : Qh( 19/9 , 703/81 )

dx=[Qh[k][2] for k in range(9)]
df=[Qh[k][1]-P[1] for k in range(9)]
print dx
print df

[1, 1/2, 1/3, 1/4, 1/5, 1/6, 1/7, 1/8, 1/9] [7, 13/4, 19/9, 25/16, 31/25, 37/36, 43/49, 49/64, 55/81]

M=[df[k]/dx[k] for k in range(8)]
print M

[7, 13/2, 19/3, 25/4, 31/5, 37/6, 43/7, 49/8]

h=var('h')
ra=(f(x0+h)-f(x0))/h
limit(ra,h=0)

6

x=var('x')
Gf=plot(f(x),xmin=1,xmax=4,ymin=7,ymax=16)
R=[]
for k in range(8):
    R.append(plot(M[k]*(x-x0)+f(x0), hue=k/10,xmin=1,xmax=4,ymin=7,ymax=16)+Gf)
a=animate(R,xmin=1,ymin=7,ymax=16)
show(a)