- Email: [email protected]
Linear iterative solvers based on perturbation techniques
C. R. Acad. Sci. Paris, t. 329, Série II b, p. 457–462, 2001
Linear iterative solvers based on perturbation techniques Jean-Marc CADOU a , Naima MOUSTAGHFIR b , El Hassan MALLIL b , Noureddine DAMIL b , Michel POTIER-FERRY a a
Laboratoire de physique et mécanique des matériaux, I.S.G.M.P., université de Metz, Ile du Saulcy, 57045 Metz, France b Laboratoire de calcul scientifique en mécanique, faculté des sciences Ben M’Sik, université Hassan II – Mohammedia, Casablanca, Maroc E-mail: [email protected] (Reçu le 10 février 2001, accepté le 9 avril 2001)
Abstract.
In this paper, we propose a new class of algorithms to solve large scale linear algebraic equations. This method is based on homotopy, perturbation technique and Padé approximants. 2001 Académie des sciences/Éditions scientifiques et médicales Elsevier SAS computational solid mechanics / solids and structures
Des solveurs itératifs basés sur les techniques de perturbation Résumé.
Dans cette Note nous proposons une nouvelle classe d’algorithmes pour résoudre les systèmes linéaires de grande taille. Cette méthode est basée sur des techniques d’homotopie et de perturbation et sur les approximants de Padé. 2001 Académie des sciences/Éditions scientifiques et médicales Elsevier SAS mécanique des solides numérique / solides et structures
Version française abrégée Cette Note concerne la solution de l’équation linéaire (1), où A est une matrice régulière n × n, b ∈ Rn et x ∈ Rn est l’inconnue. En général la solution de (1) est calculée soit par une méthode directe soit par une méthode itérative [1]. L’objet de ce travail est de proposer une variante à ces méthodes classiques. L’algorithme proposé repose sur trois techniques de base : homotopie, perturbation et approximants de Padé. Premièrement, on déforme le problème (1) en introduisant un paramètre et une matrice A∗ choisie arbitrairement, comme dans (2). De cette manière, la solution y() de (2) passe de manière continue de la solution de (3) pour = 0 à la solution de (1) pour = 1. La seule restriction porte sur la régularité de A∗ . Remarquons que la solution de (2), y() = A−1 · b, est une fraction rationnelle en . Deuxièmement y() est cherché sous la forme d’une série tronquée à l’ordre p comme dans (4). Finalement l’approximation polynomiale (4) est remplacée par une approximation rationnelle [2]. Une version itérative de cet algorithme s’obtient en cherchant la solution de (5) ou (6) sous la forme (4). Les termes de la série (4) s’obtiennent en résolvant les problèmes (8) pour k allant de 0 jusqu’à p. La matrice S est similaire à la matrice d’itérations introduite dans les méthodes itératives classiques [3]. En posant = 0 dans (7), on retrouve les méthodes de Jacobi et de Gauss–Seydel selon le choix de la matrice Note présentée par Georges D UVAUT. S1620-7742(01)01357-5/FLA 2001 Académie des sciences/Éditions scientifiques et médicales Elsevier SAS. Tous droits réservés
457
J.M. Cadou et al.
de préconditionnement A∗ . Pour obtenir une représentation rationnelle de la série vectorielle, nous avons suivi la procédure détaillée dans [5]. Dans cette Note, trois exemples de mécanique des structures minces et élastiques sont considérés pour tester les performances de l’algorithme. Le premier exemple est un calcul de plaque carrée encastrée sur les bords et chargée en son centre. Le tableau 1 présente les approximations obtenues par les représentations polynomiales et rationnelles pour différents ordres de troncature. La représentation polynomiale converge vers la solution exacte (calculée par la méthode de Crout) mais plus lentement que l’approximant de Padé correspondant. On considère ensuite deux cas de coques minces : un toît chargé en son centre (exemple 2) et un cylindre troué (exemple 3). On présente sur le tableau 2, l’ordre des séries nécessaire pour converger avec l’algorithme proposé (PM) et le nombre d’itérations avec le gradient conjugué préconditionné (PCG) (l’erreur admissible est 10−10 ). Ce tableau montre que les deux algorithmes sont très sensibles au choix de A∗ : on considère ici des décompositions incomplètes de Choleski. Plus le niveau du préconditionneur est élevé, meilleure est la convergence des deux algorithmes. Par exemple, pour le niveau 3 et pour l’exemple 2 (respectivement l’exemple 3), l’algorithme PM converge pour un ordre de troncature égal à 16 (respectivement 69) alors que l’algorithme PCG converge en 11 itérations (respectivement 32). Ces exemples, même si ils ne sont pas représentatifs des problèmes à très grand nombre d’inconnues, montrent que l’algorithme proposé converge toujours, généralement deux fois moins vite que le PCG. D’autres informations peuvent être obtenues grâce à l’algorithme PM, à partir des racines du dénominateur Q(), qui sont les pôles des approximations rationnelles. En effet, si l’approximation rationnelle était exacte, ces pôles seraient exactement les racines de det(A ) = 0, c’est-à-dire ses valeurs singulières s caractérisées par ∃x = 0 tel que As · x = 0. Or ces valeurs singulières peuvent être reliées aux valeurs propres λs = −(1 − s )/s du problème A · x = λs A∗ · x. Dans les exemples ci-dessus, A∗ et A sont symétriques définies positives. Donc les valeurs singulières sont réelles. De plus A est aussi symétrique définie positive pour entre 0 et 1, donc les valeurs singulières s ne sont pas dans cet intervalle. Or l’algorithme de perturbation donne beaucoup de pôles complexes et quelques uns ont un module entre 0 et 1 : cela signifie que les pôles ne sont que des approximations des valeurs singulières, mais nous avons observé que ces approximations sont cohérentes lorsque est assez petit et que l’ordre de troncature est assez grand. Un préconditionneur efficace sous-entend que A∗ est proche de A, c’est-à-dire que les valeurs propres λs sont voisines de 1. Dans ce cas, les valeurs singulières s sont grandes devant 1 et la technique de perturbation doit converger rapidement. C’est bien ce qui a été observé dans nos tests numériques, qui sont rapportés au tableau 3. Dans ce tableau, nous avons reporté tous les pôles complexes de modules voisins de 1 et tous les pôles réels proches de 1, pour l’ensemble des itérés : pour des niveaux élevés de IC, il n’y a pas de pôles proches de 1 ; au contraire pour les petits niveaux le nombre de ces valeurs critiques est élevé. Cette technique semble donc être un bon outil pour tester la qualité d’un préconditionneur. On a ensuite fixé l’ordre de troncature et effectué plusieurs itérations avec l’algorithme PM. On a reporté sur le tableau 4 le nombre d’itérations nécessaires à la convergence de l’algorithme (pour l’exemple 3, et pour plusieurs choix de A∗ ). Il faut noter que l’algorithme PM tronqué à l’ordre p et avec n itérations conduit à p × n résolutions de problèmes linéaires demandant un seul traitement de A∗ et le calcul de p × n second membres. Dans tous les cas, il paraît préférable d’utiliser un processus itératif avec un ordre de troncature entre 30 et 60. Dans cet article, nous avons montré que l’association d’une homotopie et d’une technique de perturbation peut conduire à un algorithme efficace de résolution des systèmes linéaires. D’après ces premiers tests, le nombre de produits matrice-vecteurs pour obtenir la convergence avec cette méthode est du même ordre qu’avec les méthodes itératives de référence. L’un des principaux intérêts de cet algorithme est qu’il permet d’estimer rapidement la qualité d’un préconditionneur et la vitesse de convergence de l’algorithme, par l’analyse des pôles des approximants de Padé.
458
Linear iterative solvers based on perturbation techniques
1. Introduction This article is concerned with the solution of the following linear equation: A·x=b
(1)
where A is a regular n × n matrix, b ∈ Rn and x ∈ Rn is the unknown. Generally the solution of equation (1) is performed by using either a direct method or an iterative one [1]. The aim of this work is to propose a variant to these classical methods. The proposed algorithm stands on three basic techniques: homotopy, perturbation and Padé approximants. First, we deform the initial linear problem (1) by introducing a homotopy parameter and an arbitrarily chosen matrix A∗ as follows: A · y = (1 − )A∗ + A · y = b
(2)
By this way, the solution y() of equation (2) passes continuously from the solution of the problem: A∗ · y = b
(3)
for = 0 to the solution of the initial problem (1) for = 1. A priori the only restriction is that A∗ must be a regular matrix. Let us remark that the solution of the modified problem (2), y() = A−1 · b, is a rational fraction with respect to . Second y() is searched by using a perturbation technique. So the solution is represented in the form of truncated series with respect to : y() = y0 + y1 + 2 y2 + 3 y3 + · · · + p yp
(4)
Third this polynomial approximation is transformed into a rational one by the classical technique of Padé approximants [2]. In what follows, an iterative process is tested by using a similar technique to compute a sequence of approximated solutions. 2. An iterative algorithm based on perturbation technique Let xi be a trial solution of (1) at the iteration i. Thus the increment xi+1 − xi satisfies: A · xi+1 − xi = −A · xi + b
(5)
In order to compute the next iterate xi+1 from the equation (5), we introduce another homotopy transformation: (6) A∗ · y − xi + A − A∗ · y − xi = −A · xi + b The problem (6) can be rewritten under the form: −1 y − S · y − xi = S · xi + A∗ b,
−1 S = I − A∗ ·A
(7)
The matrix S is similar to the iteration matrix introduced in classical iterative methods [3]. Let us remark that iterative schemes (Jacobi if A∗ is the diagonal of A or Gauss–Seidel if A∗ is the lower or the upper part of A) are deduced from (7) when the parameter is zero. The equation (6) is solved by perturbation technique. By inserting the series (4) into (6), we get linear equations for the vectors yi : A∗ · y0 = −A · xi + b,
A∗ · yk = A∗ − A · yk−1
for k = 1, . . . , p
(8)
459
J.M. Cadou et al.
The same iterative process, based on homotopy and perturbation, has been recently presented in [4] to solve nonlinear partial differential equations. Thus the present algorithm is nothing but the linear version of a technique which proved efficient, especially when it is coupled with Padé approximants. There are many ways to deduce a rational approximation from a series [2]. Here we shall use exactly the same one as in many previous works, see for instance [5]. First the set of vectors (y1 , . . . , yp ) is orthonormalized by a Gram–Schmidt procedure. In this manner, some polynoms fi () appear as coordinates of the orthonormalised basis. Next these polynoms are replaced by rational fractions Pi ()/Q() of degrees [p − 1 − i, p − 1]. Note that all these fractions have the same denominator Q(). The last vector of the orthonormalized basis yp∗ is dropped. We refer to [5] for the details of the algorithm, for a more complete discussion and many practical applications. 3. Numerical results and discussion We first consider a clamped square plate, loaded at its center. We used 32 triangular DKT elements, this leads to a mesh with 150 degrees of freedom (d.o.f.). First we look at the convergence of the proposed technique, in the case where the matrix A∗ is the diagonal part of A. We discuss the convergence of the transverse displacement at the center, the exact value being 0.35348011963248 (computed with the help of the Crout method). The approximations obtained by the series representation and by the Padé representation are reported in table 1. The series converge to the exact solution, but rather slowly as compared with the Padé representation. In the following examples, the proposed method (PM) will always be the one of the Padé approximants. Now we compare the proposed method (PM) with a traditional iterative one: the preconditioned conjugate gradient (PCG). The matrix (A∗ )−1 is computed with an Incomplete Choleski Triangulation (IC) [6] and it is the so-called preconditioning matrix. The level of IC is directly connected to the number of non-zero terms in the matrix A∗ . We refer to [6] for a more general presentation of this algorithm. Table 1. Example 1, comparisons of solutions by series and Padé, with A∗ = diag(A). Tableau 1. Exemple 1, comparaisons des solutions par séries et Padé. Order
Polynomial approximation
Padé approximant
6
0.19827490206583
0.35348449255070
10
0.24510814546692
0.35348011963248
20
0.30972533984282
–
50
0.35067085337659
–
Table 2. Number of computed vectors with the proposed method and the preconditioned conjugate gradient, according to the preconditioner. Tableau 2. Nombre de vecteurs calculés avec la methode proposée et le gradient conjugué préconditionné, selon le préconditionneur. Level of IC Example 2: PM
422
1
2
3
4
5
6
87
57
16
10
7
4
64
29
22
11
6
4
3
Example 3: PM
2866
1275
189
69
24
12
9
Example 3: PCG
199
106
55
32
12
9
6
Example 2: PCG
460
0
Linear iterative solvers based on perturbation techniques
We consider two numerical tests. The first one is a cylinder shell loaded by a single force (example 2) [4]. The second test concerns a cylindrical shell with two diametrically opposed cutouts subjected to a uniform axial compression (example 3) [5]. The element used for both tests is the classical triangular DKT element. The mesh of example 2 is made of 200 elements (726 d.o.f.), while the mesh corresponding to example 3 is composed of 1608 elements (5190 d.o.f.). In table 2, we compare the number of iterations of the PCG algorithm with the order of truncature of the PM algorithm to converge to the solution. We consider here that the methods have converged when the ratio rk /r0 is lower than 10−10 , · is the Euclidian norm, rk and r0 are respectively the residual vector at iteration k and 0 (for the PCG) or order p and 0 (for the PM). In table 2, one can see that the two methods have the same behaviour: the higher is the level of IC, the better is the convergence of the two methods. For example, for the level 3 and example 2 (respectively example 3), the PM algorithm converges for a truncature order 16 (respectively 69) while the PCG scheme needs 11 iterations (respectively 32). According to these two examples, that are representative of moderately large problems, the PM algorithms always converge, but two times less quickly than the PCG. The convergence of PM depends greatly on the preconditioner as with any iterative method. The PM algorithm yields further informations from the roots of the denominator Q(), that are the poles of the rational approximation. Indeed, if the rational approximation was exact, these poles would coincide exactly with the zeros of the determinant of A , i.e. to its singular values s . The latter are characterized by: ∃x ∈ Rn ,
As · x = 0
(9)
These singular values are connected with the eigenvalues λs = −(1 − s )/s of the following problem: A · x = λs A∗ · x
(10)
In the above examples, A∗ and A are symmetric positive definite. Hence the singular values s are real. Further A is also symmetric positive definite for ∈ [0, 1], thus the values s do not lie in this interval. In fact, we have found many complex poles and a few ones in [0, 1]. This contradiction is due to the fact that the poles are only approximations of the singular values. Nevertheless the smallest poles remain more or less constant, whatever be the order, and can be considered as a reliable information. An efficient preconditioner implies that A∗ is close to A, this leads to λs ≈ 1. In this case the singular values s are far from 1 and the perturbation technique must converge rapidly. This is what we have observed in the numerical tests reported in table 3. In this table, we have reported all the complex poles with a modulus close to 1 and all the real poles close to 1: for large levels of IC, there are no poles closed to 1 but there are a lot of poles in the critical area for lower levels. Thus this technique seems to be a good tool to evaluate the efficiency of a preconditioner. The previous tests enabled us to evaluate the performance of the PM algortihm in the case of one iteration. In the following, we fix the truncature order and we report the number of iterations at convergence (table 4). Let us note that the PM algorithm truncated at order p with n iterations needs the resolution of n × p Table 3. Number of Padé’s poles in the critical area. Tableau 3. Nombre de pôles des Padé dont le module est dans la zone critique. Level of IC Example 2: real roots
0
1
137
23
2 9
3
4
5
6
2
0
0
0
2327
42
69
0
0
0
0
Example 3: real roots
–
–
62
18
3
0
0
Example 3: complex roots
–
–
412
62
3
0
0
Example 2: complex roots
461
J.M. Cadou et al. Table 4. Example 3, number of iterations and computed right-hand sides (in parenthesis). Tableau 4. Exemple 3, nombre d’itérations et de seconds membres calculés (entre parenthèses). Order
10
20
30
40
50
60
∞
IC(0)
303 (3021)
140 (2781)
81 (2404)
46 (1826)
34 (1651)
29 (1681)
1 (2866)
IC(1)
121 (1203)
56 (1101)
49 (1467)
26 (1001)
18 (894)
17 (1010)
1 (1275)
IC(2)
36 (351)
13 (242)
6 (179)
6 (224)
4 (163)
4 (224)
1 (189)
IC(3)
13 (121)
5 (84)
3 (72)
2 (63)
2 (62)
2 (61)
1 (69)
problems which requires one incomplete triangulation of A∗ and the computation of at most n × p righthand sides. From the tests, it seems preferable to perform several iterations with an order of truncature between 30 and 60. 4. Conclusion In this article, it has been established that efficient algorithms to solve linear problems can be built up from a homotopy and a perturbation technique. In the presented numerical tests, the number of required vector-matrix products is of the same order as with the reference iterative algorithms. The new perturbation algorithm quickly yields a clear estimate of the quality of the preconditioner by looking at the poles of the Padé approximants: this is one of the main interest of this algorithm. To our knowledge, this kind of solver had not yet been studied in this way. References [1] Ciarlet P., Introduction à l’analyse numérique matricielle et à l’optimisation, Masson, 1982. [2] Brezinski C., van Iseghem J., Padé approximants, in: Ciarlet P.G., Lions J.L. (Eds.), Handbook of Numerical Analysis, Vol. 3, North-Holland, Amsterdam, 1994. [3] Wesseling P., An Introduction to Multigrid Methods, John Wiley & Sons, 1992. [4] Mallil E., Lahmam H., Damil N., Potier-Ferry M., An iterative process based on homotopy and perturbation techniques, Comput. Meth. Appl. Mech. Engng. 190 (2000) 1845–1858. [5] Elhage-Hussein A., Potier-Ferry M., Damil N., A numerical continuation method based on Padé approximants, Int. J. Solids Struct. 37 (2000) 6981–7001. [6] Joly P., Vidrascu M., Quelques méthodes classiques de résolution de systèmes linéaires, in: Joly P., Vidrascu M. (Eds.), INRIA – Collection Didactique, 1994.
462
Linear iterative solvers based on perturbation techniques Jean-Marc CADOU a , Naima MOUSTAGHFIR b , El Hassan MALLIL b , Noureddine DAMIL b , Michel POTIER-FERRY a a
Laboratoire de physique et mécanique des matériaux, I.S.G.M.P., université de Metz, Ile du Saulcy, 57045 Metz, France b Laboratoire de calcul scientifique en mécanique, faculté des sciences Ben M’Sik, université Hassan II – Mohammedia, Casablanca, Maroc E-mail: [email protected] (Reçu le 10 février 2001, accepté le 9 avril 2001)
Abstract.
In this paper, we propose a new class of algorithms to solve large scale linear algebraic equations. This method is based on homotopy, perturbation technique and Padé approximants. 2001 Académie des sciences/Éditions scientifiques et médicales Elsevier SAS computational solid mechanics / solids and structures
Des solveurs itératifs basés sur les techniques de perturbation Résumé.
Dans cette Note nous proposons une nouvelle classe d’algorithmes pour résoudre les systèmes linéaires de grande taille. Cette méthode est basée sur des techniques d’homotopie et de perturbation et sur les approximants de Padé. 2001 Académie des sciences/Éditions scientifiques et médicales Elsevier SAS mécanique des solides numérique / solides et structures
Version française abrégée Cette Note concerne la solution de l’équation linéaire (1), où A est une matrice régulière n × n, b ∈ Rn et x ∈ Rn est l’inconnue. En général la solution de (1) est calculée soit par une méthode directe soit par une méthode itérative [1]. L’objet de ce travail est de proposer une variante à ces méthodes classiques. L’algorithme proposé repose sur trois techniques de base : homotopie, perturbation et approximants de Padé. Premièrement, on déforme le problème (1) en introduisant un paramètre et une matrice A∗ choisie arbitrairement, comme dans (2). De cette manière, la solution y() de (2) passe de manière continue de la solution de (3) pour = 0 à la solution de (1) pour = 1. La seule restriction porte sur la régularité de A∗ . Remarquons que la solution de (2), y() = A−1 · b, est une fraction rationnelle en . Deuxièmement y() est cherché sous la forme d’une série tronquée à l’ordre p comme dans (4). Finalement l’approximation polynomiale (4) est remplacée par une approximation rationnelle [2]. Une version itérative de cet algorithme s’obtient en cherchant la solution de (5) ou (6) sous la forme (4). Les termes de la série (4) s’obtiennent en résolvant les problèmes (8) pour k allant de 0 jusqu’à p. La matrice S est similaire à la matrice d’itérations introduite dans les méthodes itératives classiques [3]. En posant = 0 dans (7), on retrouve les méthodes de Jacobi et de Gauss–Seydel selon le choix de la matrice Note présentée par Georges D UVAUT. S1620-7742(01)01357-5/FLA 2001 Académie des sciences/Éditions scientifiques et médicales Elsevier SAS. Tous droits réservés
457
J.M. Cadou et al.
de préconditionnement A∗ . Pour obtenir une représentation rationnelle de la série vectorielle, nous avons suivi la procédure détaillée dans [5]. Dans cette Note, trois exemples de mécanique des structures minces et élastiques sont considérés pour tester les performances de l’algorithme. Le premier exemple est un calcul de plaque carrée encastrée sur les bords et chargée en son centre. Le tableau 1 présente les approximations obtenues par les représentations polynomiales et rationnelles pour différents ordres de troncature. La représentation polynomiale converge vers la solution exacte (calculée par la méthode de Crout) mais plus lentement que l’approximant de Padé correspondant. On considère ensuite deux cas de coques minces : un toît chargé en son centre (exemple 2) et un cylindre troué (exemple 3). On présente sur le tableau 2, l’ordre des séries nécessaire pour converger avec l’algorithme proposé (PM) et le nombre d’itérations avec le gradient conjugué préconditionné (PCG) (l’erreur admissible est 10−10 ). Ce tableau montre que les deux algorithmes sont très sensibles au choix de A∗ : on considère ici des décompositions incomplètes de Choleski. Plus le niveau du préconditionneur est élevé, meilleure est la convergence des deux algorithmes. Par exemple, pour le niveau 3 et pour l’exemple 2 (respectivement l’exemple 3), l’algorithme PM converge pour un ordre de troncature égal à 16 (respectivement 69) alors que l’algorithme PCG converge en 11 itérations (respectivement 32). Ces exemples, même si ils ne sont pas représentatifs des problèmes à très grand nombre d’inconnues, montrent que l’algorithme proposé converge toujours, généralement deux fois moins vite que le PCG. D’autres informations peuvent être obtenues grâce à l’algorithme PM, à partir des racines du dénominateur Q(), qui sont les pôles des approximations rationnelles. En effet, si l’approximation rationnelle était exacte, ces pôles seraient exactement les racines de det(A ) = 0, c’est-à-dire ses valeurs singulières s caractérisées par ∃x = 0 tel que As · x = 0. Or ces valeurs singulières peuvent être reliées aux valeurs propres λs = −(1 − s )/s du problème A · x = λs A∗ · x. Dans les exemples ci-dessus, A∗ et A sont symétriques définies positives. Donc les valeurs singulières sont réelles. De plus A est aussi symétrique définie positive pour entre 0 et 1, donc les valeurs singulières s ne sont pas dans cet intervalle. Or l’algorithme de perturbation donne beaucoup de pôles complexes et quelques uns ont un module entre 0 et 1 : cela signifie que les pôles ne sont que des approximations des valeurs singulières, mais nous avons observé que ces approximations sont cohérentes lorsque est assez petit et que l’ordre de troncature est assez grand. Un préconditionneur efficace sous-entend que A∗ est proche de A, c’est-à-dire que les valeurs propres λs sont voisines de 1. Dans ce cas, les valeurs singulières s sont grandes devant 1 et la technique de perturbation doit converger rapidement. C’est bien ce qui a été observé dans nos tests numériques, qui sont rapportés au tableau 3. Dans ce tableau, nous avons reporté tous les pôles complexes de modules voisins de 1 et tous les pôles réels proches de 1, pour l’ensemble des itérés : pour des niveaux élevés de IC, il n’y a pas de pôles proches de 1 ; au contraire pour les petits niveaux le nombre de ces valeurs critiques est élevé. Cette technique semble donc être un bon outil pour tester la qualité d’un préconditionneur. On a ensuite fixé l’ordre de troncature et effectué plusieurs itérations avec l’algorithme PM. On a reporté sur le tableau 4 le nombre d’itérations nécessaires à la convergence de l’algorithme (pour l’exemple 3, et pour plusieurs choix de A∗ ). Il faut noter que l’algorithme PM tronqué à l’ordre p et avec n itérations conduit à p × n résolutions de problèmes linéaires demandant un seul traitement de A∗ et le calcul de p × n second membres. Dans tous les cas, il paraît préférable d’utiliser un processus itératif avec un ordre de troncature entre 30 et 60. Dans cet article, nous avons montré que l’association d’une homotopie et d’une technique de perturbation peut conduire à un algorithme efficace de résolution des systèmes linéaires. D’après ces premiers tests, le nombre de produits matrice-vecteurs pour obtenir la convergence avec cette méthode est du même ordre qu’avec les méthodes itératives de référence. L’un des principaux intérêts de cet algorithme est qu’il permet d’estimer rapidement la qualité d’un préconditionneur et la vitesse de convergence de l’algorithme, par l’analyse des pôles des approximants de Padé.
458
Linear iterative solvers based on perturbation techniques
1. Introduction This article is concerned with the solution of the following linear equation: A·x=b
(1)
where A is a regular n × n matrix, b ∈ Rn and x ∈ Rn is the unknown. Generally the solution of equation (1) is performed by using either a direct method or an iterative one [1]. The aim of this work is to propose a variant to these classical methods. The proposed algorithm stands on three basic techniques: homotopy, perturbation and Padé approximants. First, we deform the initial linear problem (1) by introducing a homotopy parameter and an arbitrarily chosen matrix A∗ as follows: A · y = (1 − )A∗ + A · y = b
(2)
By this way, the solution y() of equation (2) passes continuously from the solution of the problem: A∗ · y = b
(3)
for = 0 to the solution of the initial problem (1) for = 1. A priori the only restriction is that A∗ must be a regular matrix. Let us remark that the solution of the modified problem (2), y() = A−1 · b, is a rational fraction with respect to . Second y() is searched by using a perturbation technique. So the solution is represented in the form of truncated series with respect to : y() = y0 + y1 + 2 y2 + 3 y3 + · · · + p yp
(4)
Third this polynomial approximation is transformed into a rational one by the classical technique of Padé approximants [2]. In what follows, an iterative process is tested by using a similar technique to compute a sequence of approximated solutions. 2. An iterative algorithm based on perturbation technique Let xi be a trial solution of (1) at the iteration i. Thus the increment xi+1 − xi satisfies: A · xi+1 − xi = −A · xi + b
(5)
In order to compute the next iterate xi+1 from the equation (5), we introduce another homotopy transformation: (6) A∗ · y − xi + A − A∗ · y − xi = −A · xi + b The problem (6) can be rewritten under the form: −1 y − S · y − xi = S · xi + A∗ b,
−1 S = I − A∗ ·A
(7)
The matrix S is similar to the iteration matrix introduced in classical iterative methods [3]. Let us remark that iterative schemes (Jacobi if A∗ is the diagonal of A or Gauss–Seidel if A∗ is the lower or the upper part of A) are deduced from (7) when the parameter is zero. The equation (6) is solved by perturbation technique. By inserting the series (4) into (6), we get linear equations for the vectors yi : A∗ · y0 = −A · xi + b,
A∗ · yk = A∗ − A · yk−1
for k = 1, . . . , p
(8)
459
J.M. Cadou et al.
The same iterative process, based on homotopy and perturbation, has been recently presented in [4] to solve nonlinear partial differential equations. Thus the present algorithm is nothing but the linear version of a technique which proved efficient, especially when it is coupled with Padé approximants. There are many ways to deduce a rational approximation from a series [2]. Here we shall use exactly the same one as in many previous works, see for instance [5]. First the set of vectors (y1 , . . . , yp ) is orthonormalized by a Gram–Schmidt procedure. In this manner, some polynoms fi () appear as coordinates of the orthonormalised basis. Next these polynoms are replaced by rational fractions Pi ()/Q() of degrees [p − 1 − i, p − 1]. Note that all these fractions have the same denominator Q(). The last vector of the orthonormalized basis yp∗ is dropped. We refer to [5] for the details of the algorithm, for a more complete discussion and many practical applications. 3. Numerical results and discussion We first consider a clamped square plate, loaded at its center. We used 32 triangular DKT elements, this leads to a mesh with 150 degrees of freedom (d.o.f.). First we look at the convergence of the proposed technique, in the case where the matrix A∗ is the diagonal part of A. We discuss the convergence of the transverse displacement at the center, the exact value being 0.35348011963248 (computed with the help of the Crout method). The approximations obtained by the series representation and by the Padé representation are reported in table 1. The series converge to the exact solution, but rather slowly as compared with the Padé representation. In the following examples, the proposed method (PM) will always be the one of the Padé approximants. Now we compare the proposed method (PM) with a traditional iterative one: the preconditioned conjugate gradient (PCG). The matrix (A∗ )−1 is computed with an Incomplete Choleski Triangulation (IC) [6] and it is the so-called preconditioning matrix. The level of IC is directly connected to the number of non-zero terms in the matrix A∗ . We refer to [6] for a more general presentation of this algorithm. Table 1. Example 1, comparisons of solutions by series and Padé, with A∗ = diag(A). Tableau 1. Exemple 1, comparaisons des solutions par séries et Padé. Order
Polynomial approximation
Padé approximant
6
0.19827490206583
0.35348449255070
10
0.24510814546692
0.35348011963248
20
0.30972533984282
–
50
0.35067085337659
–
Table 2. Number of computed vectors with the proposed method and the preconditioned conjugate gradient, according to the preconditioner. Tableau 2. Nombre de vecteurs calculés avec la methode proposée et le gradient conjugué préconditionné, selon le préconditionneur. Level of IC Example 2: PM
422
1
2
3
4
5
6
87
57
16
10
7
4
64
29
22
11
6
4
3
Example 3: PM
2866
1275
189
69
24
12
9
Example 3: PCG
199
106
55
32
12
9
6
Example 2: PCG
460
0
Linear iterative solvers based on perturbation techniques
We consider two numerical tests. The first one is a cylinder shell loaded by a single force (example 2) [4]. The second test concerns a cylindrical shell with two diametrically opposed cutouts subjected to a uniform axial compression (example 3) [5]. The element used for both tests is the classical triangular DKT element. The mesh of example 2 is made of 200 elements (726 d.o.f.), while the mesh corresponding to example 3 is composed of 1608 elements (5190 d.o.f.). In table 2, we compare the number of iterations of the PCG algorithm with the order of truncature of the PM algorithm to converge to the solution. We consider here that the methods have converged when the ratio rk /r0 is lower than 10−10 , · is the Euclidian norm, rk and r0 are respectively the residual vector at iteration k and 0 (for the PCG) or order p and 0 (for the PM). In table 2, one can see that the two methods have the same behaviour: the higher is the level of IC, the better is the convergence of the two methods. For example, for the level 3 and example 2 (respectively example 3), the PM algorithm converges for a truncature order 16 (respectively 69) while the PCG scheme needs 11 iterations (respectively 32). According to these two examples, that are representative of moderately large problems, the PM algorithms always converge, but two times less quickly than the PCG. The convergence of PM depends greatly on the preconditioner as with any iterative method. The PM algorithm yields further informations from the roots of the denominator Q(), that are the poles of the rational approximation. Indeed, if the rational approximation was exact, these poles would coincide exactly with the zeros of the determinant of A , i.e. to its singular values s . The latter are characterized by: ∃x ∈ Rn ,
As · x = 0
(9)
These singular values are connected with the eigenvalues λs = −(1 − s )/s of the following problem: A · x = λs A∗ · x
(10)
In the above examples, A∗ and A are symmetric positive definite. Hence the singular values s are real. Further A is also symmetric positive definite for ∈ [0, 1], thus the values s do not lie in this interval. In fact, we have found many complex poles and a few ones in [0, 1]. This contradiction is due to the fact that the poles are only approximations of the singular values. Nevertheless the smallest poles remain more or less constant, whatever be the order, and can be considered as a reliable information. An efficient preconditioner implies that A∗ is close to A, this leads to λs ≈ 1. In this case the singular values s are far from 1 and the perturbation technique must converge rapidly. This is what we have observed in the numerical tests reported in table 3. In this table, we have reported all the complex poles with a modulus close to 1 and all the real poles close to 1: for large levels of IC, there are no poles closed to 1 but there are a lot of poles in the critical area for lower levels. Thus this technique seems to be a good tool to evaluate the efficiency of a preconditioner. The previous tests enabled us to evaluate the performance of the PM algortihm in the case of one iteration. In the following, we fix the truncature order and we report the number of iterations at convergence (table 4). Let us note that the PM algorithm truncated at order p with n iterations needs the resolution of n × p Table 3. Number of Padé’s poles in the critical area. Tableau 3. Nombre de pôles des Padé dont le module est dans la zone critique. Level of IC Example 2: real roots
0
1
137
23
2 9
3
4
5
6
2
0
0
0
2327
42
69
0
0
0
0
Example 3: real roots
–
–
62
18
3
0
0
Example 3: complex roots
–
–
412
62
3
0
0
Example 2: complex roots
461
J.M. Cadou et al. Table 4. Example 3, number of iterations and computed right-hand sides (in parenthesis). Tableau 4. Exemple 3, nombre d’itérations et de seconds membres calculés (entre parenthèses). Order
10
20
30
40
50
60
∞
IC(0)
303 (3021)
140 (2781)
81 (2404)
46 (1826)
34 (1651)
29 (1681)
1 (2866)
IC(1)
121 (1203)
56 (1101)
49 (1467)
26 (1001)
18 (894)
17 (1010)
1 (1275)
IC(2)
36 (351)
13 (242)
6 (179)
6 (224)
4 (163)
4 (224)
1 (189)
IC(3)
13 (121)
5 (84)
3 (72)
2 (63)
2 (62)
2 (61)
1 (69)
problems which requires one incomplete triangulation of A∗ and the computation of at most n × p righthand sides. From the tests, it seems preferable to perform several iterations with an order of truncature between 30 and 60. 4. Conclusion In this article, it has been established that efficient algorithms to solve linear problems can be built up from a homotopy and a perturbation technique. In the presented numerical tests, the number of required vector-matrix products is of the same order as with the reference iterative algorithms. The new perturbation algorithm quickly yields a clear estimate of the quality of the preconditioner by looking at the poles of the Padé approximants: this is one of the main interest of this algorithm. To our knowledge, this kind of solver had not yet been studied in this way. References [1] Ciarlet P., Introduction à l’analyse numérique matricielle et à l’optimisation, Masson, 1982. [2] Brezinski C., van Iseghem J., Padé approximants, in: Ciarlet P.G., Lions J.L. (Eds.), Handbook of Numerical Analysis, Vol. 3, North-Holland, Amsterdam, 1994. [3] Wesseling P., An Introduction to Multigrid Methods, John Wiley & Sons, 1992. [4] Mallil E., Lahmam H., Damil N., Potier-Ferry M., An iterative process based on homotopy and perturbation techniques, Comput. Meth. Appl. Mech. Engng. 190 (2000) 1845–1858. [5] Elhage-Hussein A., Potier-Ferry M., Damil N., A numerical continuation method based on Padé approximants, Int. J. Solids Struct. 37 (2000) 6981–7001. [6] Joly P., Vidrascu M., Quelques méthodes classiques de résolution de systèmes linéaires, in: Joly P., Vidrascu M. (Eds.), INRIA – Collection Didactique, 1994.
462