Splines cubiques — Conditions frontières et résolution

import numpy as np

def resoudre_tridiagonal(a, b, c, d):
  """
  Résout un système tridiagonal Ax = d par l'algorithme de Thomas.

  Paramètres:
      a : sous-diagonale (a[0] non utilisé)
      b : diagonale principale
      c : sur-diagonale (c[n-1] non utilisé)
      d : second membre

  Retourne:
      x : solution
  """
  n = len(d)
  c_prime = np.zeros(n)
  d_prime = np.zeros(n)

  # Forward sweep
  c_prime[0] = c[0] / b[0]
  d_prime[0] = d[0] / b[0]

  for i in range(1, n):
      denom = b[i] - a[i] * c_prime[i - 1]
      c_prime[i] = c[i] / denom if i < n - 1 else 0
      d_prime[i] = (d[i] - a[i] * d_prime[i - 1]) / denom

  # Back substitution
  x = np.zeros(n)
  x[n - 1] = d_prime[n - 1]

  for i in range(n - 2, -1, -1):
      x[i] = d_prime[i] - c_prime[i] * x[i + 1]

  return x

def spline_naturelle(x, y):
  """
  Calcule une spline cubique naturelle.

  Paramètres:
      x : abscisses (n points)
      y : ordonnées

  Retourne:
      S : dérivées secondes aux noeuds
      a, b, c, d : coefficients des polynômes
  """
  n = len(x)
  h = np.diff(x)

  # Construire le système pour S_2, ..., S_{n-1}
  # (S_1 = S_n = 0 pour spline naturelle)
  m = n - 2  # nombre d'inconnues intérieures

  if m == 0:
      # Cas de 2 points : segment linéaire
      S = np.zeros(n)
      a = np.zeros(n - 1)
      b = np.zeros(n - 1)
      c = (y[1:] - y[:-1]) / h
      d = y[:-1].copy()
      return S, a, b, c, d

  # Diagonales du système
  diag_a = np.zeros(m)
  diag_b = np.zeros(m)
  diag_c = np.zeros(m)
  rhs = np.zeros(m)

  for i in range(m):
      idx = i + 1  # indice dans le tableau original
      diag_b[i] = 2 * (h[idx - 1] + h[idx])
      if i > 0:
          diag_a[i] = h[idx - 1]
      if i < m - 1:
          diag_c[i] = h[idx]
      rhs[i] = 6 * ((y[idx + 1] - y[idx]) / h[idx] - (y[idx] - y[idx - 1]) / h[idx - 1])

  # Résoudre
  S_int = resoudre_tridiagonal(diag_a, diag_b, diag_c, rhs)

  # Reconstruire S complet
  S = np.zeros(n)
  S[1:n-1] = S_int

  # Calculer les coefficients
  a = np.zeros(n - 1)
  b = np.zeros(n - 1)
  c = np.zeros(n - 1)
  d = np.zeros(n - 1)

  for i in range(n - 1):
      a[i] = (S[i + 1] - S[i]) / (6 * h[i])
      b[i] = S[i] / 2
      c[i] = (y[i + 1] - y[i]) / h[i] - h[i] * (S[i + 1] + 2 * S[i]) / 6
      d[i] = y[i]

  return S, a, b, c, d

def evaluer_spline(x_data, a, b, c, d, x_eval):
  """Évalue la spline en un ou plusieurs points."""
  x_eval = np.atleast_1d(x_eval)
  result = np.zeros_like(x_eval)

  for j, x in enumerate(x_eval):
      # Trouver l'intervalle
      i = np.searchsorted(x_data, x) - 1
      i = max(0, min(i, len(a) - 1))

      dx = x - x_data[i]
      result[j] = a[i] * dx**3 + b[i] * dx**2 + c[i] * dx + d[i]

  return result if len(result) > 1 else result[0]

# Exemple
x = np.array([1, 2, 3, 4], dtype=float)
y = np.array([4, -2, 3, 1], dtype=float)

S, a, b, c, d = spline_naturelle(x, y)

print("Dérivées secondes S :")
print(f"  S = {S}")

print("\nCoefficients des polynômes :")
for i in range(len(a)):
  print(f"  Q_{i+1}(x) : a={a[i]:.2f}, b={b[i]:.2f}, c={c[i]:.2f}, d={d[i]:.2f}")

# Évaluation
x_test = np.linspace(1, 4, 50)
y_test = evaluer_spline(x, a, b, c, d, x_test)

print(f"\nSpline en x=1.5 : {evaluer_spline(x, a, b, c, d, 1.5):.4f}")
print(f"Spline en x=2.5 : {evaluer_spline(x, a, b, c, d, 2.5):.4f}")