Série 5 : Méthodes itératives pour la résolution de systèmes linéaires

import numpy as np
import matplotlib.pyplot as plt

n = 5  # ou tout autre nombre approprié
A = 2 * np.diag(np.ones(n)) - np.diag(np.ones(n - 1), -1) - np.diag(np.ones(n - 1), 1)
np.min(np.linalg.eigvals(A))

D = np.diag(np.diag(A))
E = -np.tril(A, -1)
F = -np.triu(A, 1)
NMAX = 10000
TOL = 1e-9
X0 = np.zeros(n)
b = A @ np.ones(n)

# Note: La fonction sor() doit être définie ou importée
count = 0
omega_range = np.arange(0, 2.1, 0.1)
rho = []
iter = []
for omega in omega_range:
    count += 1
    B = np.linalg.inv(np.eye(n) - omega * np.linalg.inv(D) @ E) @ ((1 - omega) * np.eye(n) + omega * np.linalg.inv(D) @ F)
    rho.append(max(abs(np.linalg.eigvals(B))))
    # X, ITER = sor(A, b, X0, NMAX, TOL, omega)
    # iter.append(ITER)

plt.figure(1)
plt.plot(omega_range, rho)
plt.figure(2)
plt.plot(omega_range, iter)
plt.show()

m, argm = min((val, idx) for (idx, val) in enumerate(rho))
omega_opt = omega_range[argm]
iter_at_omega_opt = iter[argm]

arg1 = np.where(omega_range == 1)[0][0]
iter_at_omega_1 = iter[arg1]

import numpy as np

# Définir A et b
A = np.array([[2, -1], [-1, 2]])
b = np.array([1, 1])

# Vérifier si A est symétrique définie positive
if np.allclose(A, A.T) and np.all(np.linalg.eigvals(A) > 0):
    print("A est symétrique définie positive.")
else:
    print("A n'est pas symétrique définie positive.")

# Méthode du gradient conjugué
def conjugate_gradient(A, b, x0, tol=1e-10, max_iter=1000):
    x = x0
    r = b - np.dot(A, x)
    d = r.copy()
    rs_old = np.dot(r.T, r)
    x_list = [x.copy()]

    for i in range(max_iter):
        Ad = np.dot(A, d)
        alpha = rs_old / np.dot(d.T, Ad)
        x += alpha * d
        r -= alpha * Ad
        rs_new = np.dot(r.T, r)
        if np.sqrt(rs_new) < tol:
            break
        d = r + (rs_new / rs_old) * d
        rs_old = rs_new
        x_list.append(x.copy())

    return x, x_list

# Appliquer la méthode
x0 = np.array([1.5, 2])
x_sol, x_list = conjugate_gradient(A, b, x0, max_iter=2)
print( f"x_list={x_list}, x_sol={x_sol}" )

# Fonctionnelle d'énergie Phi
def phi(x):
    return 0.5 * np.dot(x, A @ x) - np.dot(b, x)

# Tracer les lignes de niveau de Phi
x_vals = np.linspace(0, 2, 400)
y_vals = np.linspace(0, 2, 400)
X, Y = np.meshgrid(x_vals, y_vals)
Z = np.array([[phi(np.array([x, y])) for x, y in zip(x_row, y_row)] for x_row, y_row in zip(X, Y)])

phi_levels = sorted([phi(x) for x in x_list])

import plotly.graph_objs as go


# Créer la figure avec les lignes de niveau
fig = go.Figure(data =
    go.Contour(
        z=Z,
        x=x_vals,
        y=y_vals,
        contours=dict(
            start=min(phi_levels),
            end=max(phi_levels),
            size=int(max(phi_levels)/len(phi_levels))
        ),
        colorscale='Blues'
    )
)

# Ajouter les points et les annotations
for i, point in enumerate(x_list):
    fig.add_trace(go.Scatter(x=[point[0]], y=[point[1]], mode='markers+text',
                             text=[f'x^{i}'],
                             textposition="bottom center",
                             marker=dict(color='red', size=10)))

# Ajouter x_sol
fig.add_trace(go.Scatter(x=[x_sol[0]], y=[x_sol[1]], mode='markers',
                         marker=dict(color='green', size=10)))

# Tracer les segments entre les points
for i in range(len(x_list) - 1):
    fig.add_trace(go.Scatter(x=[x_list[i][0], x_list[i+1][0]],
                             y=[x_list[i][1], x_list[i+1][1]],
                             mode='lines',
                             line=dict(color='red')))
fig.add_trace(go.Scatter(x=[x_list[-1][0], x_sol[0]],
                         y=[x_list[-1][1], x_sol[1]],
                         mode='lines',
                         line=dict(color='red')))

# Mise en forme finale de la figure
fig.update_layout(title='Lignes de niveau de Phi et points de la méthode du gradient conjugué',
                  xaxis_title='x1',
                  yaxis_title='x2')
fig.show()