2026 Contrôle Continu – M64 – Data Science & Scientific Computing 6 : Schémas numériques pour les problèmes de Cauchy

2026 Contrôle Continu

Auteur·rice

Gloria Faccanoni

Date de publication

24 avril 2026

⏱ 2h, déposer le notebook dans Moodle (si 1/3 temps additionnel, 2h40)

📚 Autorisés : notebooks de cours, notes personnelles.

Compléter ici :

  • NOM :
  • Prénom :
  • N° étudiant :
Code 
import numpy as np
import matplotlib.pyplot as plt

from scipy.optimize import fsolve
from scipy.integrate import solve_ivp, odeint

import sympy as sp

from IPython.display import display, Markdown, Math

EXERCICE 1 : Questionnaire

Justifier chaque réponse en une ligne ou deux.

Question 1

Une méthode est d’ordre \(\omega\). Si on divise le pas de discrétisation \(h\) par 2, comment évolue l’erreur globale ?

Réponse 1

L’erreur globale est divisée par \(2^{\omega}\). En effet, si \(e(h)\approx C h^{\omega}\) alors \(e(h/2)\approx C (h/2)^{\omega} = C h^{\omega} / 2^{\omega} = e(h) / 2^{\omega}\).

Une autre manière de le voir est de se souvenir que, en échelle log-log, la pente de la courbe de l’erreur est égale à l’ordre \(\omega\). Donc, si on divise \(h\) par 2, on se déplace de \(\ln(2)\) vers la gauche, et l’erreur diminue de \(\omega \ln(2) = \ln(2^{\omega})\), soit une division par \(2^{\omega}\).

Question 2

Sur l’équation test \(y'(t)=-\beta y(t)\), \(\beta>0\), quelles méthodes peuvent produire des oscillations numériques pour des pas de discrétisation \(h\) trop grands parmi Euler Explicite (EE), Euler Implicite (EI), Crank-Nicolson (CN), Heun (H), Euler Modifié (EM) ?

Réponse 2

Les méthodes d’Euler explicite (EE), d’Euler Modifié (EM) et de Heun (H) peuvent produire des oscillations numériques pour des pas trop grands car elles ne sont A-stable que pour des pas suffisamment petits. Les méthodes d’Euler implicite (EI) et de Crank-Nicolson (CN) sont inconditionnellement A-stables et ne produisent pas d’oscillations numériques quel que soit le pas.

Question 3 - sujet A

Une méthode multi-pas linéaire est consistante mais non convergente. Quelle propriété manque nécessairement pour avoir la convergence ?

Réponse 3 - sujet A

La zéro-stabilité.

Question 3 - sujet B

Une méthode multi-pas linéaire est zéro-stable mais non convergente. Quelle propriété manque nécessairement pour avoir la convergence ?

Réponse 3 - sujet B

La consistance.

Question 4

La méthode multi-pas linéaire \(u_{n+1} = - u_{n} + u_{n-1} + u_{n-2} + h \varphi_{n}\) est-elle convergente ?

Réponse 4

Une méthode MPL est convergente ssi elle est consistante et zéro-stable. Or, la méthode donnée n’est ni l’une ni l’autre. En effet :

  • elle n’est pas zéro-stable, car le polynôme caractéristique associé est \(\rho(r) = r^3 + r^2 - r - 1 = (r-1)(r+1)^2\), dont une racine (\(-1\)) est sur le cercle unité avec multiplicité \(2\) ;
  • elle n’est pas consistante non plus, car \(\xi(0)=0\) mais \(\xi(1) = -a_1 - 2a_2 + b_0 + b_1 + b_2 + b_{-1} = -1 -2 + 1 = -2 \neq 1\).

Question 5

Quelle est la principale différence entre une méthode de Runge–Kutta et une méthode multi-pas linéaire ?

Réponse 5

Les méthodes de Runge-Kutta sont des méthodes à un pas et reconstruisent l’information dans des points qui ne sont pas dans la grille via les étages internes.

Les méthodes multi-pas n’utilisent que des points de la grille mais nécessitent plusieurs valeurs précédentes pour effectuer le calcul.

EXERCICE 2 - sujet A : Analyse d’un schéma RK explicite à 2 étages

On vous donne cette matrice de Butcher pour une méthode de Runge-Kutta explicite à deux étages : \( \def\arraystretch{2} \begin{array}{c|cc} 0 & 0 & 0 \\ 1/2 & 1/3 & 0 \\ \hline & 1/2 & 1/2 \end{array} \)

  • Quelle condition de consistance n’est pas satisfaite ? Corrigez la matrice \(\mathbb{A}\) pour que la méthode soit consistante (la méthode doit rester explicite).
  • La méthode corrigée est-elle d’ordre 2 ? Si non, corrigez les poids \(\mathbf{b}\) pour que la méthode soit d’ordre 2.
  • La méthode finale est-elle inconditionnellement A-stable ?

  • La condition de consistance qui n’est pas satisfaite est \(\sum_{j=1}^{s} a_{ij} = c_i\) pour \(i=2\). En effet, \(a_{21} + a_{22} = 1/3 + 0 = 1/3 \neq c_2 = 1/2\). Pour corriger la matrice \(\mathbb{A}\), on doit remplacer \(a_{21}\) par \(1/2\) : \( \def\arraystretch{2} \begin{array}{c|cc} 0 & 0 & 0 \\ 1/2 & \mathbf{1/2} & 0 \\ \hline & 1/2 & 1/2 \end{array} \)

  • La méthode corrigée n’est pas d’ordre 2 car la condition \(\sum_{i=1}^{s} b_i c_i = 1/2\) n’est pas satisfaite : \(b_1 c_1 + b_2 c_2 = 1/2 \cdot 0 + 1/2 \cdot 1/2 = 1/4 \neq 1/2\). Pour corriger les poids \(\mathbf{b}\), on doit remplacer \(b_2\) par \(1\) et donc \(b_1\) par \(0\): \( \def\arraystretch{2} \begin{array}{c|cc} 0 & 0 & 0 \\ 1/2 & 1/2 & 0 \\ \hline & \mathbf{0} & \mathbf{1} \end{array} \)

  • La méthode finale n’est pas inconditionnellement A-stable car aucune méthode explicite ne peut l’être. Elle est seulement stable pour des pas \(h\) suffisamment petits.

    NON DEMANDÉ : En effet, si on applique la méthode à l’équation test \(y' = -\beta y\), on obtient : \( \begin{cases} K_0 = \varphi(t_n, u_n) = -\beta u_n\\ K_1 = \varphi\left(t_n + \frac{h}{2}, u_n + \frac{h}{2} K_0\right) = -\beta \left( u_n + \frac{h}{2} K_0 \right) = K_0 - \frac{h \beta}{2} K_0 \\ u_{n+1} = u_n + h K_1 = u_n + h K_0 - \frac{h^2 \beta}{2} K_0 = \left( 1 - h \beta + \frac{h^2 \beta^2}{2} \right) u_n \end{cases} \) La fonction d’amplification est donc \(q(x) = 1 - x + \frac{x^2}{2}\) avec \(x = h \beta\). Le critère de stabilité impose \(|q(x)| \leq 1\), ce qui n’est vérifié que pour \(h\beta<2\). Pour ces valeurs de \(h\beta\), la convergence est monotone.

    Vérifions ces calculs avec sympy :

Code 
def calculer_q_rk_general_v2(s, A=None, b=None):
    if A is None:  A = sp.Matrix(s, s, lambda i, j: sp.symbols(f'a_{i+1}{j+1}'))
    if b is None:  b = sp.Matrix(s, 1, lambda i, j: sp.symbols(f'b_{i+1}'))
    ones = sp.Matrix.ones(s, 1)
    I = sp.eye(s)
    q = sp.det(I + x * A - x * ones * b.T) / sp.det(I + x * A)
    return sp.simplify(q)


A = sp.Matrix([[0, 0], [sp.Rational(1, 2), 0]])
b = sp.Matrix([0, 1])
display(Markdown(f"\(A = {sp.latex(A)}, \\quad b^T = {sp.latex(b.T)}\)"))

x = sp.symbols('x')
q = calculer_q_rk_general_v2(s = len(b), A=A, b=b)
display(Markdown(f"\(q(x) = {sp.latex(q)}\)"))


q_np = sp.lambdify(x, q, 'numpy')

x_max = 2.5
xx = np.linspace(0, x_max, 500)
plt.figure(figsize=(10, 4))
plt.plot(xx, q_np(xx), label='q(x)', color='blue')

# Tracé des limites de stabilité (|q(x)| = 1)
plt.axhline(y=1, color='black', linestyle='-', alpha=0.5)
plt.axhline(y=-1, color='black', linestyle='-', alpha=0.5)
plt.fill_between(xx, -1, 1, color='gray', alpha=0.1, label='Zone de A-Stabilité')

# Configuration des axes
plt.xlabel(r'$x = \beta h$')
plt.ylabel(r'$q(x)$')
plt.ylim(-2, 2) # On zoome sur la zone intéressante [-1, 1]
plt.xlim(0, x_max)
plt.grid(True, which='both', linestyle='--', linewidth=0.5)
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left');


plt.tight_layout()
plt.show()

\(A = \left[\begin{matrix}0 & 0\\\frac{1}{2} & 0\end{matrix}\right], \quad b^T = \left[\begin{matrix}0 & 1\end{matrix}\right]\)

\(q(x) = \frac{x^{2}}{2} - x + 1\)

EXERCICE 2 - sujet B : Analyse d’un schéma RK explicite à 2 étages

On vous donne cette matrice de Butcher pour une méthode de Runge-Kutta explicite à deux étages : \( \def\arraystretch{2} \begin{array}{c|cc} 0 & 0 & 0 \\ 1/3 & 1/2 & 0 \\ \hline & 1/2 & 1/2 \end{array} \)

  • Quelle condition de consistance n’est pas satisfaite ? Corrigez la matrice \(\mathbb{A}\) pour que la méthode soit consistante (la méthode doit rester explicite).
  • La méthode corrigée est-elle d’ordre 2 ? Si non, corrigez les poids \(\mathbf{b}\) pour que la méthode soit d’ordre 2.
  • La méthode finale est-elle inconditionnellement A-stable ?

  • La condition de consistance qui n’est pas satisfaite est \(\sum_{j=1}^{s} a_{ij} = c_i\) pour \(i=2\). En effet, \(a_{21} + a_{22} = 1/2 + 0 = 1/2 \neq c_2 = 1/3\). Pour corriger la matrice \(\mathbb{A}\), on doit remplacer \(a_{21}\) par \(1/3\) : \( \def\arraystretch{2} \begin{array}{c|cc} 0 & 0 & 0 \\ 1/3 & \mathbf{1/3} & 0 \\ \hline & 1/2 & 1/2 \end{array} \)

  • La méthode corrigée n’est pas d’ordre 2 car la condition \(\sum_{i=1}^{s} b_i c_i = 1/2\) n’est pas satisfaite : \(b_1 c_1 + b_2 c_2 = 1/2 \cdot 0 + 1/3 \cdot 1/2 = 1/6 \neq 1/2\). Pour corriger les poids \(\mathbf{b}\), on doit remplacer \(b_2\) par \(3/2\) et donc \(b_1\) par \(-1/2\): \( \def\arraystretch{2} \begin{array}{c|cc} 0 & 0 & 0 \\ 1/3 & 1/3 & 0 \\ \hline & \mathbf{-1/2} & \mathbf{3/2} \end{array} \)

  • La méthode finale n’est pas inconditionnellement A-stable car aucune méthode explicite ne peut l’être. Elle est seulement stable pour des pas \(h\) suffisamment petits.

    NON DEMANDÉ : En effet, si on applique la méthode à l’équation test \(y' = -\beta y\), on obtient : \( \begin{cases} K_0 = \varphi(t_n, u_n) = -\beta u_n\\ K_1 = \varphi\left(t_n + \frac{h}{3}, u_n + \frac{h}{3} K_0\right) = -\beta \left( u_n + \frac{h}{3} K_0 \right) = K_0 - \frac{h \beta}{3} K_0 \\ u_{n+1} = u_n - \frac{h}{2} K_0 + \frac{3h}{2} K_1 = u_n - \frac{h}{2} K_0 + \frac{3h}{2} \left( K_0 - \frac{h \beta}{3} K_0 \right) = u_n + h K_0 - \frac{h^2 \beta}{2} K_0 = \left( 1 - h \beta + \frac{h^2 \beta^2}{2} \right) u_n \end{cases} \) La fonction d’amplification est donc \(q(x) = 1 - x + \frac{x^2}{3}\) avec \(x = h \beta\). Le critère de stabilité impose \(|q(x)| \leq 1\), ce qui n’est vérifié que pour \(h\beta<3\). Pour ces valeurs de \(h\beta\), la convergence est monotone.

    Vérifions ces calculs avec sympy :

Code 
def calculer_q_rk_general_v2(s, A=None, b=None):
    if A is None:  A = sp.Matrix(s, s, lambda i, j: sp.symbols(f'a_{i+1}{j+1}'))
    if b is None:  b = sp.Matrix(s, 1, lambda i, j: sp.symbols(f'b_{i+1}'))
    ones = sp.Matrix.ones(s, 1)
    I = sp.eye(s)
    q = sp.det(I + x * A - x * ones * b.T) / sp.det(I + x * A)
    return sp.simplify(q)


A = sp.Matrix([[0, 0], [sp.Rational(1, 3), 0]])
b = sp.Matrix([sp.Rational(-1, 2), sp.Rational(3, 2)])
display(Markdown(f"\(A = {sp.latex(A)}, \\quad b^T = {sp.latex(b.T)}\)"))

x = sp.symbols('x')
q = calculer_q_rk_general_v2(s = len(b), A=A, b=b)
display(Markdown(f"\(q(x) = {sp.latex(q)}\)"))

q_np = sp.lambdify(x, q, 'numpy')

x_max = 2.5
xx = np.linspace(0, x_max, 500)
plt.figure(figsize=(10, 4))
plt.plot(xx, q_np(xx), label='q(x)', color='blue')

# Tracé des limites de stabilité (|q(x)| = 1)
plt.axhline(y=1, color='black', linestyle='-', alpha=0.5)
plt.axhline(y=-1, color='black', linestyle='-', alpha=0.5)
plt.fill_between(xx, -1, 1, color='gray', alpha=0.1, label='Zone de A-Stabilité')

# Configuration des axes
plt.xlabel(r'$x = \beta h$')
plt.ylabel(r'$q(x)$')
plt.ylim(-2, 2) # On zoome sur la zone intéressante [-1, 1]
plt.xlim(0, x_max)
plt.grid(True, which='both', linestyle='--', linewidth=0.5)
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left');


plt.tight_layout()
plt.show()

\(A = \left[\begin{matrix}0 & 0\\\frac{1}{3} & 0\end{matrix}\right], \quad b^T = \left[\begin{matrix}- \frac{1}{2} & \frac{3}{2}\end{matrix}\right]\)

\(q(x) = \frac{x^{2}}{2} - x + 1\)

EXERCICE 3 - sujets A, B, C, D : Système de poursuite

Le système suivant décrit l’évolution de la position d’un chien \((x(t), y(t))\) en fonction de celle de son maître \((x_m(t), y_m(t))\) :

\( \begin{cases} x'(t) = kv \dfrac{X(t)}{r(t)}, \\ y'(t) = kv \dfrac{Y(t)}{r(t)}, \end{cases} \qquad \text{avec } \qquad \begin{cases} X(t) = x_m(t) - x(t), \\ Y(t) = y_m(t) - y(t), \\ r(t) = \sqrt{X^2(t) + Y^2(t)}, \end{cases} \) où \(k\) et \(v\) sont des constantes positives.

Le système est fermé par le point de départ du chien \((x(0),y(0))\) et la trajectoire du maître \((x_m(t),y_m(t))\).

L’exercice consiste à analyser ce système dans l’une des quatre configurations suivantes (où \(a\) est une constante positive).
Identifiez votre sujet parmi les quatre configurations suivantes : si votre numéro étudiant est de la forme \(4n\) : sujet A ; \(4n+1\) : sujet B ; \(4n+2\) : sujet C ; \(4n+3\) : sujet D.

Sujet Position initiale du chien \((x(0),y(0))\) Trajectoire du maître \((x_m(t),y_m(t))\) \(E(t)\)
A \((0, 0)\) \((a, vt)\) \(k \ln(r(t) + Y(t)) - (k-1) \ln(X(t))\)
B \((a, 0)\) \((0, vt)\) \(k \ln(r(t) + Y(t)) - (k-1) \ln(X(t))\)
C \((0, 0)\) \((vt, a)\) \(k \ln(r(t) + X(t)) - (k-1) \ln(Y(t))\)
D \((0, a)\) \((vt, 0)\) \(k \ln(r(t) + X(t)) - (k-1) \ln(Y(t))\)
Code 
num_etudiant = 123457  # Remplacez par votre numéro étudiant

# --- CALCUL DU SUJET ---
sujet = "ABCD"[num_etudiant % 4]  # Choix automatique du sujet en fonction du numéro étudiant
display(Markdown(f"**Sujet <mark> {sujet} </mark> sélectionné pour le numéro étudiant {num_etudiant}**"))
# ----------------------

Sujet B sélectionné pour le numéro étudiant 123457

Travail à effectuer

Partie 1 : Propriétés analytiques

  1. Norme de la vitesse : démontrer que \(\big(x'(t)\big)^2 + \big(y'(t)\big)^2 = \text{constante}\).
  2. Invariant : à l’aide de sympy, vérifier que la quantité \(E(t)\) associée à votre sujet est un invariant du système.

Partie 2 : Simulations numériques

Utiliser les paramètres suivants :

  • Pas de temps \(h = 0.01\)
  • Vitesse du maître \(v = 1\)
  • Vitesse du chien : \(kv\) avec \(k=1.5\)
  • Distance caractéristique : \(a=5\)

Condition d’arrêt : stopper la simulation quand \(r(t) < \varepsilon = 10^{-2}\).

  1. Méthode d’Euler explicite :
    • Implémenter le schéma.
    • Graphique 1 : tracer les trajectoires \((x, y)\) du chien et \((x_m,y_m)\) du maître
    • Graphique 2 : tracer \(t\mapsto E(t)\) pour \(t \in [0,5]\). Analyser la qualité de la conservation.
  2. Méthode de Heun (si num. étudiant impair) OU Euler Modifiée (si num. étudiant pair) :
    • Répéter les mêmes études.

Norme de la vitesse

En utilisant les définitions de \(x'(t)\) et \(y'(t)\), on obtient :

\( (x'(t))^2 + (y'(t))^2 = \Bigg(kv \frac{X(t)}{r(t)}\Bigg)^2 + \Bigg(kv \frac{Y(t)}{r(t)}\Bigg)^2 = (kv)^2 \, \frac{X^2(t) + Y^2(t)}{r^2(t)} = (kv)^2 \quad \text{pour tout } t \geq 0. \)

Ainsi, la norme de la vitesse du chien est constante et égale à \((kv)^2\).

Energie :

\(E(t)\) est un invariant exact du système : sa dérivée est nulle pour tout \(t\).

Je conseille d’utiliser sympy pour vérifier que \(E'(t) = 0\) en effectuant les calculs symboliques. Mais il est bien sûr possible de faire les calculs à la main.

Pour les sujets A et B, on a \( \begin{aligned} X'(t)&=-x'(t)=-\frac{kv}{r(t)}X(t), \\ Y'(t)&=v-y'(t)=v-\frac{kv}{r(t)}Y(t), \\ r'(t) &= \frac{X(t)X'(t) + Y(t)Y'(t)}{r(t)} = v \frac{Y(t)}{r(t)} - kv. \end{aligned} \) Ainsi \( \begin{aligned} E'(t) &= k \frac{r'(t) + Y'(t)}{r(t)+Y(t)} - (k-1) \frac{X'(t)}{X(t)} = k \frac{v \frac{Y(t)}{r(t)} - kv + v - \frac{kv}{r(t)}Y(t)}{r(t)+Y(t)} + (k-1) \frac{kv}{r(t)} \\ &= k \frac{(1-k)v \frac{Y(t)}{r(t)} + (1-k)v}{r(t)+Y(t)} + (k-1) \frac{kv}{r(t)} = k v (1-k) \frac{ \frac{Y(t)}{r(t)} + 1}{r(t)+Y(t)} + (k-1) \frac{kv}{r(t)} \\ &= \frac{kv}{r(t)} (1-k) + (k-1) \frac{kv}{r(t)} =0. \end{aligned} \)

Pour les sujets C et D, on a \( \begin{aligned} X'(t) &= v - x'(t) = v - \frac{kv}{r(t)}X(t), \\ Y'(t) &= -y'(t) = -\frac{kv}{r(t)}Y(t), \\ r'(t) &= \frac{X(t)X'(t) + Y(t)Y'(t)}{r(t)} = v \frac{X(t)}{r(t)} - kv. \end{aligned} \) Ainsi \( \begin{aligned} E'(t) &= k \frac{r'(t) + X'(t)}{r(t)+X(t)} - (k-1) \frac{Y'(t)}{Y(t)} = k \frac{v \frac{X(t)}{r(t)} - kv + v - \frac{kv}{r(t)}X(t)}{r(t)+X(t)} + (k-1) \frac{kv}{r(t)} \\ &= k \frac{(1-k)v \frac{X(t)}{r(t)} + (1-k)v}{r(t)+X(t)} + (k-1) \frac{kv}{r(t)} = k v (1-k) \frac{ \frac{X(t)}{r(t)} + 1}{r(t)+X(t)} + (k-1) \frac{kv}{r(t)} \\ &= \frac{kv}{r(t)} (1-k) + (k-1) \frac{kv}{r(t)} =0. \end{aligned} \)

Code 
t, v_sp, k_sp, a_sp = sp.symbols('t v k a', positive=True, real=True)
x = sp.Function('x')(t)
y = sp.Function('y')(t)

# --- Sélection de la configuration du maître en fonction du sujet ---
xm_sp = v_sp * t if sujet in ['C', 'D'] else (a_sp if sujet=='A' else 0)
ym_sp = v_sp * t if sujet in ['A', 'B'] else (a_sp if sujet=='C' else 0)
X_sp = xm_sp - x
Y_sp = ym_sp - y
R_sp = sp.sqrt(X_sp**2 + Y_sp**2)
E_sp = k_sp * sp.ln(R_sp + Y_sp) - (k_sp - 1) * sp.ln(X_sp) if sujet in ['A', 'B'] else (k_sp * sp.ln(R_sp + X_sp) - (k_sp - 1) * sp.ln(Y_sp))

# --- La fonction E ---
E_sp = sp.simplify(E_sp)
display(sp.Eq(sp.Symbol('E(t)'), E_sp))

# Dérivée par rapport au temps
E_dt_sp = sp.diff(E_sp, t)

# RHS des EDOs
phi_x_sp = k_sp * v_sp * X_sp / R_sp
phi_y_sp = k_sp * v_sp * Y_sp / R_sp

# Substitution des dérivées dx/dt et dy/dt pour vérifier l'invariant
E_dt_sub_sp = sp.simplify(E_dt_sp.subs({
    sp.diff(x, t): phi_x_sp, 
    sp.diff(y, t): phi_y_sp
}))

display(sp.Eq(sp.Symbol("E'(t)"), E_dt_sub_sp))

\(\displaystyle E(t) = k \log{\left(t v + \sqrt{\left(t v - y{\left(t \right)}\right)^{2} + x^{2}{\left(t \right)}} - y{\left(t \right)} \right)} - \left(k - 1\right) \log{\left(- x{\left(t \right)} \right)}\)

\(\displaystyle E'(t) = 0\)

Simulations

Notation : on note \(\boldsymbol{\varphi}(t,x,y) = (\varphi_1(t,x,y), \varphi_2(t,x,y))\) la fonction qui apparaît dans le système de poursuite, c’est à dire

\( \varphi_1(t,x,y) = kv \dfrac{x_m(t) - x}{\sqrt{(x_m(t) - x)^2 + (y_m(t)-y)^2}}, \qquad \varphi_2(t,x,y) = kv \dfrac{y_m(t)-y}{\sqrt{(x_m(t) - x)^2 + (y_m(t)-y)^2}}, \)

ainsi le système s’écrit

\( \begin{cases} x'(t) = \varphi_1(t,x(t),y(t)), \\ y'(t) = \varphi_2(t,x(t),y(t)). \end{cases} \)

  • On définit d’abord les paramètres globaux : \(k\) le facteur de vitesse du chien, \(v\) la vitesse du maître, \(a\) la distance caractéristique (selon le sujet), \(t_0\) le temps initial, \(dt\) le pas de temps, \(\varepsilon\) la tolérance pour arrêter la simulation, et éventuellement max\_iter le nombre maximum d’itérations pour éviter les boucles infinies.

  • On initialise ensuite la position du chien à \(t=0\) et les fonctions de la trajectoire du maître en fonction du sujet.

  • On implémente ensuite les différentes méthodes numériques (Euler explicite, Heun ou Euler modifié) pour simuler la trajectoire du chien. Nota bene : on ne passera pas le vecteur tt à la fonction de calcul de la trajectoire du chien, car on ne connaît pas à l’avance le temps d’arrêt de la simulation (il dépend de la condition \(r(t) < \varepsilon\)). On peut simplement faire une boucle while qui s’arrête quand \(r(t) < \varepsilon\) (ou quand on atteint un nombre maximum d’itérations pour éviter les boucles infinies en cas d’erreur). Chaque fonction de calcul de la trajectoire du chien reçoit la fonction \(\boldsymbol{\varphi}\), le temps initial \(t_0\), le pas de temps \(h\), le vecteur initial yy0, et doit retourner les vecteurs tt et uu correspondants à la trajectoire du chien.

  • On lance la simulation pour les différentes méthodes et on trace les graphiques demandés.

Code 
# Paramètres globaux
# ===============================

k = 1.5        # facteur de vitesse du chien
v = 1.0        # vitesse du maître
dt = 0.01      # pas de temps
tol = 1e-2     # tolérance pour arrêter la simulation
a_val = 5      # constante de position du maître ou du chien selon le sujet

t0 = 0.0

# --- Configuration du Sujet ---
# ===============================
if sujet == 'A':
    # Maître vertical à droite, Chien à l'origine
    x0, y0 = 0.0, 0.0
    xm = lambda t: np.full_like(t, a_val)
    ym = lambda t: v * t
elif sujet == 'B':
    # Maître vertical à gauche, Chien à droite
    x0, y0 = a_val, 0.0
    xm = lambda t: np.zeros_like(t)
    ym = lambda t: v * t
elif sujet == 'C':
    # Maître horizontal en haut, Chien à l'origine
    x0, y0 = 0.0, 0.0
    xm = lambda t: v * t
    ym = lambda t: np.full_like(t, a_val)
elif sujet == 'D':
    # Maître horizontal en bas, Chien en haut
    x0, y0 = 0.0, a_val
    xm = lambda t: v * t
    ym = lambda t: np.zeros_like(t)

# phi du système
# ===============================
def pphi(t, zz):
    dx = xm(t) - zz[0]
    dy = ym(t) - zz[1]
    r = np.sqrt(dx**2 + dy**2)
    return k * v * np.array([dx / r, dy / r])


# ===============================
# Schémas
# ===============================

# On ne sait pas à l'avance combien d'itérations seront nécessaires, 
# car le critère d'arrêt dépend de la trajectoire.
# On cosntruit donc le tableau tt au fur et à mesure de la simulation.

# Schéma Euler explicite
# ===============================
def Euler_explicit_system(pphi, t0, h, yy0, max_iter, tol=1e-2):
    tt = np.array([t0])
    dog_traj = [np.array(yy0, dtype=float)]
    for i in range(max_iter):
        dog_traj.append(dog_traj[-1] + h * pphi(tt[-1], dog_traj[-1]))
        tt = np.append(tt, tt[-1] + h)
        # Condition d'arrêt
        dx = xm(tt[-1]) - dog_traj[-1][0]
        dy = ym(tt[-1]) - dog_traj[-1][1]
        r = np.sqrt(dx**2 + dy**2)
        if r < tol: break
    return tt, np.vstack(dog_traj)


# Schéma Heun
# ===============================
def Heun_system(pphi, t0, h, yy0, max_iter, tol=1e-2):
    tt = np.array([t0])
    dog_traj = [np.array(yy0, dtype=float)]
    for _ in range(max_iter):
        z_n = dog_traj[-1]
        t_n = tt[-1]
        # Prédicteur
        z_pred = z_n + h * pphi(t_n, z_n)
        # Correcteur (moyenne des pentes)
        z_next = z_n + h/2 * (pphi(t_n, z_n) + pphi(t_n + h, z_pred))
        dog_traj.append(z_next)
        tt = np.append(tt, t_n + h)
        # Condition d’arrêt
        dx = xm(tt[-1]) - z_next[0]
        dy = ym(tt[-1]) - z_next[1]
        r = np.sqrt(dx**2 + dy**2)
        if r < tol: break
    return tt, np.vstack(dog_traj)


# Schéma Euler Modifie
# ===============================
def Euler_modified_system(pphi, t0, h, yy0, max_iter, tol=1e-2):
    tt = np.array([t0])
    dog_traj = [np.array(yy0, dtype=float)]
    for _ in range(max_iter):
        z_n = dog_traj[-1]
        t_n = tt[-1]
        # Calcul du point intermédiaire
        z_mid = z_n + (h/2) * pphi(t_n, z_n)
        # Mise à jour avec la pente au point intermédiaire
        z_next = z_n + h * pphi(t_n + h/2, z_mid)
        dog_traj.append(z_next)
        tt = np.append(tt, t_n + h)
        # Condition d’arrêt
        dx = xm(tt[-1]) - z_next[0]
        dy = ym(tt[-1]) - z_next[1]
        r = np.sqrt(dx**2 + dy**2)
        if r < tol: break
    return tt, np.vstack(dog_traj)


# Schéma Crank-Nicolson 
# ===============================
# def Crank_Nicolson_system(pphi, t0, h, yy0, max_iter, tol=1e-2):
#     from scipy.optimize import fsolve
#     tt = np.array([t0])
#     dog_traj = [np.array(yy0, dtype=float)]
#     for _ in range(max_iter):
#         z_n = dog_traj[-1]
#         t_n = tt[-1]
#         # Estimation initiale (Euler explicite)
#         z_init = z_n + h * pphi(t_n, z_n)
#         # Résolution de l'équation implicite
#         func = lambda z_next: z_next - z_n - (h/2) * (pphi(t_n, z_n) + pphi(t_n + h, z_next))
#         z_next = fsolve(func, z_init)
#         dog_traj.append(z_next)
#         tt = np.append(tt, t_n + h)
#         # Condition d’arrêt
#         dx = xm(tt[-1]) - z_next[0]
#         dy = ym(tt[-1]) - z_next[1]
#         r = np.sqrt(dx**2 + dy**2)
#         if r < tol: break
#     return tt, np.vstack(dog_traj)


# Schéma odeint de SciPy 
# ===============================
def odeint_system(pphi, t0, h, yy0, max_iter, tol=1e-2):
    # On génère une grille de temps suffisamment longue pour que le chien puisse attraper le maître
    t_eval = np.arange(t0, t0 + max_iter * h, h)
    sol = odeint(lambda z, t: pphi(t, z), yy0, t_eval)    
    # On cherche le point où le chien attrape le maître à la tolérance fixée
    for i in range(len(t_eval)):
        dx = xm(t_eval[i]) - sol[i, 0]
        dy = ym(t_eval[i]) - sol[i, 1]
        if np.sqrt(dx**2 + dy**2) < tol:
            return t_eval[:i+1], sol[:i+1]
    return t_eval, sol


# ===============================
# Test 
# ===============================

# "max_iter" n'est pas nécéssaire, mais je préfère fixer une limite pour éviter 
# les boucles infinies en cas de problème dans les schémas numériques.
# Je connais le temps théorique de capture, je peux donc fixer un nombre d'itérations maximum en conséquence.
# Vous pouvez ne pas l'utiliser ou juste mettre un nombre très grand pour max_iter, 
# la condition d'arrêt dans les schémas fera le travail de stopper la simulation au bon moment.

t_F = a_val*k/(v*(k**2-1)) 
max_iter = int(t_F/dt)+1 


# Initialisation du chien
yy0 = np.array([x0, y0])

# Simulation avec les différents schémas
tt_euler, uu_euler = Euler_explicit_system(pphi, t0, dt, yy0, max_iter)
tt_heun, uu_heun   = Heun_system(pphi, t0, dt, yy0, max_iter)
tt_em, uu_em       = Euler_modified_system(pphi, t0, dt, yy0, max_iter)
# tt_cn, uu_cn       = Crank_Nicolson_system(pphi, t0, dt, yy0, max_iter)
tt_odeint, uu_odeint = odeint_system(pphi, t0, dt, yy0, max_iter)

# Affichage résultats
print("Résultats de la simulation :")
print(f"h = {dt}")
print("")

# Calcul de la distance parcourue par le chien pour chaque méthode 
# (pour vérifier que c'est cohérent avec la distance théorique)
# def distance_calc(tt, uu, method_name):
#     # Méthode 1 : somme des segments
#     d1 = sum(np.linalg.norm(np.diff(uu, axis=0), axis=1))
#     # Méthode 2 : vitesse constante * temps total
#     d2 = k * v * tt[-1]
#     print(f"{method_name} Temps capture : {tt[-1]:.3f}")
#     print(f"{method_name} Distance (somme segments) : {d1:.3f}")
#     print(f"{method_name} Distance (k*v*t)          : {d2:.3f}")
#     print("")

# distance_calc(tt_euler, uu_euler, "Euler")
# distance_calc(tt_heun, uu_heun, "Heun")
# distance_calc(tt_em, uu_em, "Euler Modifié")
# # distance_calc(tt_cn, uu_cn, "Crank-Nicolson")
# distance_calc(tt_odeint, uu_odeint, "odeint (LSODA)")


# Tracés
fig, axx = plt.subplots(figsize=(10,5))

axx.plot(xm(tt_euler), ym(tt_euler), label="Maître")
axx.plot(uu_euler[:,0], uu_euler[:,1],'.', label="Chien (Euler)")
axx.plot(uu_heun[:,0], uu_heun[:,1],'-.', label="Chien (Heun)")
axx.plot(uu_em[:,0], uu_em[:,1], '--', label="Chien (Euler Modifié)")
# axx.plot(uu_cn[:,0], uu_cn[:,1], label="Chien (Crank-Nicolson)")
axx.plot(uu_odeint[:,0], uu_odeint[:,1], ':', label="Chien (odeint)")

axx.set_xlabel("x")
axx.set_ylabel("y")
axx.grid(True)
axx.legend()

plt.tight_layout()
plt.show()
Résultats de la simulation :
h = 0.01

Comparons la conservation de l’énergie \(E\) du système à chaque itération avec les méthodes d’approximation.

Un invariant est une quantité qui reste constante malgré le mouvement. Ici, \(E(t)\) est une relation géométrique entre les positions relatives. Si le schéma est précis, la courbe de \(E(t)\) doit être une ligne parfaitement horizontale. Sa déviation est une mesure directe de l’erreur numérique du schéma d’approximation.

Code 
# Dans la suite, on limite le temps de la simulation à t_max << t_F 
# pour éviter les divergences numériques et se concentrer 
# sur la conservation de l'invariant.

t_max = 5

def E_exact_num(x, y, t_vals, a=a_val, v=v, k=k, sujet=sujet):
    X = (v * t_vals - x) if sujet in ['C', 'D'] else ( (a - x) if sujet=='A' else (0 - x) )
    Y = (v * t_vals - y) if sujet in ['A', 'B'] else ( (a - y) if sujet=='C' else (0 - y) )
    R = np.sqrt(X**2 + Y**2)
    E = k * np.log(R + Y) - (k-1) * np.log(np.abs(X)) if sujet in ['A', 'B'] else (k * np.log(R + X) - (k-1) * np.log(np.abs(Y)))
    return E


# ============================
# CALCUL DES SIMU POUR t < t_max
# ============================

def get_safe_data(tt, uu, t_max):
    mask = tt <= t_max
    return tt[mask], uu[mask]

tt_euler_safe, uu_euler_safe = get_safe_data(tt_euler, uu_euler, t_max)
tt_heun_safe, uu_heun_safe   = get_safe_data(tt_heun, uu_heun, t_max)
tt_em_safe, uu_em_safe       = get_safe_data(tt_em, uu_em, t_max)
tt_odeint_safe, uu_odeint_safe = get_safe_data(tt_odeint, uu_odeint, t_max)

# Calcul de E(t) pour chaque méthode
E_euler = E_exact_num(uu_euler_safe[:,0], uu_euler_safe[:,1], tt_euler_safe)
E_heun  = E_exact_num(uu_heun_safe[:,0], uu_heun_safe[:,1], tt_heun_safe)
E_em    = E_exact_num(uu_em_safe[:,0], uu_em_safe[:,1], tt_em_safe)
E_odeint = E_exact_num(uu_odeint_safe[:,0], uu_odeint_safe[:,1], tt_odeint_safe)

# Affichage des erreurs de conservation
print(f"--- Analyse de l'invariant (Sujet {sujet}) ---")
methods = [('Euler', tt_euler_safe, E_euler), 
           ('Heun', tt_heun_safe, E_heun), 
           ('Euler Modifié', tt_em_safe, E_em), 
           ('odeint', tt_odeint_safe, E_odeint)]

for label, tt_s, E_s in methods:
    err_rel = (E_s[-1] - E_s[0]) / E_s[0]
    print(f"{label:>15} : E(0) = {E_s[0]:.6f}, E({tt_s[-1]:.3f}) = {E_s[-1]:.6f}, err_rel = {err_rel:.2e}")

# ============================
# Tracé de la conservation
# ============================
plt.figure(figsize=(10,6))
plt.plot(tt_euler_safe,  E_euler,       label=r'$\mathcal{E}(t)$ Euler')
plt.plot(tt_heun_safe,   E_heun,        label=r'$\mathcal{E}(t)$ Heun')
plt.plot(tt_em_safe,     E_em, '--',    label=r'$\mathcal{E}(t)$ Euler Modifié')
plt.plot(tt_odeint_safe, E_odeint, ':', label=r'$\mathcal{E}(t)$ odeint')
plt.axhline(E_euler[0], color='black', linestyle='-', alpha=0.3, label='Valeur exacte théorique')

plt.xlabel(r'$t$')
plt.ylabel(r'$\mathcal{E}(t)$')
plt.title(f'Conservation de l\'invariant - Sujet {sujet}')
plt.grid(True, which='both', linestyle=':', alpha=0.5)
plt.legend()
plt.show()
--- Analyse de l'invariant (Sujet B) ---
          Euler : E(0) = 1.609438, E(5.000) = 1.622954, err_rel = 8.40e-03
           Heun : E(0) = 1.609438, E(5.000) = 1.609338, err_rel = -6.21e-05
  Euler Modifié : E(0) = 1.609438, E(5.000) = 1.609375, err_rel = -3.89e-05
         odeint : E(0) = 1.609438, E(5.000) = 1.609438, err_rel = -6.64e-09

REMARQUE 1 : champ de vecteur et lignes de courant

Dans l’étude des systèmes différentiels, le portrait de phase est un outil puissant pour visualiser le comportement qualitatif des solutions d’un système autonome, c’est-à-dire un système où le champ ne dépend pas explicitement du temps. Dans une telle configuration, les trajectoires réelles des solutions se confondent exactement avec les lignes de courant du champ. Cependant, notre problème de la poursuite constitue un système non autonome. Ici, \(\varphi_1\) et \(\varphi_2\) dépendent non seulement de la position \((x,y)\), mais aussi de la position du maître, laquelle évolue avec le temps \(t\). Si l’on trace le champ de vecteurs, celui-ci n’est pas statique mais évolue au cours du temps. Ce n’est donc pas un outil pertinent pour analyser qualitativement les trajectoires du chien.

Contexte et Interprétation

Cet exercice étudie un problème de cinématique de poursuite pure.

  • Le Maître se déplace à vitesse constante \(v\) sur une ligne droite (horizontale ou verticale selon le sujet).
  • Le Chien court à vitesse constante \(kv\) (\(k>1\)). À chaque instant, son vecteur vitesse pointe vers la position actuelle du maître.

Si \(k>1\), le chien rattrape son maître à l’instant \(t_F = \dfrac{ak}{v(k^2-1)}\).

Le point de rencontre \((x_F, y_F)\) théorique pour chaque sujet est

Sujet A : \((a, \frac{ak}{k^2-1})\), Sujet B : \((0, \frac{ak}{k^2-1})\), Sujet C : \((\frac{ak}{k^2-1}, a)\), Sujet D : \((\frac{ak}{k^2-1}, 0)\)

Remarque : on ne peut pas intégrer numériquement jusqu’à \(t_F\) car le dénominateur des équations s’annule au moment de la rencontre (distance \(r=0\)). C’est pour cette raison qu’on définit un seuil de distance \(\varepsilon\) pour arrêter le calcul.

REMARQUE 2 : solutions analytiques

Pour chacune des quatre configurations, il est possible d’obtenir une expression analytique de la trajectoire du chien \(y(x)\) ou \(x(y)\) et du temps \(t\) en fonction de la position.

La trajectoire du chien est pilotée par la distance \(u\) à la droite parcourue pas le maître (et non pas à la distance au maître qui est \(r(t)\)) :

  • maître vertical : \(u\) mesure l’écart horizontal et coïncide avec \(X(t)\)
  • maître horizontal : \(u\) mesure l’écart vertical et oïncide avec \(Y(t)\)

Cette variable permet de définir deux fonctions universelles pour toutes les configurations :

  • \(F_k(u)\) : position le long de la direction du maître
  • \(T_k(u)\) : temps nécessaire pour atteindre cette position

\( \begin{array}{|c|c|c|} \hline \text{Régime} & \text{Trajectoire } F_k(u) & \text{Temps } T_k(u) \\ \hline k \neq 1 & \displaystyle \frac{ak}{k^2-1} +\frac{k}{2}\left( \frac{a^{-1/k}}{k+1}u^{\frac{k+1}{k}} -\frac{a^{1/k}}{k-1}u^{\frac{k-1}{k}} \right) & \displaystyle \frac{ak}{v(k^2-1)} -\frac{1}{2v}\left( \frac{a^{-1/k}}{k+1}u^{\frac{k+1}{k}} +\frac{a^{1/k}}{k-1}u^{\frac{k-1}{k}} \right) \\ \hline k = 1 & \displaystyle \frac{a}{4}\left( \left(\frac{u}{a}\right)^2 -2\ln\left(\frac{u}{a}\right) -1 \right) & \displaystyle \frac{a}{4v}\left( 1-\left(\frac{u}{a}\right)^2 -2\ln\left(\frac{u}{a}\right) \right) \\ \hline \end{array} \)

La trajectoire exacte s’obtient alors en substituant \(u\) et la coordonnée dépendante selon le tableau suivant :

\( \begin{array}{|c|c|c|c|} \hline \text{Sujet} & \text{Configuration du Maître} & \text{Variable } u & \text{Solution exacte} \\ \hline \mathbf{A} & \text{Verticale en } x=a & u = X(t) = a-x & y(x)=F_k(a-x)\ \text{et}\ t(x)=T_k(a-x) \\ \hline \mathbf{B} & \text{Verticale en } x=0 & u = -X(t)=x & y(x)=F_k(x)\ \text{et}\ t(x)=T_k(x) \\ \hline \mathbf{C} & \text{Horizontale en } y=a & u = Y(t)=a-y & x(y)=F_k(a-y)\ \text{et}\ t(y)=T_k(a-y) \\ \hline \mathbf{D} & \text{Horizontale en } y=0 & u = -Y(t)=y & x(y)=F_k(y)\ \text{et}\ t(y)=T_k(y) \\ \hline \end{array} \)

Code 
from ipywidgets import interact, FloatSlider

# 1. Calcul de la borne maximale de l'affichage 
d_capture = (a_val * k) / (k**2 - 1) if k > 1 else a_val * 3 
borne_fixe = max(a_val, d_capture) + 1.0

# 2. Fonction de génération 
def get_exact_solution(sujet, a, k, v):
    u_vals = np.linspace(a, 1e-6, 100) 
    if k == 1:
        pos_long = (a/4) * ((u_vals/a)**2 - 2*np.log(u_vals/a) - 1)
        tau = (a/(4*v)) * (1 - (u_vals/a)**2 - 2*np.log(u_vals/a))
    else:
        term1 = (a**(-1/k) / (k+1)) * u_vals**((k+1)/k)
        term2 = (a**(1/k) / (k-1)) * u_vals**((k-1)/k)
        pos_long = (a*k)/(k**2-1) + (k/2)*(term1 - term2)
        tau = (a*k)/(v*(k**2-1)) - (1/(2*v))*(term1 + term2)
    
    if sujet == 'A': return a*np.ones_like(u_vals), v*tau, a-u_vals, pos_long, tau
    if sujet == 'B': return np.zeros_like(u_vals), v*tau, u_vals, pos_long, tau
    if sujet == 'C': return v*tau, a*np.ones_like(u_vals), pos_long, a-u_vals, tau
    if sujet == 'D': return v*tau, np.zeros_like(u_vals), pos_long, u_vals, tau

# 3. Pré-calcul des 4 trajectoires
data = {s: get_exact_solution(s, a_val, k, v) for s in ['A', 'B', 'C', 'D']}
t_fin = data['A'][4][-1] # Temps de capture (identique pour tous)

# 4. Fonction d'affichage 
def plot_all_simu(t_cur):
    fig, axes = plt.subplots(1, 4, figsize=(20, 5)) 
    sujets = ['A', 'B', 'C', 'D']
    
    for i, s in enumerate(sujets):
        ax = axes[i]
        xm_all, ym_all, xc_ex, yc_ex, tt_ex = data[s]
        
        idx = np.searchsorted(tt_ex, t_cur)
        if idx >= len(tt_ex): idx = len(tt_ex) - 1
        
        # --- TRAJECTOIRE CHIEN ---
        ax.plot(xc_ex[:idx+1], yc_ex[:idx+1], 'b-', lw=2, label='Chien')
        ax.plot(xc_ex[idx], yc_ex[idx], 'bo', markersize=6)
        
        # --- TRAJECTOIRE MAÎTRE ---
        if s == 'A': 
            ax.plot([a_val, a_val], [0, v*t_cur], 'r-', lw=2, label='Maître')
            xmt, ymt = a_val, v*t_cur
        elif s == 'B': 
            ax.plot([0, 0], [0, v*t_cur], 'r-', lw=2, label='Maître')
            xmt, ymt = 0, v*t_cur
        elif s == 'C': 
            ax.plot([0, v*t_cur], [a_val, a_val], 'r-', lw=2, label='Maître')
            xmt, ymt = v*t_cur, a_val
        elif s == 'D': 
            ax.plot([0, v*t_cur], [0, 0], 'r-', lw=2, label='Maître')
            xmt, ymt = v*t_cur, 0
        
        ax.plot(xmt, ymt, 'ro', markersize=6)
        ax.plot([xc_ex[idx], xmt], [yc_ex[idx], ymt], 'k--', alpha=0.2)
        
        ax.set_xlim(-0.5, borne_fixe)
        ax.set_ylim(-0.5, borne_fixe)
        ax.set_aspect('equal', adjustable='box')
        ax.set_title(f"Sujet {s}")
        ax.grid(True, linestyle=':', alpha=0.6)
        ax.legend(loc='upper right')
    
    plt.suptitle(f"Comparaison des trajectoires (t = {t_cur:.2f}s)", fontsize=14, y=1.05) 
    plt.show()

interact(plot_all_simu, 
         t_cur=FloatSlider(min=0, max=t_fin, step=0.1, value=0, description='Temps'));