Runge Kutta d'ordre 4

Version imprimable

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:

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:

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

10/06/2019, 09h23
lg_53

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.
10/06/2019, 09h28
Zaza682904

Citation:

Envoyé par lg_53

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.
10/06/2019, 09h36
lg_53

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 !
10/06/2019, 09h42
Zaza682904

D'accord, je vais regarder ça.
10/06/2019, 09h49
Zaza682904

Z2 est un tableau pour python, pourtant je l'initialise à z2=0. Sauriez-vous pourquoi c'est un tableau ?
10/06/2019, 09h57
lg_53

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.

Je n'arrive pas à resserrer l'étau... Voilà mon 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
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])))

10/06/2019, 11h24
lg_53

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à ...

1 pièce(s) jointe(s)

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:

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))
 
#rk4Pièce jointe 482256
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)

10/06/2019, 14h30
lg_53

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:

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.