Geek Answers Handbook

What happened to my GMRES implementation?

I am trying to implement GMRES in a Jupyter Notebook, which (in case you don't know):

enter image description here

This is my code:

import numpy as np

def GMRes(A, b, x0, e, nmax_iter, restart=None):
    r = b - np.asarray(np.dot(A, x0)).reshape(-1)

    x = []
    q = [0] * (nmax_iter)

    x.append(r)

    q[0] = r / np.linalg.norm(r)

    h = np.zeros((nmax_iter + 1, nmax_iter))

    for k in range(nmax_iter):
        y = np.asarray(np.dot(A, q[k])).reshape(-1)

        for j in range(k):
            h[j, k] = np.dot(q[j], y)
            y = y - h[j, k] * q[j]
        h[k + 1, k] = np.linalg.norm(y)
        if (h[k + 1, k] != 0 and k != nmax_iter - 1):
            q[k + 1] = y / h[k + 1, k]

        b = np.zeros(nmax_iter + 1)
        b[0] = np.linalg.norm(r)

        result = np.linalg.lstsq(h, b)[0]

        x.append(np.dot(np.asarray(q).transpose(), result) + x0)

    return x

For me, this should be correct, but when I do:

A = np.matrix('1 1; 3 -4')
b = np.array([3, 2])
x0 = np.array([1, 2])

e = 0
nmax_iter = 5

x = GMRes(A, b, x0, e, nmax_iter)

print(x)

Note . At the moment, edoes nothing.

I get this:

[array([0, 7]), array([ 1.,  2.]), array([ 1.35945946,  0.56216216]), array([ 1.73194463,  0.80759216]), array([ 2.01712479,  0.96133459]), array([ 2.01621042,  0.95180204])]

x[k]should approach (32/7, -11/7), since this is the result, but instead it approaches (2, 1), what am I doing wrong?

+4
source share
2 answers

I think the algorithm gives the correct result.

You are trying to solve Ax = b where:

A = \ begin {bmatrix} 1 & 1 \ 3 & -4 \ end {bmatrix}, b = \ begin {bmatrix} 3 \ 2 \ end {bmatrix}

, , , , .

\ begin {cases} x_1 + x_2 = 3 \ 3x_1-4x_2 = 2 \ end {cases}

, , :

x_1 = 2, x_2 = 1

, .

, GMRES scipy:

import scipy.sparse.linalg as spla
import numpy as np

A = np.matrix('1 1; 3 -4')
b = np.array([3, 2])
x0 = np.array([1, 2])
spla.gmres(A,b,x0)


array([ 2.,  1.])
+6

, , . GMRES A. A n, (n + 1) 'th Arnoldi , . n . , , :

-    for k in range(nmax_iter):
+    for k in range(min(nmax_iter, A.shape[0])):
         y = np.asarray(np.dot(A, q[k])).reshape(-1)

-        for j in range(k):
+        for j in range(k + 1):

:

 [array([ 1.        ,  0.35294118]), array([ 2.,  1.])]

. , , n = 2.

+1

More articles:


All Articles