NewtonMethod.sagews

x=var("x")
f=(x^2-2).function(x)
ll=-0.5
pic=f.plot(xmin=-5,xmax=5,title="f(x)=x^2-2",figsize=[7,7])
root2=point((sqrt(2),0),color="red",pointsize=50)
root2_text=text("root2",(sqrt(2),ll))
show(pic+root2+root2_text)

bg=f.plot(xmin=1,xmax=2,title="f between 1 and 2",figsize=[7,7]);
show(bg+root2+root2_text)

pts,bds=root_finding(f,1,2,5,True)
bg=f.plot(xmin=1,xmax=2,title="keep finding the middle point",figsize=[7,7])
all_points=point([(pts[i],0) for i in range(len(pts))],color="green",pointsize=50);
all_text=[text("a%s"%i,(pts[i],ll)) for i in range(len(pts))];
show(bg+all_points+sum(all_text))

n=5;
lbd,upd=bds[n-1];
bg=f.plot(xmin=1,xmax=2,title="step %s"%n,figsize=[7,7])
all_points=point([(pts[i],0) for i in range(len(pts))],color="green",pointsize=50);
all_text=[text("a%s"%i,(pts[i],ll)) for i in range(len(pts))];
red_points=point([(lbd,0),(upd,0)], color="red", pointsize=50);
orange_point=point(((lbd+upd)/2,0),color="orange",pointsize=50)
show(bg+all_points+sum(all_text)+red_points+orange_point)
print "Root is between %s and %s"%(N(lbd),N(upd))
print "Try changing n=? in the code"

Root is between 1.37500000000000 and 1.43750000000000 Try changing n=? in the code

n_steps=50;
pts=root_finding(f,1,2,n_steps);
for i in range(n_steps):
    print i, N(pts[i]);
    
print "The real answer is %s. It is getting closer!"%N(sqrt(2));
print "Try changing n_steps=?"

0 1.00000000000000 1 2.00000000000000 2 1.50000000000000 3 1.25000000000000 4 1.37500000000000 5 1.43750000000000 6 1.40625000000000 7 1.42187500000000 8 1.41406250000000 9 1.41796875000000 10 1.41601562500000 11 1.41503906250000 12 1.41455078125000 13 1.41430664062500 14 1.41418457031250 15 1.41424560546875 16 1.41421508789062 17 1.41419982910156 18 1.41420745849609 19 1.41421127319336 20 1.41421318054199 21 1.41421413421631 22 1.41421365737915 23 1.41421341896057 24 1.41421353816986 25 1.41421359777451 26 1.41421356797218 27 1.41421355307102 28 1.41421356052160 29 1.41421356424689 30 1.41421356238425 31 1.41421356145293 32 1.41421356191859 33 1.41421356215142 34 1.41421356226783 35 1.41421356232604 36 1.41421356235514 37 1.41421356236970 38 1.41421356237697 39 1.41421356237333 40 1.41421356237151 41 1.41421356237242 42 1.41421356237288 43 1.41421356237311 44 1.41421356237299 45 1.41421356237305 46 1.41421356237308 47 1.41421356237309 48 1.41421356237310 49 1.41421356237310 The real answer is 1.41421356237310. It is getting closer! Try changing n_steps=?

a=2;
bg=f.plot(xmin=2*sqrt(2)-a,xmax=a,title="Newton's Method",figsize=[7,7])
aa=newton_method(f,a);
newton_points=point([(a,0),(a,f(a)),(aa,0)],color="red",pointsize=50)
v_line=line([(a,0),(a,f(a))],color="green");
t_line=line([(a,f(a)),(aa,0)],color="orange");
show(bg+newton_points+v_line+t_line)
print "Initial point: %s"%a;
print "Tangent line:", tangent_line(f,a);
print "New point: %s=%s"%(aa,N(aa));
print "Root 2=%s"%N(sqrt(2));

Initial point: 2 Tangent line: x |--> 4*x - 6 New point: 3/2=1.50000000000000 Root 2=1.41421356237310

print N(exp(0.1))

1.10517091807565

ff=exp(x).function(x);
bg=ff.plot(xmin=-1,xmax=1,title="Linear approximation of exp^x at x=0",figsize=[7,7]);
tl=tangent_line(ff,0).plot(color="orange");
blue_point=point((0.1,exp(0.1)));
orange_point=point((0.1,1.1),color="orange");
show(bg+tl+blue_point+orange_point)

###Run this cell first to provide required functions

def root_finding(f,a,b,steps,details=False):
    """
    Input:
        f: a function
        a: lower bound
        b: upper bound
        steps: number of steps to run the algorithm
    Output:
        A list of "steps" points approaching the root of f.
        If details==True, return ([points], [(lbd,upd) for each steps]).
    """
    if f(a)*f(b)>0:
        return "Not sure if there is a root"
    lbd=a;
    upd=b;
    bounds=[(a,b)];
    pts=[a,b];
    for i in range(steps):
        mid=(lbd+upd)/2;
        if mid==0:
            return mid;
        if f(lbd)*f(mid)<0:
            upd=mid;
            bounds.append((lbd,upd));
            pts.append(mid);
        elif f(mid)*f(upd)<0:
            lbd=mid;
            bounds.append((lbd,upd));
            pts.append(mid);
    if details:
        return (pts, bounds);
    else:
        return pts;
        
def tangent_line(f,x0):
    fp=f.derivative();
    m=fp(x0);
    return (m*(x-x0)+f(x0)).function(x);

def newton_method(f,a):
    """
    Input:
        f: function
        a: inital x
    Output:
        a new x that is closer to the root by applying Newton's method.
    """
    tl=tangent_line(f,a);
    new_a=solve(tl==0,x,solution_dict=True)[0][x];
    return new_a;