Метод сопряженных градиентов

[Urung]

Попытался реализовать алгоритм описанный у Юсуф Саад в книге итерационные методы для разреженных линейных систем

x1 = x0 − tau0*r0,
to ensure that r1 is orthogonal to r0. This means that the two-term recurrence can be
started with alpha0 = 1, and by setting x−1= 0. Putting these relations and definitions
together gives the following algorithm.
ALGORITHM 6.19 CG – Three-Term Recurrence Variant
1. Compute r0 = b − B*x0. Set x−1= 0 and p0 = 1.
2. For j = 0, 1, . . ., until convergence Do:
3. Compute B*rj and gammaj = (rj,rj)/(B*rj,rj)
4. If (j > 0) compute pj=1/(1-(gammaj/gammaj-1)*(rj,rj)/(rj-1,rj-1)/pj-1)
5. Compute xj+1 = pj *(xj −gammaj*rj)+(1 - pj)*xj-1
6. Compute rj+1 =p j*(rj −gammaj* Arj) + (1 −p j)*rj−1
7. EndDo
This algorithm requires slightly more storage than the standard formulation: in addition
to A, the vectors rj ,Arj , rj−1, xj and xj−1 must be kept. It is possible to avoid
keeping rj−1 by computing the residual rj+1 directly as rj+1 = b − Axj+1 in line 6
of the algorithm, but this would entail an additional matrix-vector product.
6
Вот модуль который реализует данный алгоритм

_author__ = 'kardashevskiy'
import numpy as np
from scipy import sparse
def sg2209(B,b,tol=1.0e-8):
    """
    Solve the linear system Ax=b by conjugate gradient method
    """
    n=len(b)
    x0=np.zeros((n),'float')
    r=-B*x0+b
    for k in range(20):
        z=B*r
        gamma=np.dot(r,r)/np.dot(z,r)
        if k<1:
               r1=r
               x1=x0
               p=1.
               r=r-gamma*z
               x=x0-gamma*r
               gamma1=gamma
        else:
            p=1./(1.-(r*r)/(r1*r1)*(gamma/gamma1)/p)
            x2=x
            x=p*(x-gamma*r)+(1.-p)*x1
            x1=x2
            r2=r
            r=p*(r-gamma*z)+(1-p)*r1
            r1=r2
            gamma1=gamma
        if np.dot(r,r)<tol**2:
            break
    return x,k

вроде как мне кажется все правильно а вот тестовая программа в которой этот модуль должна работать

__author__ = 'kardashevskiy'
import numpy as np
from scipy import sparse
from scipy.sparse import linalg
from sgvas import sgvas
from sg2209 import sg2209
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import time
def f(x, y):
    return 2*np.sin(x)*np.sin(y)
    #return 100000.*h*h
def poisson(n1, n2, h, tau, r):
    V = []
    I = []
    J = []
    b = []
    for j in range(n2):
        for i in range(n1):
            k = i+n1*j
            if i==n1-1 :
                b.append(np.sin(j*h))
                V.append(1.)
                I.append(k)
                J.append(k)
            elif j==n2-1 and i<n2-1:
                b.append(np.sin(i*h))
                V.append(1.)
                I.append(k)
                J.append(k)
            else:
                b.append(h*h*f(i*h+h, j*tau+tau))
                V.append(2.-2*r)
                I.append(k)
                J.append(k)
                if i>0:
                   V.append(-1.)
                   I.append(k)
                   J.append(k-1)
                if i<n1-1:
                    V.append(-1.)
                    I.append(k)
                    J.append(k+1)
                if j>0:
                    V.append(r)
                    I.append(k)
                    J.append(k-n1)
                if j<n2-1:
                    V.append(r)
                    I.append(k)
                    J.append(k+n1)
    return sparse.coo_matrix((V, (I,J))), b
N1 = 50
N2 = 50
l1 = np.pi/2
l2 = np.pi/2
h = l1/N1
tau = l2/N2
r = -1.#4.5*h*h/(tau*tau)
tst = time.clock()
A, b = poisson(N1, N2, h, tau, r)
C=A.toarray()
print C
print len(b)
B = A.tocsr()
#v, info = linalg.cg(B, b, tol=1.e-08)
#v, k v, k = sg5crs(B, b, tol=1.e-08)
v,k = sg2209(B, b, tol=1.e-08)
dt = time.clock() - tst
if k <> 0:
    print "N1 = %i, time solution = %1.3e" % (N1, dt)
    print "N2 = %i" % (k)
fig = plt.figure()
ax = Axes3D(fig)
x = np.linspace(h/2, l1-h/2, N1)
y = np.linspace(tau/2, l2-tau/2, N2)
X,Y = np.meshgrid(x, y)
Z = np.reshape(v, (N1, N2))
ax.plot_surface(X, Y, Z, rstride=2, cstride=2, cmap='spectral')
plt.xlabel('$x$')
plt.ylabel('$t$')
plt.show()

помогите я уже засомневался в книге а так ошибка много чего выходит

[Славян]

Что-то я так и не увидел в алгоритме 6.19 как вычисляется r₁. В п.1 есть расчёт r₀, а в п.6 уже происходят вычисления r₂ и больших.

[Urung]

может быть а я то использовал формулу для j=0 а тогда как

[Славян]

Постарайтесь расставлять знаки препинания, а то так не очень понятно, что именно вас озадачивает.

[Urung]

Вы книгу Ю.Саад Итерационные методы для разреженных линейных систем видели? Там на странице 218 вот этот алгоритм я может что не так понял.

[amk]

Скорее всего книгу вообще никто не видел.
Книга пока ещё свежая, издана в прошлом году. В электронные библиотеки пока никто её не клал. В обычных сейчас какое-то странное отношение к читателям, похоже к ним относятся как к помехе, по крайней мере ходят туда только по большой нужде. А много ли людей купят ненужную им в общем-то книгу только для того, чтобы отвечать на задаваемые по ней вопросы? Думаю "мало" будет слишком оптимистичным ответом.