Interpolation de surfaces

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

def newton_gregory_1d(x_data, y_data, x_eval):
  """
  Interpolation 1D par Newton-Gregory.
  """
  n = len(x_data)
  h = x_data[1] - x_data[0]

  # Table de différences
  diff = [y_data.copy()]
  for k in range(1, n):
      new_diff = []
      for i in range(n - k):
          new_diff.append(diff[-1][i + 1] - diff[-1][i])
      diff.append(np.array(new_diff))

  # Évaluation
  s = (x_eval - x_data[0]) / h
  result = diff[0][0]
  prod = 1.0

  for k in range(1, n):
      prod *= (s - k + 1) / k
      result += prod * diff[k][0]

  return result

def interpolation_2d_etapes(x_data, y_data, z_grid, x_eval, y_eval):
  """
  Interpolation 2D par la méthode des étapes successives.
  """
  n = len(x_data)

  # Étape 1 : Interpoler en y pour chaque x
  z_at_y = np.zeros(n)
  for i in range(n):
      z_at_y[i] = newton_gregory_1d(y_data, z_grid[i, :], y_eval)

  # Étape 2 : Interpoler en x
  result = newton_gregory_1d(x_data, z_at_y, x_eval)

  return result

def interpolation_2d_lagrange(x_data, y_data, z_grid, x_eval, y_eval):
  """
  Interpolation 2D par polynôme de Lagrange bidimensionnel.
  """
  n, m = len(x_data), len(y_data)
  result = 0.0

  for i in range(n):
      Lx = 1.0
      for k in range(n):
          if k != i:
              Lx *= (x_eval - x_data[k]) / (x_data[i] - x_data[k])

      for j in range(m):
          Ly = 1.0
          for l in range(m):
              if l != j:
                  Ly *= (y_eval - y_data[l]) / (y_data[j] - y_data[l])

          result += z_grid[i, j] * Lx * Ly

  return result

# Données de l'exemple
x_data = np.array([1.0, 1.5, 2.0])
y_data = np.array([0.2, 0.3, 0.4, 0.5])

z_grid = np.array([
  [0.640, 1.003, 1.359, 1.703],
  [0.990, 1.524, 2.045, 2.549],
  [1.568, 2.384, 3.177, 3.943]
])

x_eval, y_eval = 1.6, 0.33

# Calculs
result_etapes = interpolation_2d_etapes(x_data, y_data, z_grid, x_eval, y_eval)
result_lagrange = interpolation_2d_lagrange(x_data, y_data, z_grid, x_eval, y_eval)
exact = np.exp(x_eval) * np.sin(y_eval) + y_eval - 0.1

print(f"Méthode par étapes : f({x_eval}, {y_eval}) = {result_etapes:.4f}")
print(f"Méthode Lagrange 2D : f({x_eval}, {y_eval}) = {result_lagrange:.4f}")
print(f"Valeur exacte : f({x_eval}, {y_eval}) = {exact:.4f}")

# ============================================
# VISUALISATION 3D
# ============================================

# Grille fine pour la surface interpolée
x_fine = np.linspace(x_data[0], x_data[-1], 50)
y_fine = np.linspace(y_data[0], y_data[-1], 50)
X_fine, Y_fine = np.meshgrid(x_fine, y_fine)

# Calculer la surface interpolée
Z_interp = np.zeros_like(X_fine)
for i in range(len(x_fine)):
  for j in range(len(y_fine)):
      Z_interp[j, i] = interpolation_2d_lagrange(
          x_data, y_data, z_grid, x_fine[i], y_fine[j]
      )

# Surface exacte pour comparaison
Z_exact = np.exp(X_fine) * np.sin(Y_fine) + Y_fine - 0.1

# Grille des points de données
X_grid, Y_grid = np.meshgrid(x_data, y_data)

# Création de la figure 3D
fig = plt.figure(figsize=(14, 5))

# Subplot 1 : Surface interpolée avec points de données
ax1 = fig.add_subplot(131, projection='3d')
ax1.plot_surface(X_fine, Y_fine, Z_interp, cmap='viridis', alpha=0.7)
ax1.scatter(X_grid.flatten(), Y_grid.flatten(), z_grid.T.flatten(),
         color='red', s=50, label='Points de données')
ax1.scatter([x_eval], [y_eval], [result_lagrange],
         color='yellow', s=100, marker='*', label='Point évalué')
ax1.set_xlabel('x')
ax1.set_ylabel('y')
ax1.set_zlabel('z')
ax1.set_title('Surface interpolée')
ax1.legend()

# Subplot 2 : Surface exacte
ax2 = fig.add_subplot(132, projection='3d')
ax2.plot_surface(X_fine, Y_fine, Z_exact, cmap='plasma', alpha=0.7)
ax2.scatter(X_grid.flatten(), Y_grid.flatten(), z_grid.T.flatten(),
         color='red', s=50)
ax2.set_xlabel('x')
ax2.set_ylabel('y')
ax2.set_zlabel('z')
ax2.set_title('Surface exacte f(x,y) = eˣsin(y) + y - 0.1')

# Subplot 3 : Erreur d'interpolation
ax3 = fig.add_subplot(133, projection='3d')
Z_error = np.abs(Z_interp - Z_exact)
surf = ax3.plot_surface(X_fine, Y_fine, Z_error, cmap='hot', alpha=0.8)
ax3.set_xlabel('x')
ax3.set_ylabel('y')
ax3.set_zlabel('Erreur')
ax3.set_title(f'Erreur (max = {Z_error.max():.4f})')
fig.colorbar(surf, ax=ax3, shrink=0.5)

plt.tight_layout()
plt.savefig('interpolation_2d_3d.png', dpi=150)
plt.show()