Soutenez-nous
Publicité
+ Répondre à la discussion
Affichage des résultats 1 à 3 sur 3
  1. #1
    Invité de passage
    Inscrit en
    mai 2008
    Messages
    6
    Détails du profil
    Informations forums :
    Inscription : mai 2008
    Messages : 6
    Points : 2
    Points
    2

    Par défaut Runge Kutta 4 à pas variables

    Bonjour à tous,
    Je voudrais resoudre les équations du pendule double grâce à la méthode de Runge Kutta à pas variable.
    Mais je n'arrive pas à definir et à utiliser correctement le pas pour chaque étapes.
    je m'explique: j'ai d'abord créé une boucle pour un pas constant ici 0.001 qui marche très bien.
    L'utilisation d'un pas variable me ferai gagner du temps de calcul qui me permettrais de résoudre l'équation sur un temps plus long.

    Comment prendre en compte le pas variable?

    L'idée est la suivante : après avoir obtenus om1 et om2 avec un pas h, on les compares avec om1_1 et om2_1 qui on été calculés avec un pas h/2 et qui sont donc plus précis.
    Si la difference est inférieur à la précision demandée, on garde le pas, sinon on le reduit et on recommence.
    Voir
    théorie :http://media4.obspm.fr/public/m2r/co...on3_2_4_6.html
    ou théorie et code: http://www.physique.usherbrooke.ca/p...HQ405_ipad.pdf (p61)
    ou encore http://users-phys.au.dk/fedorov/Numeric/10/odes.pdf
    et enfin un chapitre de Numerical recipes ici: http://www.pit.physik.uni-tuebingen....mrec/c16-2.pdf.

    J'ai donc tenter de rajouter une boucle Runge Kutta plus précise dans la première, de comparer les resultats et ensuite de modifier le pas si necassaire.

    Voici le code :
    Code :
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    67
    68
    69
    70
    71
    72
    73
    74
    75
    76
    77
    78
    79
    80
    81
    82
    83
    84
    85
    86
    87
    88
    # -*- coding: utf-8 -*-
    from numpy import*
    from scipy import*
    import matplotlib.pylab as plt
     
    #Condition initiales
    l1,l2,m1,m2=1,1,1,1                         #Longueur et masses des pendules 
    t0,tf= 0,2                                  #Temps de l'experience
    th1_0,th2_0,om1_0,om2_0=0,3,0,0             #Angles et vitesses initiales
    g=10
     
    #Les deux equations du pendule double
    dom1=lambda t,om1,om2,th1,th2 : (-m2*l1*sin(th1-th2)*cos(th1-th2)*om1**2-m2*l2*sin(th1-th2)*om2**2-m1*g*sin(th1-th2)*cos(th2))/l1*(m1+m2*(sin(th1-th2))**2)
    dom2=lambda t,om1,om2,th1,th2 : ((m1+m2)*l1*sin(th1-th2)*om1**2+m2*l2*sin(th1-th2)*cos(th1-th2)*om2**2+(m1+m2)*g*sin(th1-th2)*cos(th1))/l2*(m1+m2*(sin(th1-th2))**2)
     
    hmin,h,hmax=0.01,0.1,0.5                                             #Pas initial
    t,th1,th2,om1,om2 = t0,th1_0,th2_0,om1_0,om2_0      #Initialisation des variables
    th1_1,th2_1,om1_1,om2_1=th1_0,th2_0,om1_0,om2_0                                                
    p=0.001                                             #Precision demandée
     
    #Tableaux des données
    T=array([t])                                  
    Th1=array([th1])
    Th2=array([th2])
    Om1=array([om1])
    Om2=array([om2])
     
    #Boucle Runge-Kutta 4
     
    while t<tf:                                  
     
        if t+h>tf: h=tf-t
     
        th1 += h*om1
        th2 += h*om2
     
        k1 = h*dom1 (t,om1,om2,th1,th2)
        j1 = h*dom2 (t,om1,om2,th1,th2)
        k2 = h*dom1 (t+h/2,om1+h*k1/2,om2+h*j1/2,th1,th2)
        j2 = h*dom2 (t+h/2,om1+h*k1/2,om2+h*j1/2,th1,th2)
        k3 = h*dom1 (t+h/2,om1+h*k2/2,om2+h*j2/2,th1,th2)
        j3 = h*dom2 (t+h/2,om1+h*k2/2,om2+h*j2/2,th1,th2)
        k4 = h*dom1 (t+h,om1+h*k3,om2+h*j2,th1,th2)
        j4 = h*dom2 (t+h,om1+h*k3,om2+h*j2,th1,th2)
     
        om1 += k1/6 + k2/3 + k3/3 + k4/6
        om2 += j1/6 + j2/3 + j3/3 + j4/6
     
        h=h/2                                  #Division du pas par 2  
     
        for k in range (2):                      #Double boucle Runge Kutta avec un pas divisé par 2
     
            th1_1 += h*om1_1                     #===> Meilleur precision 
            th2_1 += h*om2_1
     
            k1_1 = h*dom1 (t,om1_1,om2_1,th1_1,th2_1)
            j1_1 = h*dom2 (t,om1_1,om2_1,th1_1,th2_1)
            k2_1 = h*dom1 (t+h/2,om1_1+h*k1_1/2,om2_1+h*j1_1/2,th1_1,th2_1)
            j2_1 = h*dom2 (t+h/2,om1_1+h*k1_1/2,om2_1+h*j1_1/2,th1_1,th2_1)
            k3_1 = h*dom1 (t+h/2,om1_1+h*k2_1/2,om2_1+h*j2_1/2,th1_1,th2_1)
            j3_1 = h*dom2 (t+h/2,om1_1+h*k2_1/2,om2_1+h*j2_1/2,th1_1,th2_1)
            k4_1 = h*dom1 (t+h,om1_1+h*k3_1,om2_1+h*j2_1,th1_1,th2_1)
            j4_1 = h*dom2 (t+h,om1_1+h*k3_1,om2_1+h*j2_1,th1_1,th2_1)
     
            om1_1 += k1_1/6 + k2_1/3 + k3_1/3 + k4_1/6
            om2_1 += j1_1/6 + j2_1/3 + j3_1/3 + j4_1/6
     
        h=h*2                                        #On retablit le pas
        err=array([abs(om1-om1_1),abs(om2-om2_1)])   #Calcule de l'erreur
        err=max(err)
     
        if err<=p:                          #Si l'erreur est acceptable
            T=append(T,t)                #Stockage des données
            Th1=append(Th1,th1)
            Th2=append(Th2,th2)
            Om1=append(Om1,om1)
            Om2=append(Om2,om2)
                               #Augmentation du pas 
            t+=h                        #et on passe a l'etape suivante
     
        else:                               #Sinon diminution du pas 
                                   #Et on recommence la meme etape
     
            h = h *min( max( 0.9 * (p/err) **0.2,0.1),4.0) 
            if h > hmax:
                h = hmax
            if h < hmin:
                h=hmin
    Mais le code n'est pas bon, il s'arrête rapidement est les valeurs sont très loin de celle obtenus sans faire varier h.
    J'ai l'impression que le problème vient du fait que même avec un pas plus petit, l'erreur calculée ne diminue pas, elle augmente. Dès lors, on ne se trouve plus jamais dans le cas 'err<=p' et le programme ne s'arrête plus.
    Est ce que quelqu'un sait pourquoi l'erreur augment alors que le pas diminue ?
    Est-ce que vous pensez que le problème est purement mathématique ou qu'il vient du programme?

    Si vous avez des remarques ou des commentaires sur le code, je prend!

    Merci d'avance pour vos réponses.

  2. #2
    Membre actif

    Homme Profil pro Guillaume
    Développeur informatique
    Inscrit en
    janvier 2012
    Messages
    27
    Détails du profil
    Informations personnelles :
    Nom : Homme Guillaume
    Localisation : France

    Informations professionnelles :
    Activité : Développeur informatique
    Secteur : Aéronautique - Marine - Espace - Armement

    Informations forums :
    Inscription : janvier 2012
    Messages : 27
    Points : 176
    Points
    176

    Par défaut

    Penses à factoriser ton code tu fais un copier/coller de deux fois la même portion de code ce qui n'est pas bien ^^.
    Pour ton problème, si tu passes par une erreur acceptable, il faut mettre à jour om1 et om2 pour qu'ils prennent om1_1 et om2_1 comme valeur sinon cela ne sert à rien. De même pour Om1 et Om2 qui doivent prendre les valeurs om1_1 et om2_1.
    Après je ne connais pas la méthode, mais d'après le document pdf, le pas de temps est modifié aussi p66 Ligne 42 du document code: http://www.physique.usherbrooke.ca/p...HQ405_ipad.pdf

  3. #3
    Invité de passage
    Inscrit en
    mai 2008
    Messages
    6
    Détails du profil
    Informations forums :
    Inscription : mai 2008
    Messages : 6
    Points : 2
    Points
    2

    Par défaut

    Merci pour ta réponse.

    Concernant le copier coller:
    Les deux boucles sont presques les même mais elles prennent des valeurs différente à chaque pas, ce sont en fait deux boucle complètement différentes. Il est vrai que j'aurais du utiliser une seule routine dans une définition et l'appeler deux fois avec des valeurs différentes.

    Pour la mise à jour des valeurs:
    d'après ce que j'ai compris, la deuxième boucle ne sert que de vérification, et on utilise les valeurs de la première si l'erreur est acceptable. C'est pourquoi ce sont om1 et om2 qui sont ajoutés.


    Par ailleurs, j'ai demander à mon prof d'info pourquoi l'erreur augmente alors que le pas diminue qui m'a répondu: que lorsqu'on utilises une routine à pas variable il faut faire attention au fait que lorsque le pas devient très petit, l'erreur numérique tend vers zéro (en théorie), par contre l'ereur que l'on fait (due à la précision de représentation) lors des calculs augmente... il faut donc trouver un optimum.

    Un algorithme est proposé dans les numerical recipes dont on peux s'inspirer.

    http://apps.nrbook.com/c/index.html, (aller aux pages pages 714-722).

    Finalement mise en place d'une routine à pas variable est finalement un peu trop complexe pour mon niveau de programmation, initié depuis quelques semaines...

    Je n'aurais pas le temps de corriger le programme dans les prochains jours, le sujet est donc résolue.

+ Répondre à la discussion
Cette discussion est résolue.

Liens sociaux

Règles de messages

  • Vous ne pouvez pas créer de nouvelles discussions
  • Vous ne pouvez pas envoyer des réponses
  • Vous ne pouvez pas envoyer des pièces jointes
  • Vous ne pouvez pas modifier vos messages
  •