CAMA : ma41 Système matriciel non linéaire

Cours du 25 / 06

Si la matrice

$A$ dépend de

$\mathbf{x}$ (noté

$A(\mathbf{x})$), alors le système matriciel

$$A(\mathbf{x}\mathbf{)}\mathbf{x}\mathbf{=}\mathbf{b}$$
n'est pas linéaire.

Exemple :

$$\left[\begin{array}{cc}1& x_1\\ 2x_1& -x_2\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\end{array}\right]=\left[\begin{array}{c}{b}_1\\ b_2\end{array}\right]$$

donne le système suivant non linéaire puisqu'on a des
carrés :

$$\begin{array}{rl}{x}_1+x_1{x}_2& =b_1\\ 2x_1^2-x_2^2& =b_2\end{array}$$

Comment résoudre un tel système ?

La méthode du point fixe

La méthode du point fixe consiste à appliquer l'algorithme itératif suivant :

$${\mathbf{x}}^{k+1}=g({\mathbf{x}}^k)$$
pour résoudre

$g(\mathbf{x})=\mathbf{x}$.

• $x_0$ donné
• $\overline{\mathbf{x}}=g(\overline{\mathbf{x}})$ un point fixe de
  $g$
• ${\mathbf{x}}^{k+1}=g({\mathbf{x}}^k)$ avec
  $k=0,1,2,...$

Est-ce que

$g({\mathbf{x}}^k{)}^k$ converge ?

• Si
  ${\mathbf{x}}_{\mathbf{0}}<{\overline{\mathbf{x}}}_{\mathbf{2}}$ :
  $\underset{k\to +\mathrm{\infty }}{lim}={\overline{\mathbf{x}}}_{\mathbf{1}}$
• Si
  ${\mathbf{x}}_{\mathbf{0}}>{\overline{\mathbf{x}}}_{\mathbf{2}}$ :
  $\underset{k\to +\mathrm{\infty }}{lim}=+\mathrm{\infty }$

Selon le point de départ, la méthode converge ou diverge.

La méthode du point fixe pour résoudre

$A(\mathbf{x}\mathbf{)}\mathbf{x}\mathbf{=}\mathbf{b}$

On doit définir une fonction

$g$ telle que la solution de

$J(\mathbf{x})=\mathbf{x}$ soit la solution du système matriciel non linéaire :

$$g(\mathbf{x})=A^{-1}(\mathbf{x})\mathbf{b}$$

Inverser A est trop coûteux, on écrit notre algorithme itératif sous forme de problème matriciel linéaire à résoudre:

$$A({\mathbf{x}}^k){\mathbf{x}}^{k+1}=\mathbf{b}$$
Si on connait

${\mathbf{x}}^k$ on peut évaluer

$A({\mathbf{x}}^k)$, le système est linéaire et permet de trouver

${\mathbf{x}}^{k+1}$.
Le fonctionnement de l'algorithme va dépendre du type de la matrice

$A$ et de la valeur initiale

$x_0$.

Exemple :

$$\left[\begin{array}{cc}{x}_0-2x_1& x_1\\ x_0& 2x_0+x_1\end{array}\right]\left[\begin{array}{c}{x}_0\\ x_1\end{array}\right]=\left[\begin{array}{c}1\\ 9\end{array}\right]$$
Ce système a pour solutions [1, 2] et [2, 1].

def A(x):
    return np.array([[x[0] - 2*x[1], x[1]] ,
                     [x[0] , x[1] + 2*x[0]]]) / 10 

b = np.array([1, 9]) / 10         # avec normalisation grossière

x = np.array([1, 1])
for i in range(1, 10):
    x = lin.solve(A(x), b)

...
x^7 =  [1.37317932 2.74635363]
x^8 =  [0.72823777 1.45647516]
x^9 =  [1.37317809 2.74635608]

La solution oscille sans converger. La méthode du point fixe a un petit rayon de convergence.

Appliquon l'inertie

µ = 0.5  # on avance de moitié vers le prochain x

x = np.array([3, 2])
for i in range(1, 10):
    x_old = x
    x = lin.solve(A(x), b)
    x = µ * x + (1-µ) * x_old

...
x^9 =  [1.00876429 1.98253948]

La convergence est rapide (9 iterations). L'inertie augmente le rayon de convergence : plus

$\mu$ est petit plus le rayon de convergence est grand.

Pour trouver les autres solutions il faut choisir un autre point initial.

La méthode de Newton-Raphson

$${\mathbf{x}}^{k+1}=\mathbf{x}-\frac{f({\mathbf{x}}_{\mathbf{n}})}{f^{\prime }({\mathbf{x}}_{\mathbf{n}})}$$

La méthode de Newton est une méthode de point fixe.

$${\mathbf{x}}^{k+1}=g({\mathbf{x}}^k)$$
où

$g(\mathbf{x}\mathbf{)}=\mathbf{x}-\frac{f(\mathbf{x})}{f^{\prime }(\mathbf{x})}$

On souhaite résoudre notre systeme non linéaire. Grâce à la formule ci-dessus, on a en 1D:

$$f^{\prime }(x^k)(x^{k+1}-x^k)=-f(x^k)$$
Ce qui donne en

$n$ dimensions:

$$J_{\mathbf{f}}({\mathbf{x}}^k)({\mathbf{x}}^{k+1}-{\mathbf{x}}^k)=-\mathbf{f}({\mathbf{x}}^k)$$
avec

$J_{\mathbf{f}}$ la matrice Jacobienne de

$\mathbf{f}$ :

$$J_{\mathbf{f}}\left(\mathbf{x}\right)=\left(\begin{array}{ccc}{\displaystyle \frac{\partial f_1}{\partial x_1}}& \cdots & {\displaystyle \frac{\partial f_1}{\partial x_n}}\\ ⋮& \ddots & ⋮\\ {\displaystyle \frac{\partial f_n}{\partial x_1}}& \cdots & {\displaystyle \frac{\partial f_n}{\partial x_n}}\end{array}\right)$$
Notre système non linéaire avec

$f$ une fonction vectorielle:

$$\mathbf{f}(\mathbf{x})=A(\mathbf{x})\mathbf{x}-\mathbf{b}$$

Exemple

On reprend le système matriciel précèdent. La matrice Jacobienne de la fonciton

$f$ est :

$$J_{\mathbf{f}}(\mathbf{x})=\left[\begin{array}{cc}2x_0-2x_1& 2x_1-2x_0\\ 2x_0+2x_1& 2x_0+2x_1\end{array}\right]$$

def f(x):
    return A(x) @ x - b

def J_f(x):
    return 2 * np.array([[x[0] - x[1], x[1] - x[0]],
                         [x[0] + x[1], x[0] + x[1]]])

x = np.array([3, 2])
for i in range(30):
    delta = lin.solve(J_f(x), -f(x))
    x = x + delta

...
x^29 =  [2.05693134 1.05693134]

On converge (moins vite) où la methode du point fixe oscille sans converger.

• Le coût de la construction de la matrice Jacobienne peut être tres élevé.
• Pour aller plus vite on peut recalculer la matrice toutes les 3 iterations ou plus.
• Il s'agit d'une matrice pleine qui rend compliqué la resolution du systeme (une methode de gradient ne marchera pas)