Navigation

Inscrivez-vous gratuitement
pour pouvoir participer, suivre les réponses en temps réel, voter pour les messages, poser vos propres questions et recevoir la newsletter

Discussion: Runge Kutta d'ordre 4

  1. #1
    Nouveau Candidat au Club
    Homme Profil pro
    Lycéen
    Inscrit en
    juin 2019
    Messages
    6
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Âge : 16
    Localisation : France, Bas Rhin (Alsace)

    Informations professionnelles :
    Activité : Lycéen

    Informations forums :
    Inscription : juin 2019
    Messages : 6
    Points : 1
    Points
    1

    Par défaut Runge Kutta d'ordre 4

    Bonsoir, j'étudie la trajectoire d'un lanceur. Le lancer est supposé vertical, le système est soumis aux forces de frottements de l'air, à son poids (pas à la force de Coriolis). J'ai ainsi établi une équation différentielle, et je veux appliquer la méthode Runge Kutta d'ordre 4 pour la résoudre.
    Voilà le code. J'ai un message d'erreur qui s'affiche et je ne vois pas le problème :
    Code : Sélectionner tout - Visualiser dans une fenêtre à part
    1
    2
    if z2<100000:
    ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()
    Code : Sélectionner tout - Visualiser dans une fenêtre à part
    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
    alt,vit=[0],[0]
    z2,y,x = 0,0,0.6371009*(10**(-19))
    w=2*np.pi/(23*3600+59*60+4) #rad/s#vitesse de rotation de la terre(465,1 m/s )
     
    for i in range (l-3):
        if z2<100000:
            k1 = -g + (pousse[i] - 0.5*Cx*maitre_couple(tps[i])*masse_volumique(z2)*v**2)/masse[i] + x*w**2
        else:
            k1 = -g + pousse[i]/masse[i] + x*w**2
        if z2< 100000:
            k2 = -g + ((pousse[i]+pousse[i+1])/2 - 0.5*Cx*(np.pi*(5.425/2)**2 + 2*np.pi*(3.15/2)**2)*masse_volumique(z2)*((y+(pas[i]*k1/2))**2))/(masse[i]+masse[i+1])/2 + x*w**2
        else:
            k2 = -g +(pousse[i]+pousse[i+1])/2 + x*w**2
        if z2< 100000:
            k3 = -g + ((pousse[i]+pousse[i+1])/2 -0.5*Cx*(np.pi*(5.425/2)**2 + 2*np.pi*(3.15/2)**2)*masse_volumique(z2)*(y+1/2*pas[i]*k2)**2)/(masse[i]+masse[i+1])/2 + x*w**2
        else:
            k3 = -g + (pousse[i]+pousse[i+1])/2 + x*w**2
        if z2<100000:
            k4 = -g + ((pousse[i+1]+(y+pas[i]*k3)*debitm[i+1])/masse[i+1]) + x*w**2
        else:
            k4 = -g + (poussee[i+1]/masse[i+1]) + x*w**2
        y+=pas[i]*(k1/6 + k2/3 + k3/3 + k4/6)
        z2+=y*pas[i]
        x+=y*pas[i]
        alt.append(z2)
        vit.append(y)
    #plt.plot(tps[0:683],alt, label='altirkf')
    plt.plot(tps[0:684], vit)
    plt.legend
    plt.show()
    Pouvez-vous m'aider svp

  2. #2
    Membre expérimenté

    Homme Profil pro
    Ingénieur calcul scientifique
    Inscrit en
    mars 2013
    Messages
    724
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Localisation : France, Alpes Maritimes (Provence Alpes Côte d'Azur)

    Informations professionnelles :
    Activité : Ingénieur calcul scientifique

    Informations forums :
    Inscription : mars 2013
    Messages : 724
    Points : 1 432
    Points
    1 432

    Par défaut

    Le message d'erreur est clair. Vous tenter de comparer un tableau (z2) à une valeur. Vous meme vous vous y prendriez comment pour faire cela ? ...

    A vue de nez il vous manque des [i].

    2 remarques :
    1) scipy contient déjà un schéma runge kutta 4 (scipy.integrate.ode)
    2) Si vous utilisez numpy, vous n'avez justement plus besoin de parcourir vos tableaux puisque vous pouvez faire des opérations sur des tableaux entier sans avoir systématiquement à parcourir un à un les éléments de chaque array. Vous devriez vous plongez dans un tuto numpy.

  3. #3
    Nouveau Candidat au Club
    Homme Profil pro
    Lycéen
    Inscrit en
    juin 2019
    Messages
    6
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Âge : 16
    Localisation : France, Bas Rhin (Alsace)

    Informations professionnelles :
    Activité : Lycéen

    Informations forums :
    Inscription : juin 2019
    Messages : 6
    Points : 1
    Points
    1

    Par défaut

    Citation Envoyé par lg_53 Voir le message
    Le message d'erreur est clair. Vous tenter de comparer un tableau (z2) à une valeur. Vous meme vous vous y prendriez comment pour faire cela ? ...

    A vue de nez il vous manque des [i].

    2 remarques :
    1) scipy contient déjà un schéma runge kutta 4 (scipy.integrate.ode)
    2) Si vous utilisez numpy, vous n'avez justement plus besoin de parcourir vos tableaux puisque vous pouvez faire des opérations sur des tableaux entier sans avoir systématiquement à parcourir un à un les éléments de chaque array. Vous devriez vous plongez dans un tuto numpy.
    Bonjour,
    je vais essayer les schémas proposés par python alors. Le soucis, c'est que je n'arrive pas à définir proprement ma fonction f à mettre en argument, comme je travaille avec des tableaux.
    Mais ici, z2 est un floattant (il vaut zero au départ), que j'ajoute ensuite à un tabelau.

    Merci pour votre aide.

  4. #4
    Membre expérimenté

    Homme Profil pro
    Ingénieur calcul scientifique
    Inscrit en
    mars 2013
    Messages
    724
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Localisation : France, Alpes Maritimes (Provence Alpes Côte d'Azur)

    Informations professionnelles :
    Activité : Ingénieur calcul scientifique

    Informations forums :
    Inscription : mars 2013
    Messages : 724
    Points : 1 432
    Points
    1 432

    Par défaut

    Sipython vous dit que ce n'est pas un floattant, alors ce n'en est pas un ...
    Faites des print juste avant, vous verrez bien ce qu'il en retourne !

  5. #5
    Nouveau Candidat au Club
    Homme Profil pro
    Lycéen
    Inscrit en
    juin 2019
    Messages
    6
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Âge : 16
    Localisation : France, Bas Rhin (Alsace)

    Informations professionnelles :
    Activité : Lycéen

    Informations forums :
    Inscription : juin 2019
    Messages : 6
    Points : 1
    Points
    1

    Par défaut

    D'accord, je vais regarder ça.

  6. #6
    Nouveau Candidat au Club
    Homme Profil pro
    Lycéen
    Inscrit en
    juin 2019
    Messages
    6
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Âge : 16
    Localisation : France, Bas Rhin (Alsace)

    Informations professionnelles :
    Activité : Lycéen

    Informations forums :
    Inscription : juin 2019
    Messages : 6
    Points : 1
    Points
    1

    Par défaut

    Z2 est un tableau pour python, pourtant je l'initialise à z2=0. Sauriez-vous pourquoi c'est un tableau ?

  7. #7
    Membre expérimenté

    Homme Profil pro
    Ingénieur calcul scientifique
    Inscrit en
    mars 2013
    Messages
    724
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Localisation : France, Alpes Maritimes (Provence Alpes Côte d'Azur)

    Informations professionnelles :
    Activité : Ingénieur calcul scientifique

    Informations forums :
    Inscription : mars 2013
    Messages : 724
    Points : 1 432
    Points
    1 432

    Par défaut

    Votre code ne se suffit pas à lui même c-à-d que je ne peux l'éxécuter chez en l'état pour reproduire votre problème. Je n'ai pas les fonctions que vous appelez notamment. Dans ces conditions difficile de vous en dire plus.

    Faites des print de z2 un peu partout, réduisez l'étau pour savoir où est ce qu'il passe exactement d'un nombre à un tableau.

  8. #8
    Nouveau Candidat au Club
    Homme Profil pro
    Lycéen
    Inscrit en
    juin 2019
    Messages
    6
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Âge : 16
    Localisation : France, Bas Rhin (Alsace)

    Informations professionnelles :
    Activité : Lycéen

    Informations forums :
    Inscription : juin 2019
    Messages : 6
    Points : 1
    Points
    1

    Par défaut

    Je n'arrive pas à resserrer l'étau... Voilà mon code :
    Code : Sélectionner tout - Visualiser dans une fenêtre à part
    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
    g=9.81
    R=8.31 #J/K/mol
    gamma = 1.40  #gaz parfait diatomique, coeff Laplace
    m0=768789.3 #masse initiale en kg
    r=5.425 #diamètre coiffe en m
    Cx = 0.6 #supposé constant car vitesses élevées atteintes rapidement (graphe de reynolds, mach...)
     
    a=6.5*0.01 #en K/m
    T0=(15+273.15) #en kelvin
    p0= 101325 #en pascal
    M=28.966*0.01 #en kg/mol
     
    def temperature(z):
        return -a*z + T0
     
    def pression(z): #hypothèse d'un gradient de température constant
        return p0*((1-z*a/T0)**(M*g/(R*a)))
     
     
    def masse_volumique(z):
        return M*pression(z))/(R*temperature(z))
     
    pas=np.zeros(l)
    for i in range (l-2):
        pas[i]=(tps[i+1]-tps[i])
     
    #calcul du débit massique
    for i in range(l-2):
        debitm.append(abs((masse[i]-masse[i+1])/(tps[i+1]-tps[i])))

  9. #9
    Membre expérimenté

    Homme Profil pro
    Ingénieur calcul scientifique
    Inscrit en
    mars 2013
    Messages
    724
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Localisation : France, Alpes Maritimes (Provence Alpes Côte d'Azur)

    Informations professionnelles :
    Activité : Ingénieur calcul scientifique

    Informations forums :
    Inscription : mars 2013
    Messages : 724
    Points : 1 432
    Points
    1 432

    Par défaut

    Testez ce que vous présentez avant de poster. Je ne suis pas plus avancer il en manque encore.

    La fonction masse_volumique ne compile meme pas puisque contient une parenthèse de trop. La variable l, ainsi que la variable tps sont inconnues ... Je ne connais toujours pas non plus la valeur de pousse qui était dans le code de tout départ. Il me faut un code qui si je l'éxécute reproduis exactement le problème que vous dites (et pas un autre problème parce que je n'ai pas ce qu'il faut pour faire tourner votre code). Aller il faut être rigoureu, il faut se farcir l'exercice qui est certes fastidieux lorsqu'on débute mais indispensable.

    Prendre son temps de bien construire et présenter sa question est loin d'être une perte de temps. Dans 3 cas sur 4 ca vous permettera vous meme de trouver le bug tout seul. Et dans le reste des cas, ca vous fera un bon post et vous pourrez être dépanné rapidemment. Pas comme là ...

  10. #10
    Nouveau Candidat au Club
    Homme Profil pro
    Lycéen
    Inscrit en
    juin 2019
    Messages
    6
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Âge : 16
    Localisation : France, Bas Rhin (Alsace)

    Informations professionnelles :
    Activité : Lycéen

    Informations forums :
    Inscription : juin 2019
    Messages : 6
    Points : 1
    Points
    1

    Par défaut

    J'ai printé z2 dès la première boucle, et il s'agit déjà d'un tableau...
    Mon objectif est de tracer à partir de ces données, la vitesse puis la position de mon système en fonction du temps. Je dispose d'un fichier excel (ci-joint) avec des données (temps, poussée, masse). J'ai mis ces données dans des tableaux (respectivement appelés tps, pousse, masse). Je créer le tableau débit massique, debitm, en calculant une variation de masse par unité de temps. J'appelle l = len(tps).
    Comme les mesures de temps ne sont font pas à période régulière (toutes les sec puis toutes les 0.02sec, puis encore toutes les sec...), j'introduis un pas variable pour la méthode rk4, qui est la différence entre deux temps mesurés.
    Je distingue les cas z2<100000 et z2>100000, car je fais l'hypothèse que dans le second cas, la force de frottements est nulle (l'air n'est plus assez dense).

    Code : Sélectionner tout - Visualiser dans une fenêtre à part
    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
    89
    90
    91
    92
    93
    94
    95
    96
    97
     
    import numpy as np
    import scipy as sp
    import matplotlib.pyplot as plt
    from math import *
    import xlrd
    import scipy.integrate as integr
     
     
    g=9.81
    R=8.31 #J/K/mol
    gamma = 1.40  #gaz parfait diatomique, coeff Laplace
    m0=768789.3 #masse initiale en kg
    r=5.425 #diamètre coiffe en m
    Cx = 0.6 #supposé constant car vitesses élevées atteintes rapidement (graphe de reynolds, mach...)
     
     
    workbook = xlrd.open_workbook("mettre le nom du fichier")
    worksheet = workbook.sheet_by_index(0)
    tps = worksheet.col_values(0)
    pousse = worksheet.col_values(1)
    masse = worksheet.col_values(2)
    debit = worksheet.col_values(3)
    l=len(tps)
     
     
    tps, pousse, masse, debit = [], [],[], []
    for i in range (l-1):
        tps.append(float(worksheet.col_values(0)[i]))
        pousse.append(float(worksheet.col_values(1)[i]))
        masse.append(float(worksheet.col_values(2)[i]))
        debit.append(float(worksheet.col_values(3)[i]))
     
    tps1=np.zeros(l-2)
    for i in range (l-2) :
        tps1[i]=tps[i]
     
    debitm, pas, v_ejec=[],[],[]
     
    #calcul du débit massique
    for i in range(l-2):
        debitm.append(abs((masse[i]-masse[i+1])/(tps[i+1]-tps[i])))
     
    #pas variable pour equa diff
    pas=np.zeros(l)
    for i in range (l-2):
        pas[i]=(tps[i+1]-tps[i])
     
     
    #température en deg, pression en pascal et masse volumique de l'air en fonction de l'altitude z
    a=6.5*0.01 #en K/m
    T0=(15+273.15) #en kelvin
    p0= 101325 #en pascal
    M=28.966*0.01 #en kg/mol
     
    def temperature(z):
        return -a*z + T0
     
    def pression(z): #hypothèse d'un gradient de température constant
        return p0*((1-z*a/T0)**(M*g/(R*a)))
     
     
    def masse_volumique(z):
        return M*pression(z)/(R*temperature(z))
     
     
    def vitesse_son(z):
        return sqrt(gamma*R*temperature(z)/masse_volumique(z))
     
    #rk4fusee eau.xlsx
    alt,vit=[0],[0]
    z2=0
    y,x =0,0.6371009*(10**(-19))
    w=2*np.pi/(23*3600+59*60+4) #rad/s#vitesse de rotation de la terre(465,1 m/s )
     
    for i in range (l-3):
        if z2<100000:
            k1 = -g + (pousse[i] - 0.5*Cx*maitre_couple(tps[i])*masse_volumique(z2)*v**2)/masse[i] + x*w**2
        else:
            k1 = -g + pousse[i]/masse[i] + x*w**2
        if z2< 100000:
            k2 = -g + ((pousse[i]+pousse[i+1])/2 - 0.5*Cx*(np.pi*(5.425/2)**2 + 2*np.pi*(3.15/2)**2)*masse_volumique(z2)*((y+(pas[i]*k1/2))**2))/(masse[i]+masse[i+1])/2 + x*w**2
        else:
            k2 = -g +(pousse[i]+pousse[i+1])/2 + x*w**2
        if z2< 100000:
            k3 = -g + ((pousse[i]+pousse[i+1])/2 -0.5*Cx*(np.pi*(5.425/2)**2 + 2*np.pi*(3.15/2)**2)*masse_volumique(z2)*(y+1/2*pas[i]*k2)**2)/(masse[i]+masse[i+1])/2 + x*w**2
        else:
            k3 = -g + (pousse[i]+pousse[i+1])/2 + x*w**2
        if z2<100000:
            k4 = -g + ((pousse[i+1]+(y+pas[i]*k3)*debitm[i+1])/masse[i+1]) + x*w**2
        else:
            k4 = -g + (poussee[i+1]/masse[i+1]) + x*w**2
        y+=pas[i]*(k1/6 + k2/3 + k3/3 + k4/6)
        z2+=y*pas[i]
        x+=y*pas[i]
        alt.append(z2)
        vit.append(y)

  11. #11
    Membre expérimenté

    Homme Profil pro
    Ingénieur calcul scientifique
    Inscrit en
    mars 2013
    Messages
    724
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Localisation : France, Alpes Maritimes (Provence Alpes Côte d'Azur)

    Informations professionnelles :
    Activité : Ingénieur calcul scientifique

    Informations forums :
    Inscription : mars 2013
    Messages : 724
    Points : 1 432
    Points
    1 432

    Par défaut

    1) on n'a pas le fichier. Et en même temps ca complique la vie le fichier pour que les autres puissent tester. Faites nous donc une iitialisation de vos variables directement dans le code avec quelques valeurs. Par exemple :

    Code : Sélectionner tout - Visualiser dans une fenêtre à part
    tps=np.array([1,2,3,3.02,3.04,3.06,3.08,3.1,4.1,5.1])
    et de même pour pousse,masse et debit. Ca évite de vous réclamer un fichier qu'on n'a pas et ca évite égalament d'avoir des données de tailles conséquentes (là on réduit).

    2) Vous n'avez pas bien tester. Pour tester correcteemnt il faut vous remettre dans un environnement propre (effacer toute les variables globales notamment). Si vous ne savez pas faire cela, redémarrer votre IDE, mais il contiennent souvent des commande "clear" ou "reset". Là si vous le faites vous allez voir que des choses manquent encore.

Discussions similaires

  1. Runge Kutta d'ordre 4
    Par Physicien dans le forum Algorithmes et structures de données
    Réponses: 1
    Dernier message: 08/12/2006, 16h50
  2. probleme de divergence avec runge kutta d'ordre 2 pour un pendule simple
    Par fab13 dans le forum Algorithmes et structures de données
    Réponses: 1
    Dernier message: 25/11/2006, 20h19
  3. Runge-Kutta à une variable?
    Par PadawanDuDelphi dans le forum Algorithmes et structures de données
    Réponses: 18
    Dernier message: 22/11/2006, 09h20
  4. probleme avec runge kutta dimension 4
    Par fab13 dans le forum Algorithmes et structures de données
    Réponses: 3
    Dernier message: 14/11/2006, 21h47

Partager

Partager
  • Envoyer la discussion sur Viadeo
  • Envoyer la discussion sur Twitter
  • Envoyer la discussion sur Google
  • Envoyer la discussion sur Facebook
  • Envoyer la discussion sur Digg
  • Envoyer la discussion sur Delicious
  • Envoyer la discussion sur MySpace
  • Envoyer la discussion sur Yahoo