Modèle de Newton

1) Introduction

Les lois physiques fondamentales prennent vie à travers les modèles et leurs résolutions numériques. La programmation de modèles et la simulation que nous préférons appeler expérimentations exactes, constituent la deuxième face de la théorie, son aspect pratique en quelque sorte.

Le modèle est un système d'équations différentiels que l'on résoud numériquement de façon approchée par la méthode d'Euler. Cela met en oeuvre un calcul itératif qui dans la plus part des cas n'est qu'approché mais dont la précision peut être aussi grande que l'on veut. Notez que rien n'interdit qu'un calcul itératif puisse être exacte et constituer la loi proprement dite.

On propose aux lecteurs une théorie physique accompagnées d'un modèle qu'il pourra à loisir expérimenter sur son ordinateur personnel avec différentes conditions initiales.

On appelle expérimentation exacte, le calcul du modèle par l'ordinateur, en essayant un certain nombre de conditions initiales et de paramètres. Tandis qu'on appelle expérimentation réel, la comparaison du modèle avec la réalité.

Par la seul expérimentation exacte, on peut tirer certains enseignements. Et cette science basée sur l'expérimentation du calcul de modèle a l'avantage d'être abordable par tous, et d'être transmissible par simple duplication. Elle peut donc constituer à elle seul un véhicule pour transmettre un argumentaire, autrement dit un vecteur pour la propagande.

La théorie seule ne suffit pas pour convaincre, car l'utilisateur final n'est pas en mesure de vérifier ou seulement même d'en mesurer la pertinence.

Le modèle permet à l'utilisateur final d'expérimenter. Il permet d'observer la cohérence de la théorie en vérifiant que dans les simulations, les propriétés attendues se réalisent bien.

L'algorithme et la loi sont misent sur un même pied d'égalité. Ils sont transmis par copie, la copie d'une théorie et la copie d'un programme informatique.

Cette dualité offre un mécanisme de vérification d'erreur de copie de haut-niveau, l'exécution du programme vérifiant les propriétés intrinsèques de la théorie.

La charge de cette vérification revient à l'utilisateur final. Et selon la maxime de saint Thomas "Je ne crois que ce que je vois", la preuve n'est plus antérieur au programme, inévitablement soumise à la tentation d'appropriation et de rétention d'information de la part du concepteur, mais postérieur au programme, s'effectuant lors de l'exécution de celui-ci, sous l'entière responsabilité de l'utilisateur final qui satisfait ainsi à la maxime.

L'expérimentation exacte se comporte alors comme un élément du démonstrateur de théorèmes, apportant au coup par coup les indices de vérité.

La théorie que nous exposons est la théorie de la gravitation de Newton. On l'accompagne d'un modèle, d'un algorithme, ou autrement dit, d'un programme informatique. Ce programme calcul numériquement de façon approchée la trajectoire de `N` particules en interaction gravitationnelle.

2) Cinématique

Une particule est une masse ponctuelle possédant une trajectoire. La position de la particule se note `vecx`. Les différentes dérivés par rapport au temps de `vecx` définissent la vitesse `vecv`, l'accélération `veca`, la pente `vecp`, etc..

Position :
 
`vecx`
Vitesse :
 
`vecv`
`=`
`(dvecx)/(dt)`
Acceleration :
 
`veca`
`=`
`(dvecv)/(dt)`
`=`
`(d^2vecx)/(dt^2)`
Pente :
 
`vecp`
`=`
`(dveca)/(dt)`
`=`
`(d^2vecv)/(dt^2)`
`=`
`(d^3vecx)/(dt^3)`

A l'instant `0`, la particule se situe à la position `vec(x_0)` avec une vitesse `vec(v_0)`, une accélération `vec(a_0)`, une pente `vec(p_0)`, etc.. L'indice ici désignent un instant particulier.

A l'instant `t`, la particule se situe à une position `vecx` avec une vitesse `vecv`, une accéleration `veca`, une pente `vecp`, etc..

La position d'une particule que l'on note parfois `vecx(t)` pour préciser qu'elle évolue en fonction du temps `t`, se développe selon la série de Taylor :

`vecx(t) = vec(x_0) + vec(v_0)t + vec(a_0)t^2/(2!) + vec(p_0)t^3/(3!) + ...`

De même pour la vitesse et l'accélération :

`vecv(t) = vec(v_0) + vec(a_0)t + vec(p_0)t^2/(2!) + ...`

`veca(t) = vec(a_0) + vec(p_0)t + ...`

3) Dynamique

La légende dit que Newton décrit dans les pages collées de son oeuvre maîtresse Philosophiae naturalis principia mathematica, une explication symbolique des bases de la gravité, des révélations intuitives cachées pour échapper au foudre de l'inquisition religieuse toujours très présente à cette époque. La matière aime la matière. La gravité serait l'expression de l'Amour universel tel qu'on peut la conçevoir entre deux êtres. Sa force serait donc une force d'attraction, et qui serait proportionnelle à la charge d'amour de chacun des deux corps qui en constitue ainsi leur masses. Puis tout le monde sait que l'attraction amoureuse est inversement proportionnelle au carré de la distance. Et c'est Gauss qui expliquera plus-tard de façon plus rationnelle cette décroissance de la force en fonction de la distance par le principe de conservation du flux du champ.

Force de gravité selon le modèle newtonien :

Force de gravité :

`f = G(Mm)/(r^2)`


`f` : Force de gravité
`G` : Constante des forces de gravité
`M` : Masse de la première particule
`m` : Masse de la seconde particule
`r` : Distance entre les deux particules

Les forces sont attractives, centrales, égales et opposées. Et elles s'ajoutent.

On en déduit une formulation plus complète :

Force de gravité causée par la particule source `M`
et subit par la particule réceptrice `m` :

`vec(f) = - vec(u)G(Mm)/(r^2)`

`vec(r) = vec(x)-vec(X)`

`r=|vecr|`

`vec(u) = vec(r)/r`

`vecf` : Force de gravité agissant sur la particule réceptrice `m`.
`G`
: Constante des forces de gravité.
`m`
: Masse de la particule réceptrice.
`vec(x)`
: Position de la particule réceptrice.
`M`
: Masse de la particule source du champ.
`vec(X)`
: Position de la particule source du champ.
`vec(r)`: Vecteur distance de la particule source à la particule réceptrice.
`vec(u)`: Vecteur unitaire de la particule source vers la particule réceptrice.
`r` : Distance entre les deux particules.
Force de gravité exerçée par `N` particules :

`AA(i,j) in {1..N}^2`

`vec(f_j) = - Gm_jsum_(i!=j) (vec(u_(ij))m_i)/(r_(ij)^2)`

`vec(r_(ij)) = vec(x_j)-vec(x_i)`

`r_(ij) = |vec(r_(ij))|`

`vec(u_(ij)) = vec(r_(ij))/r_(ij)`

`vec(f_i)` : Force de gravité s'appliquant à la particule n°`i`.
`G` : Constante des forces de gravité.
`m_i` : Masse de la particule n°`i`.
`vec(x_i)` : Position de la particule n°`i`.
`vec(r_(ij))` : Vecteur distance de la particule n°`i` à la particule n°`j`.
`vec(u_(ij))`: Vecteur unitaire de la particule n°`i` vers la particule n°`j`.
`vec(x_i)` : Position de la particule n°`i`.
`r_(ij)`
: Distance entre la particule n°`i` et la particule n°`j`.

Les données sont la position et la vitesse pour les `N` particules :

`vec(x_1), vec(x_2),..., vec(x_N)`

`vec(v_1), vec(v_2),..., vec(v_N)`

Les indices ici désignent la particule concernée.

Le symbole `sum_(i!=j)` indique une somme pour toutes les valeur possibles de l'indice `i` mais distincte de `j`.

De même nous avons :

`sum_(i!=j)`
`=`
`sum_(i in ({1..N}-{j}))`
`sum_i`
`=`
`sum_(i in {1..N})`
`=`
`sum_(i=1)^N`
`sum_isum_(j<i)`
`=`
`sum_(i in {2..N)) sum_(j in {1..i-1})`
`=`
`sum_(i=2)^N sum_(j=1)^(i-1)`

Inertie :

Principe d'inertie :

`vecf = mveca`

`vecf` : Force s'exerçant sur la particule.
`m` : Masse de la particule.
`veca` : Accélération de la particule.

4) Energie

Energie cinétique d'une particule

Energie cinétique d'une particule :

`E_c = 1/2mv^2`

`E_c` : Energie cinétique de la particule.
`m` : Masse de la particule.
`v` : Vitesse de la particule.
Energie cinétique d'un système de `N` particules :

`E_c = 1/2sum_i m_iv_i^2`

`E_c` : Energie cinétique du système.
`m_i` : Masse de la particule n°`i`.
`v_i` : Vitesse de la particule n°`i`.

Potentiel de gravité

Potentiel de gravité à la position `vecx`
dans un système de `N` particules :

`U = G sum_i m_i/r_i`


`vec(r_i) = vec(x)-vec(x_i)`

`r_i=|vec(r_i)|`

`U`
: Potentiel de gravité à la position `vecx`.
`G`
: Constante des forces de gravité.
`m_i`
: Masse de la particule n°`i`.
`vec(r_i)`
: Vecteur distance de la particule n°`i` à la position `vecx`.
`r_i`
: Distance entre la particule n°`i` et la position `vecx`.
Energie potentiel de gravité d'un système de `N` particules :

`E_p = G sum_i sum_(j<i) m_i/r_(ij)`

`vec(r_(ij)) = vec(x_j)-vec(x_i)`

`r_ij=|vec(r_(ij))|`


`E_p` : Energie potentielle du système.
`G` : Constante des forces de gravité.
`m_i`
: Masse de la particule n°`i`.
`vec(r_(ij))` : Vecteur distance de la particule n°`i` à la particule n°`j`.
`r_(ij)` : Distance entre la particule n°`i` et la particule n°`j`.

Energie totale du système

Energie totale d'un système de `N` particules :

`E = E_c + E_p`
`E = (1/2sum_i m_iv_i^2) + G (sum_i sum_(j<i) m_i/r_(ij))`

`E_p` : Energie potentielle du système.
`E_c` : Energie cinétique du système.
`G` : Constante des forces de gravité.
`v_i` : Vitesse de la particule n°`i`.
`m_i`
: Masse de la particule n°`i`.
`r_(ij)`
: Distance entre la particule n°`i` et la particule n°`j`.

5) Système d'équations différentielles

Le modèle newtonien appliqué à `N` particules se ramène à un système d'équations différentiels du second ordre, l'accéleration étant la dérivé seconde de la position :

`AA(i,j) in {1..N}^2`
`(d^2vec(x_j))/dt^2 = - Gsum_(i!=j) m_i/(r_(ij)^2)vec(u_(ij))`
`vec(r_(ij)) = vec(x_j)-vec(x_i)`
`r_(ij) = |vec(r_(ij))|`
`vec(u_(ij)) = vec(r_(ij))/r_(ij)`

6) Méthode de résolution numérique approchée

On résoud le système d'équations différentielles numériquement de façon approchée, par la méthode d'Euler. Cela consiste à ségmenter le temps en une succession d'intervalles de duré `Deltat`. Pendant chaque intervalle de temps `Deltat` on considère que l'accélération de chaque particule évolue linéairement c'est à dire que la pente est constante.

A l'instant `0`, la particule se situe à la position `vec(x_0)` avec une vitesse `vec(v_0)`, une accélération `vec(a_0)` et une pente `vec(p_0)`. Puis à l'instant `Deltat`, la particule se situe à la position `vec(x_1)` avec une vitesse `vec(v_1)`, une accéleration `vec(a_1)` et une pente `vec(p_1)`. Le temps vaut `0` lorsque la particule est au point `vec(r_0)` et vaut `Deltat` lorsque la particule est au point `vec(x_1)`. La pente étant constante pendant toute la durée de l'intervalle, nous avons par intégration :

`vec(a_1) = vec(a_0) + vec(p_0) Deltat`
`vec(v_1) = vec(v_0) + vec(a_0)Deltat + vec(p_0)(Deltat^2)/2`
`vec(x_1) = vec(x_0) + vec(v_0)Deltat + vec(a_0)(Deltat^2)/2 + vec(p_0)(Deltat^3)/6`

Notez qu'ici l'indice désigne un instant, une date précise, un temps, un nombre entier de `Deltat`.

7) Estimation de l'erreur dû à l'approximation

Considérons la série de Taylor de la position d'une particule `vecx` appliquée à l'intervalle `Deltat` :

`vec(x_1) = vec(x_0) + vec(v_0)Deltat + vec(a_0)(Deltat^2)/2 + vec(p_0)(Deltat^3)/(3!) + (dvec(p_0))/dt(Deltat^4)/(4!) + ...`

Les indices ici désigne un instant correspondant à un nombre entier de `Deltat`.

Des considérations générales permettent de supposer que les fluctuations de la position, de la vitesse, de l'accélération, ainsi que de n'importe quelle dérivée `n`-ième de la position, sont d'un même ordre de grandeur que l'on peut majorer.

On note quelque fois la dérivée n-ième de `vecx(t)` comme suit :

`vecx^((n)) = (d^nvecx)/(dt^n)`

L'erreur `epsilon` commise en arrêtant la série de Taylor au quatrième terme, est du même ordre de grandeur que le cinquième terme qui est :

`epsilon  ≈  |(dvec(p_0))/(dt)|(Deltat^4)/24`

`epsilon  ≈  |vec(p_1)-vec(p_0)|/(Deltat) (Deltat^4)/24`

`epsilon  ≈  |vec(p_1)-vec(p_0)|(Deltat^3)/24`

Donc on choisie `Deltat` comme suit afin d'avoir un même ordre d'erreur tout au long du calcul de la trajectoire.

`Deltat = ( (24 epsilon )/|vec(p_1)-vec(p_0)|)^(1/3)`

Et dans le cas d'un système à `N` particules, on considère la valeur maximum `|vec(p_1)-vec(p_0)|` parmis toutes les particules :

`Deltat = ( (24 epsilon )/("Max"(|vec(p_1)-vec(p_0)|)))^(1/3)`

8) Algorithme

Pour chaque particule on connait : `vec(x_0)`, `vec(v_0)`, et on en déduit la force de gravité `vec(f_0)` et donc l'accélération `vec(a_0)`.

Puis on connait la pente constante pour l'intervalle de temps précédent, que l'on note `vec(p_(pre))`, et on va déterminer la pente constante pour l'intervalle de temps actuel, que l'on note `vec(p_0)`, à l'aide du calcul itératif suivant :

On initialise `vec(p_0) = vec(p_(pre))`

On calcul les positions approchées à l'instant `Deltat` suivant :

`vec(v_1) = vec(v_0) + vec(a_0)Deltat + vec(p_0)(Deltat^2)/2`

`vec(x_1) = vec(x_0) + vec(v_0)Deltat + vec(a_0)(Deltat^2)/2 + vec(p_0)(Deltat^3)/6`

On calcul la force de gravité `vec(f_1)`, l'accélération `vec(a_1)` et on recalcule de nouveau la pente `vec(p_0)`.

On calcul parmis toutes les particules, la variation de pente la plus grande : `"Max"(|vec(p_0)-vec(p_(pre))|)`

On recalcule l'intervalle de temps : `Deltat = ( (24 epsilon )/("Max"(|vec(p_0)-vec(p_(pre))|)))^(1/3)`

Puis on réitère le calcul jusqu'à ce que la variation du résultat soit inférieure à `epsilon`. Dans la pratique on répète l'itération 4 fois systématiquement.

De nombreux langages de programmation proposent des calculs avec une précision arbitrairement grande. On peut donc choisir une valeur d'`epsilon` arbitrairement petite et obtenir ainsi une précision arbitraire grande. Néanmoins la propagation des erreurs étant de nature exponentielle, on est vite épuisé par le temps de calcul.

9) Exemple

Trajectoire de 4 particules de même masse, lancées d'une façon symétrique (ce qui explique la symétrie centrale des trajectoires) :

Programme en go :

package main

import (
	"fmt"
"math" "github.com/veandco/go-sdl2/sdl"
) type Particule struct { // À l'instant présent m, // Masse x0, y0, z0, // Position vx0, vy0, vz0, // Vitesse ax0, ay0, az0, // Accélération px0, py0, pz0, // Pente // À l'instant dt suivant x1, y1, z1, // Position vx1, vy1, vz1, // Vitesse ax1, ay1, az1, // Accélération fx1, fy1, fz1, // Force s'excerçant sur la particule // Durant l'intervalle de temps dt précédent px_, py_, pz_ float64 // Pente } const N = 4 // Nombre de particules var P [N]Particule // Liste des N particules const epsilon = 1e-15 // Erreur tolérée à chaque itération var G = 100.0 // Constante des forces de gravité var dt, t, // Intervalle de temps et temps cumulé r, r2, // Distance entre deux particules et son carré rx, ry, rz, // Vecteur de distance entre deux particules ux, uy, uz, // Vecteur unitaire entre deux particules fx, fy, fz, // Force Ec, Ep, // Energie cinétique et energie potentielle mp float64 // Maximum des variations de la pente func main() { // Initialise la gestion de l'image
sdl.Init(sdl.INIT_VIDEO)
defer sdl.Quit() // initialise une fenêtre de fond noire 1000×1000 w, _ := sdl.CreateWindow("modèle", 0, 0, 1000, 1000, sdl.WINDOW_SHOWN) defer w.Destroy() ws, _ := w.GetSurface() var np, // Nombre de pas k, // Nombre d'iterations par pas i, j int // Indice des particules var p = new(Particule) // Référence d'une particule var q = new(Particule) // Référence d'une particule // Etats des particules à l'instant 0 P[0] = Particule{1, 400, 500, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0} P[1] = Particule{1, 600, 500, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0} P[2] = Particule{1, 500, 500, -150, -0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0} P[3] = Particule{1, 500, 500, 150, 0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0} t = 0.0 // Temps=0 dt = 1e-10 // intervalle initiale // Boucle sans fin, exécute une infinité de pas for np = 0; ; np++ { // Répéter 3 fois for k = 0; k < 3; k++ { // Pour chaque particule faire for i = 0; i < N; i++ { p = &P[i] // V1 = V0 + A0*dt + P0*dt^2/2 // X1 = X0 + V0*dt + A0*dt^2/2 + P0*dt^3/6 p.vx1 = p.vx0 + dt*(p.ax0+dt*p.px0/2) p.vy1 = p.vy0 + dt*(p.ay0+dt*p.py0/2) p.vz1 = p.vz0 + dt*(p.az0+dt*p.pz0/2) p.x1 = p.x0 + dt*(p.vx0+dt*(p.ax0+dt*p.px0/3)/2) p.y1 = p.y0 + dt*(p.vy0+dt*(p.ay0+dt*p.py0/3)/2) p.z1 = p.z0 + dt*(p.vz0+dt*(p.az0+dt*p.pz0/3)/2) } // Pour chaque particule for i = 0; i < N; i++ { p = &P[i] // f1 = 0 p.fx1 = 0 p.fy1 = 0 p.fz1 = 0 } // Pour chaque ensemble de deux particules {p,q} for i = 1; i < N; i++ { p = &P[i] for j = 0; j < i; j++ { q = &P[j] // R = q.X1-p.X1 // r = norm(R) // r2 = r*r // U = R/r // F = U*G*q.m*p.m/r2 // p.F = p.F + F // q.F = q.F - F rx = q.x1 - p.x1 ry = q.y1 - p.y1 rz = q.z1 - p.z1 r = math.Sqrt(rx*rx + ry*ry + rz*rz) r2 = r * r ux = rx / r uy = ry / r uz = rz / r fx = ux * G * q.m * p.m / r2 fy = uy * G * q.m * p.m / r2 fz = uz * G * q.m * p.m / r2 p.fx1 += fx p.fy1 += fy p.fz1 += fz q.fx1 -= fx q.fy1 -= fy q.fz1 -= fz } } mp = 1e-10 // mp = infiniment petit // Pour chaque particule for i = 0; i < N; i++ { p = &P[i] // A1 = F1/m // P0 = (A1-A0)/dt p.ax1 = p.fx1 / p.m p.ay1 = p.fy1 / p.m p.az1 = p.fz1 / p.m p.px0 = (p.ax1 - p.ax0) / dt p.py0 = (p.ay1 - p.ay0) / dt p.pz0 = (p.az1 - p.az0) / dt // R = P0-Ppred // r = Norm(R) // if r>mp then mp=r rx = q.px0 - p.px_ ry = q.py0 - p.py_ rz = q.pz0 - p.pz_ r = math.Sqrt(rx*rx + ry*ry + rz*rz) if r > mp { mp = r } } dt = math.Pow(24*epsilon/mp, 0.3333) // Recalcul de dt = (24*epsilon/mp)^(1/3) //--------------------------------------------------------------- } t += dt // Temps t=t+dt // Calcule de l'énergie totale Ep+Ec Ec = 0 for i = 0; i < N; i++ { p = &P[i] Ec += p.m * (p.vx1*p.vx1 + p.vy1*p.vy1 + p.vz1*p.vz1) / 2 } Ep = 0 for i = 1; i < N; i++ { p = &P[i] for j = 0; j < i; j++ { q = &P[j] rx = q.x1 - p.x1 ry = q.y1 - p.y1 rz = q.z1 - p.z1 r = math.Sqrt(rx*rx + ry*ry + rz*rz) Ep += -q.m / r } } Ep *= G // Pour chaque particule for i = 0; i < N; i++ { p = &P[i] // X0 = X1 // V0 = V1 // A0 = A1 // Ppred = P0 p.x0 = p.x1 p.y0 = p.y1 p.z0 = p.z1 p.vx0 = p.vx1 p.vy0 = p.vy1 p.vz0 = p.vz1 p.ax0 = p.ax1 p.ay0 = p.ay1 p.az0 = p.az1 p.px_ = p.px0 p.py_ = p.py0 p.pz_ = p.pz0 } // Pour chaque particule for i = 0; i < N; i++ { p = &P[i] // Afficher le point X ws.FillRect(&sdl.Rect{int32(p.x1), int32(p.y1), 1, 1}, 0xffffffff)
} x := math.Mod(float64(np), 10000.0) if x == 0 { w.UpdateSurface() // Actualiser l'affichage de la fenêtre fmt.Print("t=", t, " | ", Ep+Ec, "\n") } } sdl.Delay(60000) sdl.Quit() }

---- 27 novembre 2017 ----

 


Dominique Mabboux-Stromberg