Méthode Euler Explicite

Version imprimable

2 pièce(s) jointe(s)

Bonjour tout le monde;

Voilà, j'ai un petit projet à faire avec l'intégration numérique par plusieurs méthodes.
Je vois bien le principe je dirais analytique (crayon et papier), mais quand à programmer cela j'ai bien plus de mal.

Voilà l'équation, la méthode sur papier que j'ai trouvé et un début de code dont j'aurais besoin de la plus grande aide ^^

L'équation correspond à la première des 3.

Merci d'avance !
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 # usr/bin/python from pylab import * import numpy as np ####################################### # ORBITES AUTOUR D'UN TROU NOIR DE SCHWARZSCHILD # ####################################### # -------------------------------- # # Methode d'Euler explicite # # -------------------------------- # Rs = 2 # Rayon de Swharzschild m = 1 # Masse quelconque d'une particule c = 1 # Vitesse lumiere r = 0.8 # rayon r E = 1 # Energie particule F = ( (1-(Rs/r))**(-1) ) * E/(m*c**2) def Euler(F, y0, t):

09/11/2015, 20h36
marco056

Il y a pléthore de solutions grâce à moteur de recherche basique.
09/11/2015, 20h58
Andromedae

J'ai déjà tapé Euler sur ce forum, 5 sujets, aucuns ne m'aident réellement :/
09/11/2015, 22h25
marco056

Citation:

Envoyé par Andromedae

J'ai déjà tapé Euler sur ce forum, 5 sujets, aucuns ne m'aident réellement :/

Pas très perspicace :
https://www.google.fr/search?client=...D9Ss8wft_4SQBA
Il suffit de taper "méthode euler python" et il tombe des dizaines de réponses pertinentes (les solutions sont même données !).
10/11/2015, 09h21
Andromedae

Tu crois que je ne l'ai pas fait ?
J'attendais juste d'avoir quelques pistes en rapport avec mes équations ... qui ne sont pas y' = y ...
10/11/2015, 10h15
marco056

C'est sûr : si tu attends que quelqu'un fasse tes devoirs, tu risques d'attendre longtemps.
Propose un bout de code, montre un peu de bonne volonté et tu auras sans doute de l'aide ici.

Je comprends la discorde, mon code ne s'était pas collé en entier ..
Le code est général, je n'ai pas encore défini y0 et t0 car j'ai un peu de mal à les trouver sachant que ma fonction ne dépend pas de tau :

y'(tau) = constante

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
# usr/bin/python
 
 
from pylab import *
import numpy as np
 
 
 
##################################################
 
# ORBITES AUTOUR D'UN TROU NOIR DE SCHWARZSCHILD #
 
##################################################
 
 
 
 
# ------------------------- #
# Methode d'Euler explicite #
# ------------------------- #
 
Rs = 2                        # Rayon de Swharzschild
m = 1                        # Masse quelconque d'une particule
c = 1                        # Vitesse lumiere
r = 0.8                        # Rayon r
E = 1                        # Energie particule
 
 
F = ( (1-(Rs/r))**(-1) ) * E/(m*c**2)        # Fonction (3a) du sujet
 
def Euler (F, t0, y0, tf, N) :
 
     h = (tf - t0)/N                # Definition du pas fixe
 
     Ly, Lt = [y0],[t0]            # Stockage du point initial
 
     for k in range (N) :
 
         y0 = y0 + h*np.array(F(t0,y0))
         t0+ = h
         Ly.append(y0)
         Lt.append(t0)
 
     return np.array(Lt), np.array(Ly)

10/11/2015, 12h52
marco056

Je comprends mieux aussi :mrgreen:
10/11/2015, 13h00
Andromedae

Je pense que :

y(0) = (1- (Rs/r))^-1 * E/mc² avec r = r^n donc n=1

y(tau) = (1- (Rs/r^n))^-1 * E/mc²
10/11/2015, 19h47
marco056

Je ne comprends pas trop : ton pb est résolu ou non (il est marqué comme oui) ?
11/11/2015, 06h18
tyrtamos

Bonjour,

Sans discuter de la validité de l'algorithme, il y a des choses bizarres dans le code Python donné:

- Si le F transmis comme argument dans la fonction Euler est le même que le F du début, il y a un problème: F est écrit comme une variable, mais présenté et utilisé comme une fonction: F(t0,y0)

- je ne comprends pas l'intérêt du np.array de la ligne y0 = y0 + h*np.array(F(t0,y0))

- dans la ligne "t0+ = h", il ne doit pas y avoir d'espace entre le "+" et le "=" (erreur de syntaxe)

Erreur de manipulation : sujet non résolu

Citation:

Si le F transmis comme argument dans la fonction Euler est le même que le F du début, il y a un problème: F est écrit comme une variable, mais présenté et utilisé comme une fonction: F(t0,y0)

F correspond à mon équation différentielle. Il faudrait donc que je supprime de l'argument Euler alors ?

Citation:

je ne comprends pas l'intérêt du np.array de la ligne y0 = y0 + h*np.array(F(t0,y0))

Je pensais lister les valeurs dans un tableau pour après en faire un graphique (avec matplotlib)

Citation:

dans la ligne "t0+ = h", il ne doit pas y avoir d'espace entre le "+" et le "=" (erreur de syntaxe)

Autant pour moi, erreur de syntaxe en effet

Le code donnerait cela alors avec déjà ces quelques modifications :
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 # usr/bin/python from pylab import * import numpy as np ################################################## # ORBITES AUTOUR D'UN TROU NOIR DE SCHWARZSCHILD # ################################################## # ------------------------- # # Methode d'Euler explicite # # ------------------------- # Rs = 2 # Rayon de Swharzschild m = 1 # Masse quelconque d'une particule c = 1 # Vitesse lumiere r = 0.8 # Rayon r E = 1 # Energie particule F = ( (1-(Rs/r))**(-1) ) * E/(m*c**2) # Fonction (3a) du sujet def Euler (t0, y0, tf, N) : h = (tf - t0)/N # Definition du pas fixe Ly, Lt = [y0],[t0] # Stockage du point initial for k in range (N) : y0 = y0 + h*np.array(F(t0,y0)) t0 += h Ly.append(y0) Lt.append(t0) return np.array(Lt), np.array(Ly)

11/11/2015, 17h53
marco056

J'ai un gros doute : r ne varie pas ?
11/11/2015, 17h55
Andromedae

Ce que je pense, c'est que pour y0 ==> r^n = r donc n = 1
Mais à chaque itération, n varie donc r varie

Mais je ne vois pas comment introduire cette fonction r(tau)^n

La méthode d'Euler devant donner je pense : Fonction (r^(n+h)) - Fonction(r^(n)) / h
11/11/2015, 19h22
marco056

Le souci, à mon avis, est en physique, avant d'en effectuer la programmation.
Que cherchez-vous ? A déterminer r ?

Bonjour,

Je reviens vers vous car j'ai continué à chercher, et peut-être à légèrement progresser.

Voici le code amélioré par rapport à la dernière fois :

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
#!usr/bin/python 
 
from pylab import *
from scipy import *
from scipy .integrate import odeint		# Pour resoudre les equations differentielles
import numpy as np
import math as m
import matplotlib .pyplot as plt		# Pour tracer des graphiques
 
 
										##################################################
 
										# ORBITES AUTOUR D'UN TROU NOIR DE SCHWARZSCHILD #
 
										##################################################
 
 
 
# Definition des variables :
#---------------------------
 
Rs = 1.0						# Rayon de Swharzschild
r = 20.0*Rs						# Rayon r
E = 10.0						# Energie particule = mc**2
L = 0.0							# Moment cinetique de la particule
 
# Implicitement on a :
m = 1.0							# Masse quelconque d'une particule
c = 1.0							# Vitesse lumiere
 
 
V0 = 0.0
 
 
# Definition du systeme d'equations :
#------------------------------------
 
t_point = 1/(1-(Rs/r)) * E						# Equation differentielle (3a) du sujet / t_point = dt/dtau
 
phi_point = L/(m*r**2)							# Equation differentielle (3b) du sujet / phi_point = dphi/dtau
 
V0**2 == E**2 - (1-(Rs/r))*(1+(L**2/r**2))				# Equation differentielle (3c) du sujet
 
A0 = (Rs/(2*r**2)) + (3/2)*((Rs*L**2 - r*L**2)/(r**4))			# Equation differentielle (3c) du sujet
 
 
print ' '
print 't_point vaut :' ,t_point
print 'phi_point vaut :' ,phi_point
print 'V0**2 vaut :' ,V0**2
print 'A0 vaut :' ,A0
print ' '
 
 
# ----------------------------------------------------------------- # 
# Methode d'Euler explicite pour le mouvement de chute libe radiale #
# ----------------------------------------------------------------- # 
 
t0 = 0.0
tf = 100.0
y0 = 0.0
n = 10000.0
 
def Euler (t_point, phi_point, V0, A0, t0, y0, tf, n) :
 
	h = (tf - t0)/n				# Definition du pas
 
	t = np.arange(t0,tf,h)
 
	y = np.zeros(n)
 
	y[0] = y0
 
	for i in range (n-1) :
 
		x[i+1] = x[i] + h
		y[i+1] = y[i] + h*t_point(y[i],x[i])
 
		print x[i+1]	
	return(x,y)

==> Mes fonctions sont bien définies pour le problème que je pose.
==> Dans mon programme :

* t_point <=> dt/dtau
* phi_point <=> dphi/dtau
* Vo² <=> (dr/dtau)²
* Ao <=> d/dtau (Vo²)

J'ai donc un système de 4 équations du premier ordre

==> Par ailleurs, dans mon programme, m = c = 1 Ce qui permet d'alléger mes équations.

==> Le code fonctionne jusqu'au print de mes fonctions, mais j'ai de gros soucis pour mettre en place ma méthode d'Euler explicite.
==> La quantité à faire varier étant "tau", je bloque sur Euler.

Quelqu'un pourrait-il m'aiguiller un peu ?

Merci d'avance !

Citation:
Envoyé par Andromedae
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 #!usr/bin/python from pylab import * from scipy import * from scipy .integrate import odeint # Pour resoudre les equations differentielles import numpy as np import math as m import matplotlib .pyplot as plt # Pour tracer des graphiques ################################################## # ORBITES AUTOUR D'UN TROU NOIR DE SCHWARZSCHILD # ################################################## # Definition des variables : #--------------------------- Rs = 1.0 # Rayon de Swharzschild r = 20.0*Rs # Rayon r E = 10.0 # Energie particule = mc**2 L = 0.0 # Moment cinetique de la particule # Implicitement on a : m = 1.0 # Masse quelconque d'une particule c = 1.0 # Vitesse lumiere V0 = 0.0 # Definition du systeme d'equations : #------------------------------------ t_point = 1/(1-(Rs/r)) * E # Equation differentielle (3a) du sujet / t_point = dt/dtau phi_point = L/(m*r**2) # Equation differentielle (3b) du sujet / phi_point = dphi/dtau V0**2 == E**2 - (1-(Rs/r))*(1+(L**2/r**2)) # Equation differentielle (3c) du sujet A0 = (Rs/(2*r**2)) + (3/2)*((Rs*L**2 - r*L**2)/(r**4)) # Equation differentielle (3c) du sujet print ' ' print 't_point vaut :' ,t_point print 'phi_point vaut :' ,phi_point print 'V0**2 vaut :' ,V0**2 print 'A0 vaut :' ,A0 print ' ' # ----------------------------------------------------------------- # # Methode d'Euler explicite pour le mouvement de chute libe radiale # # ----------------------------------------------------------------- # t0 = 0.0 tf = 100.0 y0 = 0.0 n = 10000.0 def Euler (t_point, phi_point, V0, A0, t0, y0, tf, n) : h = (tf - t0)/n # Definition du pas t = np.arange(t0,tf,h) y = np.zeros(n) y[0] = y0 for i in range (n-1) : x[i+1] = x[i] + h y[i+1] = y[i] + h*t_point(y[i],x[i]) print x[i+1] return(x,y)
==> Mes fonctions sont bien définies pour le problème que je pose.
==> Dans mon programme :

* t_point <=> dt/dtau
* phi_point <=> dphi/dtau
* Vo² <=> (dr/dtau)²
* Ao <=> d/dtau (Vo²)

J'ai donc un système de 4 équations du premier ordre

==> Par ailleurs, dans mon programme, m = c = 1 Ce qui permet d'alléger mes équations.

==> Le code fonctionne jusqu'au print de mes fonctions, mais j'ai de gros soucis pour mettre en place ma méthode d'Euler explicite.
==> La quantité à faire varier étant "tau", je bloque sur Euler.

Quelqu'un pourrait-il m'aiguiller un peu ?

Merci d'avance !
Bonjour,
Je ne comprends pas grand chose, mais:

Code:

V0**2 == E**2 - (1-(Rs/r))*(1+(L**2/r**2))

Signifie tester l'égalité entre "0.0**2" et le reste. C'est le but recherché?
Cela met une erreur en mettant "=" comme pour les autres? Alors on essaye des trucs…
D'autre part:
Code:

1 2 t_point = 1/(1-(Rs/r)) * E phi_point = L/(m*r**2)…
ne sont pas des équations différentielles, mais des constantes.
Code:

1 2 x[i+1] = x[i] + h y[i+1] = y[i] + h*t_point(y[i],x[i])
D'où vient la liste "x"?
Ainsi que dit plus haut, "t_point" est une constante, pas une liste.

L'usage du shell Python est fortement conseillé!

Aiguillage: avant d'envisager pylab, scipy, numpy, matplotlib… il serait approprié de commencer par un tutoriel Python, tel que celui de G.Swinnen. Gratuit.

Clodion

PS:
-il n'y a qu'une fonction de définie et elle n'est pas appelée.
-non, il n'y a pas quatre équations de premier ordre. Il y a quatre constantes.
-le code fonctionne (vraiment?) pour les imports et les déclarations de constantes (et encore!).