Equation différentielle scalaire d'ordre 2

I) Objectifs
    1) Equation différentielle scalaire d'ordre 2
     2) Méthode
II) Description et évaluation des méthodes de résolution approchée
    1) Schémas numérique de résolution d'une équation différentielle scalaire d'ordre 2
    2) Algorithmes
III) Tests
    1) Le pendule
    2) Equation de Bessel

I) Objectifs

1) Equation différentielle scalaire d'ordre 2

  On veut résoudre le problème de Cauchy ( PC ) { y '' ( t ) = f ( t , y ( t ) , y ' ( t ) ) y ( a ) = y 0 et y ' ( a ) = y ' 0 , sur un segment [ a , b ] , c'est à dire calculer la fonction y : t y ( t ) définie sur [ a , b ] et vérifiant le système (PC) qui est composé d'une équation différentielle scalaire d'ordre 2 ( l'ED y'' (t) = f (t, y (t), y'(t)) ) et de conditions initiales ( y ( a ) = y 0 et y ' ( a ) = y ' 0 ).
  Dans ce système (PC), f est une fonction  définie sur une partie de R 3 et à valeurs dans R .

2) Méthode

  On ramène la résolution de cette équation scalaire d'ordre 2 à la résolution d'un système différentiel:

Pour résoudre ( PC2 ) { y '' ( t ) = f ( t , y ( t ) , y ' ( t ) ) y ( a ) = y 0 et y ' ( a ) = y ' 0 , où y : [ a , b ] R est la fonction inconnue:

    (1) On pose Y ( t ) = ( y ( t ) y ' ( t ) ) ,   Y 0 = ( y 0 y ' 0 ) et g ( t , Y ) = g ( t , x , y ) = ( y f ( t , x , y ) ) . Alors Y ' ( t ) = ( y ' ( t ) y '' ( t ) ) .
    (2) On ramène (PC2) au problème de Cauchy : ( PC1 ) { Y ' ( t ) = g ( t , Y ( t ) ) Y ( a ) = Y 0
    (3) On résout (PC1) et on récupère la fonction inconnue y ( t ) comme abcisse de Y ( t ) .

En effet, (PC1) ⇔ { Y ' ( t ) = g ( t , Y ( t ) ) Y ( a ) = Y 0 { ( y ' ( t ) , y '' ( t ) ) = g ( t , y ( t ) , y ' ( t ) ) ( y ( a ) , y ' ( a ) ) = ( y 0 , y ' 0 ) { ( y ' ( t ) , y '' ( t ) ) = ( y ' ( t ) , f ( t , y ( t ) , y ' ( t ) ) ) ( y ( a ) , y ' ( a ) ) = ( y 0 , y ' 0 ) ⇔   { y '' ( t ) = f ( t , y ( t ) , y ' ( t ) ) y ( a ) = y 0 et y ' ( a ) = y ' 0 ⇔ (PC2)

II) Description et évaluation des méthodes de résolution approchée

  On utilise encore, comme pour les ED scalaires, la méthode d'Euler et la méthode de Runge Kutta d'ordre 4.

1) Schémas numérique de résolution d'une équation différentielle scalaire d'ordre 2

Pour calculer sur [ a , b ] une solution approchée z : t z ( t ) de la solution y : t y ( t ) du problème de Cauchy ( PC ) { y '' ( t ) = f ( t , y ( t ) , y ' ( t ) ) y ( a ) = y 0 et y ' ( a ) = y ' 0 avec n points:
    (0) On pose g ( t , x , y ) = ( y f ( t , x , y ) ) et Y 0 = ( y 0 y ' 0 ) .

    (1) On pose   t k = a + k h avec h = b - a n pour k { 0 , ... , n } et   σ n = ( t 0 , t 1 , ... , t n ) la subdivision régulière d'ordre n du segment [ a , b ] .

    (2a) (Méthode d'Euler) On définit la suite ( Z 0 , Z 1 , ... , Z n )   par récurrence par:  (les Z k sont des points de R 2 )
Z 0 = Y 0 et k { 0 , ... , n - 1 } , Z k + 1 = Z k + h g ( t k , Z k )
    
    (2b) (Méthode RK4) On définit la suite ( Z 0 , Z 1 , ... , Z n )   par récurrence par:  (les Z k et les X k sont des points de R 2 )  
Z 0 = Y 0 et k { 0 , ... , n - 1 } , Z k + 1 = Z k + 1 6 ( X 1 + 2 X 2 + 2 X 3 + X 4 ) avec: { X 1 = h g ( t k , Z k ) X 2 = h g ( t k + h 2 , Z k + X 1 2 ) X 3 = h g ( t k + h 2 , Z k + X 2 2 ) X 4 = h g ( t k + h , Z k + X 3 )

    (3) La fonction z est alors  la fonction continue et affine par morceaux sur [ t k , t k + 1 ] (pour k { 0 , ... , n - 1 } ) définie par:
z ( t k ) = x k pour k { 0 , ... , n } , avec Z k = ( x k y k ) .

2) Algorithmes

  En plus du calculs des schémas numériques, on propose quatre algorithmes permettant de voir l'influence du choix du pas et de la continuité par rapport aux conditions initiales.

from math import *
import numpy as np, pylab as plt

a) Pour commencer, le calcul des schémas numériques Euler_ED2() et RK4_ED2() par les deux méthodes Euler puis RK4 pour  le problème de Cauchy ( PC ) { y '' ( t ) = f ( t , y ( t ) , y ' ( t ) ) y ( a ) = y 0 et y ' ( a ) = y ' 0 sur [a, b] avec n points: (on calcule le schéma vectoriel associé)

def Euler_ED2(f, Y0, a, b, n):
    """ Schéma numérique E pour le problème de Cauchy y'(t) = f(t, y(t), y'(t))
    avec [y(a), y'(a)] = Y0 sur [a,b], avec n points par la méthode d'Euler"""
# On se ramène au système différentiel Y'(t) = g(t, Y(t)) avec Y(t)=[y(t),y'(t)]
    def g(t, X):
        return np.array([X[1], f(t, X[0], X[1])])
    h, s, E = (b-a) / n, np.linspace(a, b, n + 1), np.zeros((n+1, 2))
    E[0] = Y0
    for k in range(1, n + 1):
        E[k] = E[k-1] + h * g(s[k-1], E[k - 1])
    return E

def RK4_ED2(f, Y0, a, b, n):
    """ Schéma numérique R pour le problème de Cauchy y'(t) = f(t, y(t), y'(t))
    avec [y(a), y'(a)] = Y0 sur [a,b], avec n points par la méthode RK4"""
# On se ramène au système différentiel Y'(t) = g(t, Y(t)) avec Y(t)=[y(t),y'(t)]
    def g(t, X):
        return np.array([X[1], f(t, X[0], X[1])])
    h, s, R = (b-a) / n, np.linspace(a, b, n + 1), np.zeros((n+1, 2))
    R[0] = Y0
    for k in range(1, n + 1):
        t, V = s[k-1], R[k-1]
        X1 = h * g(t, V)
        X2 = h * g(t + h/2, V + X1/2)
        X3 = h * g(t + h/2, V + X2/2)
        X4 = h * g(t + h, V + X3)
        R[k] = V + (X1 + 2*X2 + 2*X3 + X4)/6
    return R

b) Le premier algorithme ED2_pedago() permet de tracer sur un même graphe pour le problème de Cauchy ( PC ) { y '' ( t ) = f ( t , y ( t ) , y ' ( t ) ) y ( a ) = y 0 et y ' ( a ) = y ' 0 donné:
    • La solution exacte
    • la solutions approchée calculée par la méthode d'Euler ou RK4
    • les courbes intégrales passant par les points calculés de la solution approchée

def ED2_pedago(f, Y0, a, b, n, y, ED, methode):
    """ Pour le problème de Cauchy y'(t) = f(t, y(t), y'(t))
    avec [y(a), y'(a)] = Y0 sur [a,b],
    on trace sur un même graphe:
        * la solution approchée z par la méthode d'Euler ou RK4 avec n points
    calculés
        * la solution exacte du problème de cauchy
        * les courbes intégrales (les solutions exactes) passant par les points
     calculés de la solution approchée.

      La courbe intégrale (solution exacte) de y'(t) = f(t, y(t), y'(t))
      avec [y(c), y'(c)] = D est t -> y(t, c, D).
     methode = 0 pour Euler ou 1 pour RK4.
     ED = chaîne de caractères pour décrire le SD  """
# Divers
    plt.rcParams['figure.figsize'] = 12, 9   # fenêtre de tracé
    plt.grid(True)      # Grille
# Calcul et tracé de la solution exacte y
    x_exact = np.linspace(a, b, 200)
    y_exact = [y(x, a, Y0) for x in x_exact]
    plt.plot(x_exact, y_exact, linewidth = 2.0, color = "black")
# Calcul et tracé de la solution approchée Z() du système différentiel
    h, s = (b-a) / n, np.linspace(a, b, n + 1)
    if methode == 0:
        meth, col = "Euler", "blue"
        Z = Euler_ED2(f, Y0, a, b, n)
    elif methode == 1:
        meth, col = "RK4", "red"
        Z = RK4_ED2(f, Y0, a, b, n)
    else:
        print("""methode doit être 0 (Euler) ou bien 1 (RK4) """)
        return
    x = Z.transpose()[0]
    plt.plot(s, x, linewidth = 2.0, color = col, marker = ".", markersize = 10)
# Calcul et tracé des courbes intégrales passant par les points calculés
    for k in range(n + 1):
        y1 = [y(t, s[k], Z[k]) for t in x_exact]
        plt.plot(x_exact, y1, color = "grey")
# On "repasse" la solution exacte (altérée par les courbes intégrales)
    plt.plot(x_exact, y_exact, linewidth = 2.0, color = "black")
# Labels, légende, titre, enregistrement
    plt.xlabel("t")
    plt.ylabel("y(t)")
    plt.legend(("Solution exacte", "Solution approchée",
                "Courbes intégrales"), 'upper center')
    titre = "{0} sur [{1:.4}, {2:.4}] avec n = {3} pour {4} et \
[y({1}), y'({1})] = {5}".format(meth, float(a), float(b), n, ED, Y0)
    plt.title(titre)
    plt.savefig("ED2_" + meth + "_" + str(n) + ".png", dpi = 100)
    plt.show()

Par exemple, pour le problème de Cauchy (PC) { y '' ( t ) = - y ( t ) y ( 0 ) = 1 et y ' ( 0 ) = 0 sur [0, 12.6] avec 16 points par les deux méthodes Euler puis RK4: (la solution exacte est y ( t ) = cos ( t ) )

def P1():
    def f(t, y, dy):
        return -y
    def y(t, c, D):
        """ Solution exacte y(t) = t(t, c, D) telle que [y(c), y'(c)] = D"""
        return D[0] * cos(t - c) + D[1] * sin(t - c)
    ED = """y" = -y """
    ED2_pedago(f, [1, 0], 0, 12.6, 16, y, ED, 0)
    ED2_pedago(f, [1, 0], 0, 12.6, 16, y, ED, 1)
##P1()

  Comme toujours, le schéma de Runge Kutta est bien meilleur:

numerique_ED2_1.pngnumerique_ED2_2.png

c) Le second algorithme ED2_graphes_pas() permet de tracer les solutions approchées par les méthodes d'Euler ou (et) RK4 pour un problème de Cauchy ( PC ) { y '' ( t ) = f ( t , y ( t ) , y ' ( t ) ) y ( a ) = y 0 et y ' ( a ) = y ' 0 , avec un nombre n de points dans une liste plage_n donnée, dans une série de sous graphes. On peut voir ainsi l'effet du choix du nombre de points calculés (et donc du pas) sur la qualité de l'approximation.

def ED2_graphes_pas(f, Y0,  a, b, plage_n, y, ED, methode):
    """ y est la solution exacte du problème de Cauchy
    y'(t) = f(t, y(t), y'(t)) avec [y(a), y'(a)] = Y0 sur [a,b]
    e et r sont les solutions approchées par les méthodes d'Euler et RK4
    avec un pas de h = (b-a)/n pour n dans une liste plage_n.
    Tracé des trois graphes y, e, r dans des sous-tracés.
    ED = chaîne de caractères pour décrire l'ED.
    methode = 0 (Euler), 1 (RK4) ou 2 (pour les deux)."""
# Divers
    plt.rcParams['figure.figsize'] = 12, 9   # fenêtre de tracé carrée
    plt.subplots_adjust(bottom = 0.05, top = 0.9, left = 0.03, right = 0.97 )
    plt.subplots_adjust(hspace = 0.15, wspace = 0.1)
    plage_x = np.linspace(a, b, 200)
    nb_lignes = (2 + len(plage_n)) // 3
# Calcul de la solution exacte y avec 200 points
    x_exact = np.linspace(a, b, 200)
    y_exact = [y(t) for t in x_exact]
# Calculs des graphes
    for i, n in enumerate(plage_n):
        plt.subplot(nb_lignes, 3, i + 1)
# Tracé de la solution exacte
        plt.plot(x_exact, y_exact, linewidth = 2.0, color = "black")
        plt.grid(True)                          # Grille
# Calcul des solutions approchées E et R du SD par Euler et RK4 avec n points
        h, s = (b-a) / n, np.linspace(a, b, n + 1)
        E, R = Euler_ED2(f, Y0, a, b, n), RK4_ED2(f, Y0, a, b, n)
# Légendes et tracés
        if methode == 0:
            meth = "Euler"
            legende = ("Exacte","Euler")
            e = E.transpose()[0]
            plt.plot(s, e, linewidth = 2.0, color = "blue")
        elif methode == 1:
            meth = "RK4"
            legende = ("Exacte", "RK4")
            r = R.transpose()[0]
            plt.plot(s, r, linewidth = 2.0, color = "red")
        elif methode == 2:
            meth = "Euler et RK4"
            legende = ("Exacte","Euler", "RK4")
            e = E.transpose()[0]
            plt.plot(s, e, linewidth = 2.0, color = "blue")
            r = R.transpose()[0]
            plt.plot(s, r, linewidth = 2.0, color = "red")
        else:
            legende, meth = ("Exacte",""), "Méthode à choisir..."
        if i == 0:
            plt.legend(legende, 'upper left')
        plt.title("n = {0} ; pas = {1:.4}".format(n, h))
# Titre, enregistrement
    titre = "{0} pour {1} sur [{2:.4}, {3:.4}] avec \
[y({2}), y'({2})] = {4}".format(meth, ED, float(a), float(b), Y0)
    plt.suptitle(titre, fontsize = 18)
    plt.savefig("ED2_" + meth + str(plage_n) + ".png", dpi = 100)
    plt.show()

Par exemple, pour le problème de Cauchy (PC) { y '' ( t ) = - y ( t ) y ( 0 ) = 1 et y ' ( 0 ) = 0 sur [0, 12.6] avec 16 points par les deux méthodes Euler puis RK4: (la solution exacte est y ( t ) = cos ( t ) )

def P2():
    def f(t, y, dy):
        return -y
    def y(t):
        """ Solution exacte"""
        return cos(t)
    ED = """ y" = - y """
    ED2_graphes_pas(f, [1,0], 0, 12.6, [10, 20, 40, 80, 160, 320], y, ED, 0)
##P2()

  Les graphes suivants sont explicites:

numerique_ED2_3.pngnumerique_ED2_4.png

d) Le troisième algorithme ED2_erreur_maxi() permet de calculer, pour le pas h = b - a n , les quantités E ( h ) = max k { 0 , ... , n } | e ( t k ) - y ( t k ) | et RK ( h ) = max k { 0 , ... , n } | r ( t k ) - y ( t k ) | les écarts maximum entre la solution exacte y et les solutions approchées e ( Euler) ou r (RK4) aux points de la subdivison σ n de l'intervalle [ a , b ]   pour n dans une liste donnée pour le problème de Cauchy ( PC ) { y '' ( t ) = f ( t , y ( t ) , y ' ( t ) ) y ( a ) = y 0 et y ' ( a ) = y ' 0 .
  En faisant ces calculs on peut voir l'ordre des méthodes numériques : toujours 1 pour Euler et 4 pour Runge Kutta.

def ED2_erreur_maxi(f, Y0,  a, b, liste_n, y, ED):
    """ y est la solution exacte du problème de Cauchy
    y'(t) = f(t, y(t), y'(t)) et [y(a), y'(a)] = Y0 sur [a,b]
    e et r sont les solutions approchées par les méthodes d'Euler et RK4 avec un
    pas de h = (b-a)/n. On calcule pour n dans liste_n le pas h, les écarts
    maxis E(h) et RK(h) aux points de la subdivision entre e ou r et y et les
    rapports E(h) / h et RK(h) / h^4
    ED = chaîne de caractères pour décrire l'ED  """
    print()
    print("Pour l'ED {0} avec [y({1}), y'({1})] = {3} sur [{1}, {2}]".format(
          ED, float(a), float(b), Y0))
    print()
    for n in liste_n:
        sol, h = [y(a + k *(b - a)/n) for k in range(n + 1)], (b - a)/n
        e = Euler_ED2(f, Y0, a, b, n).transpose()[0]
        r = RK4_ED2(f, Y0, a, b, n).transpose()[0]
        maxE = max(abs(sol - e))
        maxRK = max(abs(sol - r))
        print("n ={0:8} | h = {1:8.3} | E(h) = {2:8.3} |  RK(h) = {3:8.3} ".
              format(n, h, maxE, maxRK))

Par exemple, pour le problème de Cauchy (PC) { y '' ( t ) = - y ( t ) y ( 0 ) = 1 et y ' ( 0 ) = 0 sur [0, 12.6] :

def P3():
    def f(t, y, dy):
        return -y
    def y(t):
        """ Solution exacte"""
        return cos(t)
    ED = """ y" = - y """
    plage_n = [10, 100, 1000, 10000]
    ED2_erreur_maxi(f, [1, 0], 0, 12.6, plage_n, y, ED)
P3() #renvoie

Pour l'ED  y” = - y  avec [y(0.0), y'(0.0)] = [1, 0] sur [0.0, 12.6]

n =      10 | h =     1.26 | E(h) = 1.07e+02 |  RK(h) =    0.207
n =     100 | h =    0.126 | E(h) =      1.2 |  RK(h) = 2.34e-05
n =    1000 | h =   0.0126 | E(h) =   0.0826 |  RK(h) = 2.32e-09
n =   10000 | h =  0.00126 | E(h) =  0.00797 |  RK(h) = 1.45e-12

  Lorsque n augmente d'un facteur 10, E(h) diminue d'un facteur 10 pour Euler (méthode d'ordre 1) et 10 4 pour RK4 (méthode d'ordre 4).

e)  On se place maintenant dans des situations où l'on ne sait pas forcément calculer les solutions exactes.

  On modifie un peu le calcul du schéma numérique RK_ED2() pour n'en conserver qu'une partie bornée (tout éventuellement)

def RK4_ED2_borne(f, Y0, a, b, n, maxi):
    """ Schéma numérique R pour le problème de Cauchy y'(t) = f(t, y(t), y'(t))
    avec [y(a), y'(a)] = Y0 sur [a,b], avec n points par la méthode RK4.
    On s'arrête si |r(t)| > maxi."""
# On se ramène au système différentiel Y'(t) = g(t, Y(t)) avec Y(t)=[y(t),y'(t)]
    def g(t, X):
        return np.array([X[1], f(t, X[0], X[1])])
    h, s, R = (b-a) / n, np.linspace(a, b, n + 1), np.zeros((n+1, 2))
    R[0] = Y0
    for k in range(1, n + 1):
        t, V = s[k-1], R[k-1]
        X1 = h * g(t, V)
        X2 = h * g(t + h/2, V + X1/2)
        X3 = h * g(t + h/2, V + X2/2)
        X4 = h * g(t + h, V + X3)
        R[k] = V + (X1 + 2*X2 + 2*X3 + X4)/6
        if abs(R[k][0]) > maxi:
            break
    return (np.array(s[:k]), np.array(R[:k]))

  Le quatrième algorithme  ED2_RK() trace sur l'intervalle [ a , b ] la solution approchée RK4 du problème de Cauchy ( PC ) : { y '' ( t ) = f ( t , y ( t ) , y ' ( t ) ) y ( a ) = y 0 et y ' ( a ) = y ' 0   avec n points en ne conservant que la partie de solution (éventuellement toute la solution) pour lesquelles |y(t)|≤maxi.

def ED2_RK(f, Y0, a, b, n, maxi, ED):
    """ graphe de la solution approchée r avec n points du problème de Cauchy
    y'(t) = f(t, y(t), y'(t)) et [y(a), y'(a)] = Y0 sur [a,b] avec RK4.
    et donc un pas de h = (b-a)/n . On s'arrête si |r(t)| > maxi.
    ED = chaîne de caractères pour décrire l'ED"""
# Divers
    plt.rcParams['figure.figsize'] = 12, 9   # fenêtre de tracé
    plt.grid(True)                          # Grille
# Calcul et tracé
    h = (b-a) / n
    (s, R) = RK4_ED2_borne(f, Y0, a, b, n, maxi)
    r = R.transpose()[0]
    plt.plot(s, r, linewidth = 2.0, color = "red")
# Titre, labels et enregistrement
    titre = "{0} et [y({1}), y'({1})] = {3} sur [{1:4.3} , {2:4.3}] \
avec pas = {4:.4}".format(ED, float(a), float(b), Y0, h)
    plt.title(titre)
    plt.xlabel("t")
    plt.ylabel("y(t)")
    plt.savefig("ED2_ " + str(n) + ".png", dpi = 100)
    plt.show()

Par exemple pour (PC) { y '' ( t ) = - cos ( y ( t ) ) - 0.3 * sin ( y ' ( t ) ) y ( 0 ) = 1 et y ' ( 0 ) = 0 sur [0, 30]

def P4():
    def f(t, y, dy):
        return cos (y) - 0.3*sin(dy)
    ED = """ y" = - cos(y) - 0.3*sin(y') """
    ED2_RK(f, [1, 0], 0, 30, 1000, 10, ED)
P4()

numerique_ED2_5.png

f)  Enfin, la fonction ED2_graphes_CI() ci-dessous trace sur l'intervalle [ a , b ] l'ensemble des solutions numériques RK4 des problèmes de Cauchy ( P Y0 ) : { y '' ( t ) = f ( t , y ( t ) , y ' ( t ) ) y ( a ) = y 0 et y ' ( a ) = y ' 0 avec Y 0 liste_Y0 , et n points en ne conservant que les parties de solutions (éventuellement toute la solution) pour lesquelles |y(t)|≤maxi
  Ces courbes permettent de voir comment évoluent les solutions du problème de Cauchy en fonction des conditions initiales.

def ED2_graphes_CI(f, liste_Y0, a, b, n, maxi, SD):
    """ graphes des  solutions r approchées par RK4 du problème de Cauchy
    y'(t) = f(t, y(t), y'(t)) et [y(a), y'(a)] = Y0 sur [a,b] avec un pas de
    (b-a)/n et Y0 dans liste_Y0. On s'arrête si |R(t)| > maxi"""
# Divers
    plt.rcParams['figure.figsize'] = 12, 9   # fenêtre de tracé carrée
    plt.grid(True)                          # Grille
# Calcul et tracés
    h = (b-a) / n
    for Y0 in liste_Y0:
        (s, R) = RK4_ED2_borne(f, Y0, a, b, n, maxi)
        r = R.transpose()[0]
        plt.plot(s, r)
# Titre, labels et enregistrement
    titre = "{4} : {0} solutions sur [{1:4.3} , {2:4.3}] \
avec pas = {3:.4}".format(len(liste_Y0), float(a), float(b), h, SD)
    plt.title(titre)
    plt.xlabel("t")
    plt.ylabel("y(t)")
    plt.savefig("ED2_CI_ " + str(n) + ".png", dpi = 100)
    plt.show()

Par exemple pour y '' ( t ) = - cos ( y ( t ) ) - 0.3 * sin ( y ' ( t ) )   sur [0, 30]

def P5():
    def f(t, y, dy):
        return cos (y) - 0.3*sin(dy)
        return -2 * y - dy
    ED = """ y" = - cos(y) - 0.3*sin(y') """
    liste_Y0 = [[1, a] for a in  np.linspace(-3, 3, 81)]
    ED2_graphes_CI(f, liste_Y0, 0, 25, 500, 10, ED)
P5()

numerique_ED2_6.png

III) Tests

1) Le pendule

L'angle y(t) (t en s, t ∈ [0, T]) que fait un pendule non amorti de longueur 1 m par rapport à la verticale lâché à t=0 avec un angle a  est la solution du problème de Cauchy

( P ( a ) ) : { y '' = - 9.81 sin ( y ) ( y ( 0 ) , y ' ( 0 ) ) = ( a , 0 )   . La solution approchée habituelle est   y 1 ( t ) = a cos ( 9.81 t ) pour de petites oscillations (a petit).

Q1) Ecrire une fonction Pendule(a, T, n) qui:
    (1) définit l'approximation mathématique usuelle y 1 ( t ) = a cos ( 9.81 t )
    (2) calcule sur la subdivision s de n points de [0,T] la solution “exacte” r dui problème de Cauchy P ( a ) par RK4
    (3) trace sur un même graphe les fonctions y1() et r() avec légendes, labels et titres.

Constater que le “faux” pendule prend du retard sur le “vrai”. Par exemple

numerique_ED2_7.png

Avec a = 0.1, au bout de combien de temps le déphasage entre les deux mouvements est-il d'environ 1 s ?

2) Equation de Bessel

L'équation de Bessel E a est l'ED : t 2 y '' ( t ) + t x ' ( t ) + ( t 2 - a 2 ) x ( t ) = 0   (où a R + * )

Q2) Ecrire une fonction Bessal(a) calculant l'ensemble des courbes intégrales de E a sur l'intervalle I = [ 1 100 , 30 ] avec les conditions initiales s = [[0, x] for x in np.linspace(-2, 2, 61)]. par exemple, avec a = 0.5:

numerique_ED2_8.png

© 2014 - Eric Obermeyer      Powered by      Corrigé