AlgeLin - Métodos Iterativos Modernos: Gradiente e Gradiente Conjugado

conjugate_gradient_method.py

import numpy as np

A = np.array([[5,1,1], [1,4,1], [1,1,3]])
b = np.array([10, 12,12])
x = np.array([0,0,0])

def conjugate_gradient_method(A, b, x=np.array([0,0,0]), n_iter=5):
    d = b - np.dot(A,x)
    r = d
    for i in range(n_iter):
        print('Iteration:', i+1)
        alpha_num = np.dot(r.transpose(), d)
        alpha_den_1 = np.dot(d.transpose(), A)
        alpha_den = np.dot(alpha_den_1, d)
        alpha = alpha_num / alpha_den
        print('alpha = {} / {} = {}'.format(alpha_num, alpha_den, alpha))
        x = x + (alpha * d)
        print('x({}) = {}'.format(i+1, x))
        r = b - np.dot(A, x)
        print('r = {}'.format(r))
        beta_num_1 = np.dot(d.transpose(), A)
        beta_num = - np.dot(beta_num_1, r)
        beta_den_1 = np.dot(d.transpose(), A)
        beta_den = np.dot(beta_den_1, d)
        beta = beta_num / beta_den
        print('beta = {} / {} = {}'.format(beta_num, beta_den, beta))
        d = r + (beta * d)
        print('d = {}'.format(d))
        r_norm = np.linalg.norm(r, 2)
        print('r_norm = {}'.format(r_norm))
        print('===========================')
    return np.round(x, 3)

gradient_method.py

import numpy as np

A = np.array([[5,1,1], [1,4,1], [1,1,3]])
b = np.array([10, 12,12])
x = np.array([0,0,0])

def gradient_method(A, b, x=np.array([0,0,0]), n_iter=5):
    r = b - np.dot(A,x)
    print('First r:', r)
    print('x(0) = {}\n'.format(x))
    for i in range(n_iter):
        print('Iteration:', i+1)
        transpose_matmul_a = np.dot(r.transpose(), A)
        alpha_num = round(np.dot(r.transpose(), r), 3)
        alpha_dem = round(np.dot(transpose_matmul_a, r), 3)
        alpha = round(alpha_num / alpha_dem, 3)
        print('alpha = {} / {} = {}'.format(alpha_num, alpha_dem, alpha))
        x = np.round(x + np.round(alpha * r, 3), 3)
        print('x({}) = {}'.format(i+1, x))
        r = np.round(b - np.round(np.dot(A,x),3), 3)
        print('r = {}'.format(r))
        r_norm = np.round(np.linalg.norm(r, 2), 3)
        print('r_norm = {}'.format(r_norm))
        print('===========================')
    return x