Наши проекты:
Журнал · Discuz!ML · Wiki · DRKB · Помощь проекту |
||
ПРАВИЛА | FAQ | Помощь | Поиск | Участники | Календарь | Избранное | RSS |
[3.129.13.201] |
|
Сообщ.
#1
,
|
|
|
Попытался реализовать алгоритм описанный у Юсуф Саад в книге итерационные методы для разреженных линейных систем
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() |
Сообщ.
#2
,
|
|
|
Что-то я так и не увидел в алгоритме 6.19 как вычисляется r1. В п.1 есть расчёт r0, а в п.6 уже происходят вычисления r2 и больших.
|
Сообщ.
#3
,
|
|
|
может быть а я то использовал формулу для j=0 а тогда как
|
Сообщ.
#4
,
|
|
|
Постарайтесь расставлять знаки препинания, а то так не очень понятно, что именно вас озадачивает.
|
Сообщ.
#5
,
|
|
|
Вы книгу Ю.Саад Итерационные методы для разреженных линейных систем видели? Там на странице 218 вот этот алгоритм я может что не так понял.
|
Сообщ.
#6
,
|
|
|
Скорее всего книгу вообще никто не видел.
Книга пока ещё свежая, издана в прошлом году. В электронные библиотеки пока никто её не клал. В обычных сейчас какое-то странное отношение к читателям, похоже к ним относятся как к помехе, по крайней мере ходят туда только по большой нужде. А много ли людей купят ненужную им в общем-то книгу только для того, чтобы отвечать на задаваемые по ней вопросы? Думаю "мало" будет слишком оптимистичным ответом. |