# Question 2 and 3
%pylab inline
def getA(alpha, c, n):
assert alpha>0.
A = c/alpha*np.eye(n)
i, j = np.indices(A.shape)
A[i==j+1] = -c
A[i==j-1] = -c
return A
def resolutionIterative(A, b, x0, omega, eps, nIterMax):
D = diag(diag(A))
if abs(det(D))<1.e-14:
print("Erreur: Matrice D non inversible.")
exit(-1)
else:
# Préparation des matrices
Dinv = diag(1./diag(A))
# Remarque: L'inversion de D est très facile dans le cas présent
# car D diagonale. Il suffit de prendre l'inverse de coefs
# diagonaux comme ci-dessus. Si D n'avait pas été diagonale,
# il aurait fallu résoudre un système linéaire à chaque itération
# avec la commande solve. Il n'aurait surtout pas fallu calculer
# l'inverse avec la commande inv(D) car elle est plus coûteuse
# et beaucoup moins robuste face aux erreurs d'arrondi
J = dot(Dinv, D - A)
J_omega = (1.-omega)*eye(A.shape[0]) + omega*J
b_omega = omega*dot(Dinv,b)
# Itérations
res = norm(dot(A,x0)-b)
x = x0
nIter = 0
while res>eps and nIter<nIterMax:
nIter+=1
x = dot(J_omega, x) + b_omega
res = norm(dot(A,x)-b)
print(res)
return [x, nIter, res]
alpha = 0.1
c = 1.
n = 10
A = getA(alpha, c, n)
omega = 0.5
eps = 1.e-6
x0 = np.zeros(n)
b = np.ones(n)
nIterMax = 100
[x, nIter, res] = resolutionIterative(A, b, x0, omega, eps, nIterMax)
print('residu = ', res)
print('nIter = ', nIter)
print('x = ', x)