TP EDO: résolution d'équations différentielles ordinaires¶

Marc BUFFAT, dpt mécanique, Université Lyon 1

In [1]:
%matplotlib inline
# bibliotheque
import numpy as np
from matplotlib import pyplot as plt
plt.rc('font', family='serif', size='14')
In [2]:
from validation.valide_markdown import test_markdown, test_code
from validation.validation import info_etudiant, bib_validation, test_function
bib_validation('cours','MGC2028L')
from Oscillation import test1,test2,test3,test4,test5,test6
from IPython.display import display, Markdown, clear_output
def printmd(string):
    display(Markdown(string))
# test si numero étudiant spécifier
try: NUMERO_ETUDIANT
except NameError: NUMERO_ETUDIANT = None 
if type(NUMERO_ETUDIANT) is not int :
    NOM, PRENOM, NUMERO_ETUDIANT = info_etudiant()
# parametres spécifiques
notebook_name='TP_EDO.ipynb'
_uid_  = NUMERO_ETUDIANT
np.random.seed(_uid_)
_U0_ = np.random.randint(10,30)
_omega0 = np.random.randint(1,15)
_m  = np.round(0.1+2*np.random.rand(),2)
_k  = _m*_omega0**2
###
printmd("**Etudiant** {} {}  id={}".format(NOM,PRENOM,NUMERO_ETUDIANT))

Etudiant NBGrader Course id=2592889824

EDO d'ordre 1¶

Rappel: résolution EDO avec RK2¶

On se propose de résoudre numériquement une EDO d'ordre 1, i.e trouver $y(t)$ t.q. $$ \frac{dy}{dt} = f(y,t) \mbox{ avec la condition initiale } y(0)= y_0 $$

Pour cela on se fixe un pas en temps $\Delta t$ et un nombre de pas en temps $n$, ce qui permet de définir un vecteur T contenant les temps où on calcule numériquement la solution, et qui est donnée par :

$$ T_i = i\Delta t \mbox{ pour } i=0,n $$

On calcule ensuite la solution numérique sous forme d'un vecteur Y correspondant à la valeur numérique de y(t) aux temps T

$$ Y[i] \approx y(t=T[i]) \mbox{ pour } i=1,n \mbox{ et } Y[0] = Y_0$$

Pour calculer les valeurs de Y, on utilise la méthode de Runge Kutta 2, dont les étapes sont données ci-dessous:

  1. la valeur pour i=0 est donnée par la CI : Y[0] = y0

  2. on calcule ensuite itérativement $Y_i$ à partir de la valeur précédente $Y_{i-1}$ pour i de 1 à n de la façon suivante:

    • sur l'intervalle $[t_{i-1}, t_i]$, on calcule la valeur moyenne $Y_m$ au milieu de l'intervalle par la méthode d'Euler:

    $$ Y_m = Y_{i-1} + \frac{(T_{i} - T_{i-1})}{2} \, f(Y_{i-1}, T_{i-1})$$

    • on utilise cette valeur moyenne pour calculer $Y_i$ par la méhode de Runge Kutta 2 (RK2):

    $$ Y_i = Y_{i-1} + (T_{i} - T_{i-1}) \, f(Y_m, T_m) \mbox{ avec } Tm = \frac{T_{i-1}+T_i}{2}$$

Cette méthode demande la création d'une fonction F(y,t) qui renvoie la valeur du second membre de l'EDO $f(y,t)$ pour des valeurs de y et t données. Le calcul de Y se fait ensuite à l'aide d'une boucle itérative en utilisant les relations ci-dessus.

Trajectoire verticale d'un projectile¶

On tire verticalement un projectile avec une vitesse initiale $v_0$, dont la valeur est affichée ci-dessous

La vitesse $v(t)$ est solution de l'équation différentielle suivante:

$$ \frac{dv}{dt} = -g - signe(v) C_f v^2 $$

avec $g=9.81$ et $C_f = 0.02$. La fonction signe(v) vaut +1 si $v>0$ et -1 sinon. Elle assure une force de frottement toujours opposée au mouvement.

L'objectif est de calculer quelle est la hauteur maximale $h_{max}$ atteinte par le projectile et à quel instant $t_{max}$.

Pour cela on résout d'abord numériquement cette équation pour calculer la vitesse. On détermine ensuite à quel instant cette vitesse s'annule, ce qui correspond à l'instant où le projectile atteint sa hauteur max. Ensuite on calcule sa trajectoire (en utilisant la méthode vue dans le TP précédent).

Pour résoudre numériquement cette EDO:

  1. écrire la fonction F(y,t) qui, pour des valeurs de y et de t données en argument, calcule le second membre de l'EDO étudiée. Pour être général, on conserve l'argument t, même si dans notre cas, F ne dépend pas explicitement de t. On définira les valeurs numériques des paramètres g et Cf dans la fonction.
  2. Dans la cellule suivante, mettre dans F0 le résultat de l'appel de la fonction au temps initial t=0. Afficher la valeur de F0 et son type. Commentez aussi le signe de F0.
In [3]:
print("Vitesse initiale v0={} m/s".format(_U0_))

Vitesse initiale v0=17 m/s
In [4]:
# fonction F(y,t)
## BEGIN SOLUTION
def F(y,t):
    g = 9.81
    Cf = 0.02
    if y < 0:
        Cf = -Cf
    return -g - Cf*y**2
## END SOLUTION
In [5]:
# verification
F0 = None
## BEGIN SOLUTION
Y0 = _U0_
F0 = F(Y0,0.0) 
print(Y0,F0)
if type(F0)== float and F0 < 0:
    print("OK")
#le projectile va décélerer
## END SOLUTION

17 -15.59
OK
In [6]:
# ne rien écrire dans cette cellule
### BEGIN HIDDEN TESTS
_Y1_ = np.random.rand()
_t1_ = np.random.rand()
assert(np.isclose(F(_Y1_,_t1_), -9.81 - 0.02*_Y1_**2))
### END HIDDEN TESTS

calcul de la solution numérique¶

On veut calculer la solution approchée dans un tableau V de taille N fixé. On se fixe aussi le temps Tmax de la simulation.

On part des valeurs de N et Tmax proposées

  • Créer un tableau numpy T contenant les valeurs de t équi-réparties entre 0 et Tmax en utilisant la fonction numpy linspace ( la fonction np.linspace(a,b,N) renvoie un tableau de N valeurs équi-réparties entre a et b)
  • En déduire la valeur de $\Delta t$ que l'on mettra dans la variable dt
  • Initialiser à zero le tableau numpy V
  • calculer itérativement la valeur de V en utilisant la méthode RK2.
  • afficher V, la valeur min et max de V, pour vérifier que V s'annule bien sur l'intervalle
  • changer si besoin les paramètres (Tmax et N) (p.e. pour augmenter la précision)
In [7]:
Tmax = 5
N = 101
T = None
V = None
dt = None
## BEGIN SOLUTION
T = np.linspace(0,Tmax,N)
V = np.zeros(N)
dt = T[1]-T[0]
V[0] = _U0_
for i in range(1,N):
    Vm   = V[i-1] + F(V[i-1],T[i-1])*dt/2.
    V[i] = V[i-1] + F(Vm,T[i-1]+dt/2.)*dt
print(V.min(), V.max())
print(V)
## END SOLUTION    

-20.27335613382166 17.0
[ 17.          16.23359959  15.49166832  14.77235966  14.07398098
  13.39497653  12.73391257  12.08946437  11.46040471  10.84559371
  10.24396985   9.65454195   9.07638205   8.50861897   7.95043263
   7.40104883   6.85973459   6.32579388   5.79856378   5.27741093
   4.76172824   4.25093194   3.74445876   3.24176331   2.74231568
   2.24559911   1.75110779   1.2583448    0.76682004   0.27604831
  -0.21445263  -0.70474133  -1.19433931  -1.68276895  -2.16955715
  -2.65423722  -3.13635061  -3.61544862  -4.09109405  -4.56286268
  -5.03034473  -5.49314616  -5.9508899   -6.40321686  -6.84978692
  -7.29027972  -7.72439531  -8.15185469  -8.57240023  -8.9857959
  -9.39182741  -9.79030226 -10.18104956 -10.56391988 -10.93878488
 -11.30553691 -11.66408851 -12.01437182 -12.3563379  -12.68995605
 -13.01521303 -13.33211225 -13.64067289 -13.94092911 -14.23292907
 -14.51673407 -14.79241767 -15.06006473 -15.31977054 -15.57163994
 -15.81578645 -16.05233141 -16.28140317 -16.50313629 -16.71767079
 -16.92515138 -17.12572683 -17.31954923 -17.50677342 -17.68755637
 -17.86205663 -18.03043382 -18.19284814 -18.34945992 -18.50042923
 -18.64591546 -18.78607701 -18.92107094 -19.05105271 -19.17617593
 -19.2965921  -19.41245041 -19.5238976  -19.63107773 -19.73413212
 -19.83319918 -19.92841435 -20.01990999 -20.10781534 -20.19225644
 -20.27335613]
In [8]:
# ne rien écrire dans cette cellule
### BEGIN HIDDEN TESTS
assert (T.size == N) and (V.size == N)
assert (np.isclose(dt, Tmax/(N-1)))
assert (np.isclose(V[1], V[0] + F(V[0]+F(V[0],0)*dt/2.,0.)*dt))
### END HIDDEN TESTS

tracer vitesse¶

  • Tracer la solution numérique V en fonction de T
  • Tracer aussi l'axe v=0
  • Mettre un titre, des labels sur les axes, et des légendes.
In [9]:
## BEGIN SOLUTION
plt.figure()
plt.plot(T,V,'b-',label="V")
plt.hlines(0.,T[0],T[-1],colors='r')
plt.title("Solution numérique")
plt.xlabel('t en sec.')
plt.ylabel('vitesse en m/s')
plt.legend()
## END SOLUTION
Out[9]:
<matplotlib.legend.Legend at 0x7f3e73b512b0>
In [10]:
# ne rien écrire dans cette cellule
### BEGIN HIDDEN TESTS
assert(test_code(notebook_name,'cell-verif12',['plt.plot','plt.title',
                                               'plt.xlabel','plt.ylabel','plt.legend']))
### END HIDDEN TESTS

Estimation du temps Tmax¶

On obtient une estimation du temps tmax $t_{max}$ en déterminant où dans le tableau des vitesses V la valeur de la vitesse devient, en valeur absolue, minimale.

Calculer cet indice k en utilisant la fonction numpy argmin: la fonction `np.argmin(X) renvoie les indices du ou des minimums du tableau X

Mettre dans la variable k l'indice et en déduire $t_{max}$ que l'on mettra dans Tmax

Affichez la valeur de Tmax et la valeur de la vitesse à cet instant (précision) et vérifiez par rapport au graphe précédent.

Optionnel: comment pourrait-on améliorer simplement la précision de ce calcul (sans le refaire) ?

Mettre dans un tableau Va les vitesses V entre 0 et k (inclus) et de même dans Ta les temps associés. Cela correspond à la partie ascendante de la trajectoire. Vérifier le type et la dimension de ces tableaux

In [11]:
k = None
Tmax = None
### BEGIN SOLUTION
k = np.argmin(np.abs(V))
Tmax = T[k]
print(k,Tmax,V[k])
# ou mieux : interpolation linéaire
#
Ta = T[0:k+1]
Va = V[0:k+1]
print("Verification : ",Ta.size==Va.size)
### END SOLUTION

30 1.5 -0.21445263153507987
Verification :  True
In [12]:
# ne rien écrire dans cette cellule
### BEGIN HIDDEN TESTS
assert((Tmax>0) and (Tmax<T[-1]))
assert(Ta.size == Va.size)
### END HIDDEN TESTS

Calcul de trajectoire¶

D'après l'étude précédente sur chaque intervalle de temps $[t_i,t_{i+1}]$, on connait la vitesse $u_i$ en $t_i$ et la vitesse $u_{i+1}$ en $t_{i+1}$. On peut donc en déduire la distance parcourue $d_i$ par l'objet entre $t_i$ et $t_{i+1}$ en utilisant la vitesse moyenne avec la relation :

$$ d_i = \frac{u_{i}+u_{i+1}}{2} (t_{i+1} - t_{i}) $$

A partir de la position $x_i$ à $t_i$, on peut alors en déduire la position $x_{i+1}$ à $t_{i+1}$

$$ x_{i+1} = x_i + d_i $$

Connaissant la position $x_0$ à l'instant initiale, on a donc une méthode itérative pour calculer la position du solide au cours du temps.

En créant un tableau numpy Xa pour stocker la position aux instants Ta, calculer les valeurs de Xa en utilisant cette méthode avec une boucle for .

Afficher les valeurs min et max de Xa et en déduire la hauteur $h_{max}$ à mettre dans Hmax

In [13]:
Xa = None
Hmax = 0
## BEGIN SOLUTION
Xa = np.zeros(Ta.size)
for i in range(Ta.size-1):
    Xa[i+1] = Xa[i] + (Ta[i+1]-Ta[i])*(Va[i+1]+Va[i])/2.
Hmax = Xa[-1]
print(min(Xa),max(Xa),Hmax)
## END SOLUTION

0.0 11.581217479245128 11.581217479245128
In [14]:
# ne rien écrire dans cette cellule
### BEGIN HIDDEN TESTS
assert(Xa.size == Ta.size)
assert(Hmax == Xa[-1])
### END HIDDEN TESTS

Etude du système simple masse ressort sans frottement et sans forcage¶

La loi de Newton appliquée au système masse ressort sans frottement donne:

\begin{equation} -k x = m \ddot{x} \end{equation}

avec les conditions initiales: $x(t=0)=x_0$ et $\dot{x}(t=0)=0$

Montrer que l'équation du mouvement peut s'écrire:

\begin{equation} \ddot{x} + \omega_0^2 x = 0 \end{equation}

On définiera le paramètre $\omega_0$ (et son interprétation) et on justifiera dans le compte rendu que dans notre cas cette équation différentielle du second ordre possède une solution analytique qui représente un mouvement harmonique simple (dont on donnera les caractéristiques période, amplitude):

$$x(t) = x_0 \cos(\omega_0 t)$$

Mise sous forme vectorielle¶

Pour résoudre numériquement une équation différentielle du second ordre, il faut tout d'abord la transformer en un ensemble de deux équations du premier ordre: dans ce cas, respectivement pour la position et la vitesse:

\begin{eqnarray} \dot{x} &=& v \nonumber\\ \dot{v} &=& -\omega_0^2 x \end{eqnarray}

Comme vu en cours, nous écrivons l'état du système comme un vecteur $Y(t)$,

\begin{equation} \mathbf{Y}(t) = \begin{bmatrix} x \\ v \end{bmatrix}, \end{equation}

et l'équation differentielle sous forme vectorielle:

\begin{equation} \dot{\mathbf{Y}} = \mathbf{F}(\mathbf{Y}(t),t) \mbox{ avec } \mathbf{F}(\mathbf{Y}(t),t) = \begin{bmatrix} v \\ -\omega_0^2 x \end{bmatrix}. \end{equation}

Écrire une fonction MasseRessort pour calculer le second membre $F(Y(t))$ de ce système d'équations différentielles masse-ressort.

    MasseRessort(Y,t)

On définira la valeur de $\omega_0$ du problème dans une variable globale omega0, pour pouvoir ensuite l'utiliser dans la fonction en le spécifiant dans la fonction

     global omega0

On vérifiera ensuite la fonction en l'appelant avec un vecteur d'état $Y_0$ correspondant à $x=0.1$ et $v=0.0$ pour t=0 et en mettant le résultat dans F0 . On affichera F0 et son type pour vérifier.

On remarque que l'on passe le temps t comme second argument de la fonction même si $Y$ ne dépend pas explicitement du temps, car on en aura besoin dans la suite pour traiter le cas avec un forcage

In [15]:
printmd("**Paramétres de l'étude:**  masse m={} [kg] raideur k={} [N/m]".format(_m,_k))

Paramétres de l'étude: masse m=1.11 [kg] raideur k=187.59 [N/m]

In [16]:
# definition de omega0 et de la fonction MasseRessort(Y,t)
omega0 = 0
### BEGIN SOLUTION
omega0 = np.sqrt(_k/_m)

def MasseRessort(Y,t):
    global omega0
    '''
    calcul le second membre du système masse ressort sans amortissement
    renvoie dY:  derivée du vecteur d'etat [v ,- ω*ω*x]
    '''
    global omega0
    dY = np.array([Y[1], -omega0**2*Y[0]])
    return dY
### END SOLUTION
In [17]:
# Verification: appel de la fonction pour la CI 
Y0 = None 
F0 = None 
### BEGIN SOLUTION 
Y0 = [0.1,0.0] 
F0 = MasseRessort(Y0,0) 
print(F0,type(F0)) 
### END SOLUTION

[  0.  -16.9] <class 'numpy.ndarray'>
In [18]:
# ne rien écrire dans cette cellule
### BEGIN HIDDEN TESTS
assert(test_function(MasseRessort,test1,omega0))
assert(test_code(notebook_name,'cell-verif1','MasseRessort('))
assert(np.abs(omega0/_omega0-1.0)<1.e-6)
assert(np.allclose(MasseRessort([1,1],0),[1,-omega0**2]))
## END HIDDEN TESTS

Méthode de Runge Kutta 2¶

Pour calculer la solution numérique du système d'EDO

$$\dot{\mathbf{Y}} = \mathbf{F}(\mathbf{Y}(t),t)$$

on utilise la méthode de Runge Kutta, mais sous forme matricielle puisque la solution est un vecteur (de dimension 2) et non plus un scalaire.

Pour cela on se fixe le temps de simulation $t_max$ et un nombre de pas en temps $n$, ce qui permet de définir un vecteur T contenant les temps où on calcule numériquement la solution, et qui est donnée par :

$$ T_i = i\Delta t \mbox{ pour } i=0,n $$

On calcule ensuite la solution numérique sous forme d'une matrice Y de dimension (n+1,2) qui pour chaque ligne i contient la valeur numérique de $\mathbf{Y}(T_i)$ au temps $T_i$

$$ Y[i,:] \approx \mathbf{Y}(t=T[i]) \mbox{ pour } i=1,n \mbox{ et } Y[0,:] = \mathbf{Y}_0$$

On calcule la valeur de ce vecteur solution $Y[i,:]$ avec la méthode numérique de Runge Kutta 2, qui s'écrit sous forme matricielle (en notant $t_m = \frac{t_i + t_{i-1}}{2}$) \begin{align} \mathbf{Y}_m & = \mathbf{Y}[i-1,:] + \frac{\Delta t}{2} \mathbf{F}(\mathbf{Y}[i-1,:],t_{i-1}) \\ \mathbf{Y}[i,:] & = \mathbf{Y}[i-1,:] + \Delta t \,\, \mathbf{F}(\mathbf{Y}_m,t_m). \end{align}

Fonction RK2¶

Pour appliquer cette méthode, on demande d'écrire une fonction RK2(F,T,Y0) qui prend pour arguments les paramètres du problème : la fonction second membre F de l'EDO, le vecteur T des pas en temps et la condition initiale Y0.

Cette fonction calcule la valeur numérique de la solution Y (qui est une matrice). En utilisant une boucle itérative, on calcule chaque ligne de Y en utilisant la méthode précédente. A la fin, la fonction renvoie
la valeur Y.

In [19]:
## Fonction RK2(F,T,Y0)
## BEGIN SOLUTION
def RK2(F,T,Y0):
    n = T.size
    Y = np.zeros((n,Y0.size))
    Y[0,:] = Y0
    for i in range(1,n):
        dt = T[i]-T[i-1]
        Ym     = Y[i-1,:] + F(Y[i-1,:],T[i-1]) * dt*0.5
        tm = (T[i]+T[i-1])/2
        Y[i,:] = Y[i-1,:] + F(Ym, tm) *dt
    return Y
## END SOLUTION
In [20]:
# ne rien écrire dans cette cellule
### BEGIN HIDDEN TESTS
assert(test_function(RK2,test6,MasseRessort))
### END HIDDEN TESTS

Vérification

En utilisant comme vecteur T=[0,0.1,0.2] , utilisez la fonction précédente pour calculer Y en ces instants. On vérifiera le type et la taille de Y.

In [21]:
# Verification
T = np.array([0,0.1,0.2])
Y0 = None
Y = None
### BEGIN SOLUTION
T = np.array([0,0.1,0.2])
Y0 = np.array([0.1,0.])
Y = RK2(MasseRessort,T,Y0)
print(Y,Y.shape,type(Y))
### END SOLUTION

[[ 0.1        0.       ]
 [ 0.0155    -1.69     ]
 [-0.1665975 -0.5239   ]] (3, 2) <class 'numpy.ndarray'>
In [22]:
# ne rien écrire dans cette cellule
### BEGIN HIDDEN TESTS
assert(Y.size == T.size*2)
assert(test_code(notebook_name,'cell-verif1b',['MasseRessort','RK2']))
#assert(test_function(iterationRK2,test3,MasseRessort))
### END HIDDEN TESTS

Paramétres et simulation¶

On étudie le mouvement pour une durée tfin égale à 10 périodes avec la valeur de la condition initiale Y0 donnée précédemment.

Pour étudier la précision de la méthode, on va faire le calcul en choisissant différents pas en temps, pour ensuite comparer les solutions calculées.

Pour cela on choisit de construire 3 tableaux de temps, pour calculer 3 tableaux solution:

  1. T1 avec 20 pas en temps par période, solution dans Y1
  2. T2 avec 40 pas en temps par période, solution dans Y2
  3. T3 avec 80 pas en temps par période, solution dans Y3

Calculer tfin et définir ces 3 tableaux T1,T2,T3 avec la fonction numpy linspace . Afficher le pas en temps pour chacun de ces tableaux et vérifier qu'ils sont bien multiples de 2.

Puis calculer la valeur de la numérique avec chacun de ces tableaux en utilisant la fonction RK2 précédente et en mettant la solution dans les matrices Y1, Y2 et Y3. Pour chaque cas, afficher la taille et le type de la solution calculée.

In [23]:
# parametres
Y0 = None
tfin = 0
T1 = None
T2 = None
T3 = None
### BEGIN SOLUTION
Y0 = np.array([0.1,0.0])
period = 2*np.pi/omega0
tfin = 10*period
T1 = np.linspace(0,tfin,10*20+1)
T2 = np.linspace(0,tfin,10*40+1)
T3 = np.linspace(0,tfin,10*80+1)
dt = period/100
print(T1[1]-T1[0],T2[1]-T2[0],T3[1]-T3[0])
Y1 = RK2(MasseRessort,T1,Y0)
print(Y1.shape,type(Y1))
Y2 = RK2(MasseRessort,T2,Y0)
print(Y2.shape,type(Y2))
Y3 = RK2(MasseRessort,T3,Y0)
print(Y3.shape,type(Y3))
### END SOLUTION

0.0241660973353061 0.01208304866765305 0.006041524333826525
(201, 2) <class 'numpy.ndarray'>
(401, 2) <class 'numpy.ndarray'>
(801, 2) <class 'numpy.ndarray'>
In [24]:
# ne rien écrire dans cette cellule
### BEGIN HIDDEN TESTS
assert(Y1.size == 2*T1.size)
assert(Y2.size == 2*T2.size)
assert(Y3.size == 2*T3.size)
### END HIDDEN TESTS

Analyse de la simulation¶

En utilisant matplotlib, comparer sur 2 figures la position $x(t)$ et la vitesse $v(t)$ calculées dans les 3 matrices Y1,Y2 et Y3.

In [25]:
# tracer
### BEGIN SOLUTION
fig = plt.figure(figsize=(12,10))
plt.subplot(1,2,1)
plt.plot(T1,Y1[:, 0], linewidth=2,  label='T1')
plt.plot(T2,Y2[:, 0], linewidth=2,  label='T2')
plt.plot(T3,Y3[:, 0], linewidth=2,  label='T3')
#plt.xlabel('Temps [s]')
plt.ylabel('$x$ [m]')
plt.title("Position : système masse ressort.\n")
plt.legend();
plt.subplot(1,2,2)
plt.plot(T1,Y1[:, 1], linewidth=2,  label='T1')
plt.plot(T2,Y2[:, 1], linewidth=2,  label='T2')
plt.plot(T3,Y3[:, 1], linewidth=2,  label='T3')
plt.xlabel('Temps [s]')
plt.ylabel('$u$ [m/s]')
plt.title("Vitesse : système masse ressort.\n")
plt.legend();
### END SOLUTION
In [26]:
# ne rien écrire dans cette cellule
### BEGIN HIDDEN TESTS
assert(test_code(notebook_name,'cell-verif22',['plt.plot','plt.title',
                                               'plt.xlabel','plt.ylabel','plt.legend']))
### END HIDDEN TESTS

Pour comparer les 3 solutions, on va calculer la norme de l'écart sur la position et la vitesse à partir des 2 tableaux d'écart: dY1 entre Y1 et Y3 et dY2 entre Y2 et Y3.

Attention Y1 , Y2 et Y3 non pas la même dimension, mais ont des points de calcul en commun. On calculera donc dY1 et dY2 en ces points uniquement.

Pour calculer la norme de l'écart sur la position et sur la vitesse à partir de dY1 et dY2, on pourra utiliser la fonction numpy suivante pour calculer la norme d'un vecteur X en normalisant par la dimension de X

np.linalg.norm(X)/X.size

Mettre les résultats dans ErrX1, ErrX2 , ErrU1, ErrU2 et afficher la valeur de ces écarts.

On calculera aussi la valeur maximum des écarts en valeur absolue sur la position dans ErrXmax1 et ErrXmax2.

Que peut on en conclure ?

In [27]:
dY1 = None
dY2 = None
ErrX1 = None
ErrU1 = None
ErrX2 = None
ErrU2 = None
ErrXmax1 = None
ErrXMax2 = None
## BEGIN SOLUTION
dY1 = Y1 - Y3[::4,:]
dY2 = Y2 - Y3[::2,:]
ErrX1 = np.linalg.norm(dY1[:,0])/T1.size
ErrU1 = np.linalg.norm(dY1[:,1])/T1.size
ErrX2 = np.linalg.norm(dY2[:,0])/T2.size
ErrU2 = np.linalg.norm(dY2[:,1])/T2.size
print("Erreur X ",ErrX1,ErrX2)
print("Erreur U ",ErrU1,ErrU2)
ErrXmax1 = np.max(np.abs(dY1[:,0]))
ErrXmax2 = np.max(np.abs(dY2[:,0]))
print("Erreur max ",ErrXmax1,ErrXmax2)
## END SOLUTION

Erreur X  0.0029687696468651274 0.00039909091186068473
Erreur U  0.03937391875118082 0.005215428233671283
Erreur max  0.10238352920032573 0.019184425427542458
In [28]:
# ne rien écrire dans cette cellule
### BEGIN HIDDEN TESTS
assert((dY1.size == Y1.size) and (dY2.size == Y2.size))
assert((ErrX1 > ErrX2 >0.0) and (ErrU1 > ErrU2 >0.0))
assert(ErrXmax1 > ErrXmax2 > 0.0)
### END HIDDEN TESTS

Comparaison de la période¶

A partir de la solution la plus précise Y3, on veut déterminer la période de la solution numérique, pour comparer avec la solution analytique.

Pour cela on remarque qu'avec la condition initiale, la vitesse est tout d'abord négative ou nulle, avant de devenir positive au bout d'une demi période. On va donc déterminer l'indice k de la première ligne du tableau Y3 correspondant à une vitesse calculée positive, puis à l'aide de T3 la valeur estimée de la période que l'on mettra dans periode_num, que l'on comparera avec la valeur exacte calculée dans periode

In [29]:
k = 0
periode_num = 0
periode = 0
## BEGIN SOLUTION
k = np.where(Y3[:,1]>0)[0][0]
periode_num = 2*T3[k]
periode = 2*np.pi/omega0
print("periode :", periode_num, periode)
## END SOLUTION

periode : 0.483321946706122 0.483321946706122
In [30]:
# ne rien écrire dans cette cellule
### BEGIN HIDDEN TESTS
assert((periode > 0.0) and (np.abs(periode-periode_num)<1.e-4))
### END HIDDEN TESTS

FIN du TP¶

In [55]:
# version
from platform import python_version,uname,platform
print("Systeme       :",uname())
print("OS            :",platform())
print("version python:",python_version())
print("version numpy :",np.__version__)

Systeme       : uname_result(system='Linux', node='l2-nbgrader', release='6.1.0-37-amd64', version='#1 SMP PREEMPT_DYNAMIC Debian 6.1.140-1 (2025-05-22)', machine='x86_64')
OS            : Linux-6.1.0-37-amd64-x86_64-with-glibc2.36
version python: 3.9.13
version numpy : 1.23.0
In [ ]: