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 un ordinateur 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 ni même d'en évaluer 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) 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 dans un espace à 3 dimensions.

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 :

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.

3) 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`.

4) 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)`

5) 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 + ...`

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`. L'approximation que nous faisons consiste à considérer que pendant chaque intervalle de temps `Deltat`, l'accélération de chaque particule évolue linéairement, c'est à dire que leurs pentes sont constantes. Autrement dit, la trajectoire est localement du 3-ième degrès. Voir l'approximation avec une trajectoire localement du 4-ième degrès au chapitre Modèle de Newton - Approximation à dévivée quatrième 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(x_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`

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!) + ...`

Des considérations générales permettent de supposer que les valeurs 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.

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

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

Et on l'appelle, le n-ième terme de la dérivée de Taylor.

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

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

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

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

Mais lorsque `dvec(p_0)"/"dt` devient nulle, il faut alors aller checher le terme suivant, le cinquième terme, qui n'est plus négligeable :

`(d^2vecp_0)/(dt^2)(Deltat^5)/(5!)`

Mais nous n'avons pas de connaissance sur `d^2vecp_0"/"dt^2` . C'est pourquoi Il est judicieux de le majorer par la valeur supérieur `zeta` rencontrée depuis le début de la trajectoire.

`epsilon  ≈  |vec(p_1)-vec(p_0)|(Deltat^3)/(4!) + zeta(Deltat^5)/(5!)`

Puis on veut pouvoir adapter `Deltat` à chaque pas, pour maintenir un ordre d'erreur constant `epsilon`. On résoud de façon approché l'équation :

` epsilon/(Deltat^3)≈ |vec(p_1)-vec(p_0)|/(4!) + zeta(Deltat^2)/(5!)`

`Deltat := root(3)((120epsilon)/(5|vec(p_1)-vec(p_0)|+zeta Deltat^2))`

Où l'on utilise l'ancienne valeur de `Deltat` pour calculer la nouvelle valeur approchée de `Deltat`.

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 := root(3)((120epsilon)/(5sf"max"(|vec(p_1)-vec(p_0)|)+zeta Deltat^2))`

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_("-"1))`, 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_("-"1))`

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 pour chaque particule son `zeta` maximum rencontré.

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

On recalcule l'intervalle de temps : `Deltat := root(3)((120epsilon)/(5sf"max"(|vec(p_0)-vec(p_("-"1))|)+zeta Deltat^2))`

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 3 fois systématiquement.

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

---- 22 mars 2021 ----

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