2250.sagews

print "Newton's method"
print ""

import math

x, s, x_1, x_2 = var('x s x_1 x_2')

f(x_1, x_2) = [x_1 + 2*x_2 - 2, x_1^2 + 4*x_2^2 - 4]
print f(x_1, x_2)
print ""

x = f.diff()[0] * 1.0

while true:
    
    s = f.diff()(x[0], x[1]).solve_right(-f(x[0], x[1]))
    x = x + s
    
    if (x - (x + s))(x_1, x_2) == vector(RR,[0, 0]):
        
        break
        
    else:
        
        print x(x_1, x_2)
        print ""
        
print "------------------------------------------------------------------------"

print "Broyden's method"
print ""

import math

B, x, y, s, x_1, x_2 = var('B x y s x_1 x_2')

f(x_1, x_2) = [x_1 + 2*x_2 - 2, x_1^2 + 4*x_2^2 - 4]
print f(x_1, x_2)
print ""

x = MatrixSpace(SR, 2, 1)(f.diff()[0] * 1.0)
B = MatrixSpace(SR, 2, 2)(f.diff()(x[0][0], x[-1][0]))

while true:
    
    s = MatrixSpace(SR, 2, 1)(B.solve_right(-f(x[0][0], x[-1][0])))
    y = MatrixSpace(SR, 2, 1)(f((x+s)[0][0], (x+s)[-1][0]) - f(x[0][0], x[-1][0]))
    
    if (s.transpose()*s)[0][0] == 0:
        
        break
    
    B = B + ((y - B*s)*s.transpose()) / (s.transpose()*s)[0][0]
    x = x + s
    print x
    print ""

Newton's method (x_1 + 2*x_2 - 2, x_1^2 + 4*x_2^2 - 4) (-0.833333333333333, 1.41666666666667) (-0.189393939393939, 1.09469696969697) (-0.0150791353020652, 1.00753956765103) (-0.000112001278299575, 1.00005600063915) (-6.27144090066707e-9, 1.00000000313572) (-2.39908370136766e-16, 1.00000000000000) (2.04180839713296e-16, 1.00000000000000) (-1.78637652117350e-17, 1.00000000000000) ------------------------------------------------------------------------ Broyden's method (x_1 + 2*x_2 - 2, x_1^2 + 4*x_2^2 - 4) [-0.833333333333333] [ 1.41666666666667] [-0.240599733103069] [ 1.12029986655153] [-0.0652258111196708] [ 1.03261290555984] [-0.00680594105921994] [ 1.00340297052961] [-0.000214245281499170] [ 1.00010712264075] [-7.26520225321093e-7] [ 1.00000036326011] [-7.78184241024192e-11] [ 1.00000000003891] [5.02800997404779e-16] [ 1.00000000000000] [-9.86312613186334e-17] [ 1.00000000000000]

print "Newton's method"
print ""

import math

x, s, x_1, x_2 = var('x s x_1 x_2')

f(x_1, x_2) = [(x_1 + 3)*(x_2^3 - 7) + 18, sin(x_2*e^x_1 - 1)]
print f(x_1, x_2)
print ""

x = f.diff()[0] * 1.0
print x(x_1, x_2)
print ""

while true:
    
    s = f.diff()(x[0], x[1]).solve_right(-f(x[0], x[1]))
    x = x + s
    
    if (x - (x + s))(x_1, x_2) == vector(RR,[0, 0]):
        
        break
        
    else:
        
        print x(x_1, x_2)
        print ""
        
print "------------------------------------------------------------------------"

print "Broyden's method"
print ""

import math

B, x, y, s, x_1, x_2 = var('B x y s x_1 x_2')

f(x_1, x_2) = [(x_1 + 3)*(x_2^3 - 7) + 18, sin(x_2*e^x_1 - 1)]
print f(x_1, x_2)
print ""

x = MatrixSpace(SR, 2, 1)(f.diff()[0] * 1.0)
B = MatrixSpace(SR, 2, 2)(f.diff()(x[0][0], x[-1][0]))

while true:
    
    s = MatrixSpace(SR, 2, 1)(B.solve_right(-f(x[0][0], x[-1][0])))
    y = MatrixSpace(SR, 2, 1)(f((x+s)[0][0], (x+s)[-1][0]) - f(x[0][0], x[-1][0]))
    
    if (s.transpose()*s)[0][0] == 0:
        
        break
    
    B = B + ((y - B*s)*s.transpose()) / (s.transpose()*s)[0][0]
    x = x + s
    print x
    print ""

Newton's method ((x_1 + 3)*(x_2^3 - 7) + 18, sin(x_2*e^x_1 - 1)) (x_2^3 - 7, 3.00000000000000*(x_1 + 3)*x_2^2)