3948.sagews

N_MAX = 10

def norm(set):
    return sqrt(sum([s^2 for s in set]))
    
def vector_to_dict(vect, vars):
    dict = {}
    n = 0
    for var in vars:
        dict[var] = vect[n]
        n = n+1
    return dict

def jacobi_matrix(equations, vars):
    n = len(equations)
    jacobi_list = [[equations[row].diff(vars[col]) for col in [0,..,n-1]] for row in [0,..,n-1]]
    return Matrix(jacobi_list)
    
def n_vector(v):
    for k in range(len(v)):
        v[k] = (v[k]).n()
    return v

def n_matrix(M):
    return Matrix([[(M[row][col]).n() for col in range(len(M.columns()))] for row in range(len(M.rows()))])

def fixpunktiteration(contraction, vars, xn, eps_min):
    n = 0
    convergence_delta = []
    eps = eps_min+1
    while eps > eps_min and n < N_MAX:
        #xn1 = gx.subs(vector_to_dict(xn, vars)).n()
        xn1 = [contraction[k].subs(vector_to_dict(xn, vars)).n() for k in range(len(vars))]
        eps = norm(vector(xn)-vector(xn1))
        xn = xn1
        convergence_delta.append(eps)
        n = n+1
    print "Fixpunkt-Lösung %d Iter.: x = %s" % (n, xn)
    return (xn, n, convergence_delta)

def newtoniteration(equation, vars, xn, eps_min):
    n = 0
    xn = vector(xn)
    convergence_delta = []
    eps = eps_min+1
    J = jacobi_matrix(equation, vars)
    while eps > eps_min and n < N_MAX:
        g = vector([equation[k].subs(vector_to_dict(xn, vars)) for k in range(len(vars))])
        Jg = n_matrix(J.subs(vector_to_dict(xn, vars)))
        xn1=n_vector(xn-Jg.inverse()*g)
        
        xn = xn1;     
        eps=norm(g)
        convergence_delta.append(eps)
        n = n+1
    
    print "Newton-Lösung %d Iter.: x = %s" % (n, xn)
    return (xn,n,convergence_delta)

var('x y z')

eps = 10^(-8)
contraction = [(cos(x)+2*y)/6,(sin(x)+x*y^2)/8]
equation = [cos(x)+2*y-6*x,sin(x)+x*y^2-8*y]

xn = [0,0]

print fixpunktiteration(contraction, [x,y], xn, eps)
print newtoniteration(equation, [x,y], xn, eps)

Fixpunkt-Lösung 10 Iter.: x = [0.171333626027552, 0.0213218191362875] ([0.171333626027552, 0.0213218191362875], 10, [0.166666666666667, 0.0208652215157842, 0.00698120969359110, 0.000906223420507779, 0.000296752961852921, 0.0000413721138978165, 0.0000128892957593713, 1.94526559327172e-6, 5.70409136155520e-7, 9.27036527016171e-8]) Newton-Lösung 4 Iter.: x = (0.171333648176476, 0.0213218141513725) ((0.171333648176476, 0.0213218141513725), 4, [1, 0.0151056334257247, 3.31282967596466e-6, 1.63667795093475e-13])

xn = [1,1]

print "Startwert: %s" % xn
fpi = fixpunktiteration(contraction, [x,y], xn, eps)
ni =newtoniteration(equation, [x,y], xn, eps)

Startwert: [1, 1] Fixpunkt-Lösung 10 Iter.: x = [0.171333757409454, 0.0213218591582174] Newton-Lösung 5 Iter.: x = (0.171333648176476, 0.0213218141513725)

def show_convergence(fpi, ni, xn):
    eco_fpi = [-log(fpi[k]-fpi[k+1]).n() for k in range(len(fpi)-1)]
    eco_ni = [-log(ni[k]-ni[k+1]).n() for k in range(len(ni)-1)]
    
    show(list_plot(eco_fpi)+list_plot(eco_ni, color='red'))
    
show_convergence(fpi[2], ni[2], ni[0])