Modèle de Newton
Approximation à dévivée quatrième constante

1) Introduction

Algorithme de calcul des trajectoires des particules avec comme approximation, des trajectoires localement du `4`-ième degrés. L'approximation, avec des trajectoires localement du `3`-ième degrès, est vue au chapitre Modèle de Newton

2) 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, qui défini l'accélération des particules :

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

`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`.
`r_(ij)`
: Distance entre la particule n°`i` et la particule n°`j`.

3) Cinématique

Une particule est une masse ponctuelle. La position de la particule se note `vecx`. Les différentes dérivés par rapport au temps de `vecx` définissent la vitesse `vec v=vec x^((1))`, l'accélération `vec a = vec x^((2))`, la troisieme dérivée `vec p = vec x^((3))` appelée pente et la quatrième dérivée `vec s=vec x^((4))`

Position :
`vec x`
`=`
`vecx`
Vitesse :
`vec v`
`=`
`(dvecx)/(dt)`
Acceleration :
`vec a`
`=`
`(d^2vecx)/(dt^2)`
3-ième dérivée :
`vec p`
`=`
`(d^3vecx)/(dt^3)`
4-ième dérivée :
`vec q`
`=`
`(d^4vecx)/(dt^4)`

A l'instant `t"="0`, la particule se situe à la position `vecx_0` avec une vitesse `vecv_0`, une accélération `veca_0`, une dérivée troisième `vecp_0`, et une dérivée quatrième `vecq_0`.

La trajectoire localement est une trajectoire du `4`-ième degrés. La dérivée quatrième est constante.

La loi de Newton nous permet de calculer l'accéleration en chaque point. Il est alors nécessaire d'avoir deux points pour calculer la pente `vecp`, et finalement il est nécessaire d'avoir `3` points pour calculer la dérivée de la pente `vecq`.

La position d'une particule fonction du temps `t` se note `vecx(t)`, et se développe en 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!) + vec(q_0)t^4/(4!) + ...`

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

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

`veca(t) = vec(a_0) + vec(p_0)t + vec(q_0)t^2/(2!) + ...`

`vecp(t) = vec(p_0) + vec(q_0)t + ...`

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 `Delta` :

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

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 quatrième terme, est du même ordre de grandeur que le cinquième terme en général qui est :

`epsilon  ≈  |(dvec(q_0))/(dt)|(Delta^5)/(5!)`

`epsilon  ≈  |vec(q_1)-vec(q_0)|/(Delta) (Delta^5)/(5!)`

`epsilon  ≈  |vec(q_1)-vec(q_0)|(Delta^4)/(5!)`

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

`(d^2vecq_0)/(dt^2)(Delta^6)/(6!)`

Mais nous n'avons pas de connaissance sur `d^2vecq_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(q_1)-vec(q_0)|Delta^4/(5!) + zeta(Delta^6)/(6!)`

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

`epsilon/(Delta^4) ≈  |vec(q_1)-vec(q_0)|/(5!) + zeta(Delta^3)/(6!)`

`Delta := root(4)(((6!)epsilon)/(6|vec(q_1)-vec(q_0)|+zetaDelta^3))`

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

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

`Delta := root(4)(((6!)epsilon)/(6sf"max"(|vec(q_1)-vec(q_0)|)+zetaDelta^3))`

5)) 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. L'approximation que nous faisons consiste à considérer que pendant chaque intervalle de temps, chaque particule suit une courbe du `4`-ième degrès.

A l'instant `0`, la particule se situe à la position `vecx_0` avec une vitesse `vecv_0`, une accélération `veca_0`, une dérivée troisième `vecp_0` et une dérivée quatrième `vecq_0`. Ces `4` paramètres définissent une courbe du `4`-ième degrés.

Puis à l'instant `Delta`, la particule se situe à la position `vecx_1` :

`vecx_1 = vecx_0 + vecv_0 Delta + veca_0 (Delta^2)/2 + vecp_0(Delta^3)/(3!)+ vecq_0(Delta^4)/(4!)`

On peut alors calculer l'accélération produites par le système d'équation différentielles, en `vecx_1`, que l'on nomme `vecA_1`. Si elles ne coïncident pas avec la valeur `veca_1`, c'est qu'il faut ajuster le paramètres constant de la trajectoire qui est la dérivée quatrième `vec q_0`.

`vecA_1 = veca_0 + vecp_0 Delta+vecq_0(Delta^2)/2`

`vecq_0 = 2/(Delta^2)(vecA_1-veca_0-vecp_0 Delta)`

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Pour cela, on applique la méthode de résolution approchée de Newton vectorielle, qui généralise la méthode de newton non vectorielle.

 

 

 

 

5) Méthode de résolution approchée de Newton

Étant donné une fonction `f`, qui appliquée à `x` produit `f(x)`, ce qui se note par l'expression `x|->f(x)`, nous voulons trouver de façon approchée une racine `x`, c'est à dire telle que `f(x)=0`. On pose une valeur `x` proche de la racine, et on répète un petit nombre de fois l'itération suivantes :

`x := x - f(x)/(f’(x))`

Le symbole `:=` indique que la valeur de la variable `x` va changer et prendre pour valeur l'évaluation du second membre. Il faut considérer cette ligne de code `x := x - f(x)"/"f’(x)` comme une instruction d'un programme et non comme une équation, où la nouvelle valeur de `x` est calculée en fonction de l'ancienne valeur de `x`.

La démonstration de la méthode de Newton se fait à l'aide de la dérivée et des expressions différentielles. On part d'une valeur `x` proche de la racine, et on cherche la variation `dx` telle que `x"+"dx` soit la racine, c'est à dire telle que `f(x"+"dx)=0`. Selon le calcul différentiel, nous avons :

`f’(x) = (f(x"+"dx) - f(x))/dx`
`f(x"+"dx) = f(x) + f’(x) dx`

On cherche `dx` tel que :

`f(x"+"dx)=0`
`f(x) + f’(x) dx =0`
`dx = - f(x)/(f’(x))`

On en déduit que la valeur `x"+"dx`, qui est égale à `x - f(x)"/"f’(x)`, constitue une approximation de la racine autant meilleur que `dx` est petit.

6) Méthode de résolution approchée de Newton vectorielle

La méthode se généralise aux fonctions vectorielles. La notation vectorielle permet de regrouper des calcules identiques s'appliquant sur plusieurs coordonnées en parallèles. Une première approche simple consiste à représenter les vecteurs sous formes de vecteur colonne. Et de considérer deux sortes de produit, le produit général qui correspond à l'application d'une fonction quelconque sur ses arguments, et le produit matriciel qui correspond à l'application d'une fonction linéaire sur ses arguments :

Néanmoins il convient de préciser quelques conventions d'écriture. Dans le texte, pour condenser l'écriture, nous écrivons les vecteurs en ligne bien que conceptuellement ils sont considérés comme étant des vecteurs colonnes. Ainsi le vecteur `vec v=(x,y)` s'écrit en équation :

`vec v = ((x),(y))`

Une fonction vectorielle `vec f`, dans notre cadre générale, est une fonction de `RR^2->RR^2`. Elle s'applique à un vecteur appartenant à `RR^2` pour produire un vecteur appartenant à `RR^2`. Et par convention on simplifie l'écriture du double niveau de parenthèse en un seul niveau de parenthèse qui s'écrit dans le texte par ` vec f{:"("(x,y)")":} = vecf"("x,y")"` et dans l'équation par :

`vec f{:(((x),(y))):} = vec f{:((x),(y)):}`

La fonction `f` se décompose en deux composantes qui sont deux formes : `A` de `RR->RR^2` et `B` de `RR-> RR^2`. La fonction vectorielle `vec f=(A,B)` appliquée au vecteur `v` produit le vecteur `vecf"("vecv")"`. Et là, il ne s'agit pas d'un produit matriciel mais d'un appel de fonction général :

`vec f{:((x),(y)):} = ((A),(B)){:((x),(y)):} = ((A{: ((x),(y)):}),(B{: ((x),(y)):}))`

`vec f(vecv) = ((A),(B))(vec v)= ((A(vecv)),(B(vecv)))`

`vec f = ((A),(B))`

ce qui se note par l'expression `(x,y)|->(F(x,y),G(x,y))`, nous voulons trouver de façon approchée une racine `(x,y)`, c'est à dire telle que `f(x,y)=(0,0)`. On pose une valeur `(x,y)` proche de la racine, et on répète un petit nombre de fois l'itération suivantes :

`(x,y) := (x,y) - f’(x)^-1 × f(x,y)

Le symbole `:=` indique que la valeur de la variable `x` va changer et prendre pour valeur l'évaluation du second membre. Il faut considérer cette ligne de code `x := x - f(x)"/"f’(x)` comme une instruction impérative d'un programme et non comme une équation.

La démonstration de la méthode de Newton se fait à l'aide des expressions différentielles et de la dérivée. On part d'une valeur `x` proche de la racine, et on cherche la variation `dx` telle que `x"+"dx` soit la racine, c'est à dire telle que `f(x"+"dx)=0`. Selon le calcul différentiel, nous avons :

`f’(x) = (f(x"+"dx) - f(x))/dx`
`f(x"+"dx) = f(x) + f’(x) dx`

On cherche `dx` tel que :

`f(x"+"dx)=0`
`f(x) + f’(x) dx =0`
`dx = - f(x)/(f’(x))`

On en déduit que la valeur `x"+"dx`, qui est égale à `x - f(x)"/"f’(x)`, constitue une approximation de la racine autant meilleur que `dx` est petit.

Etant donné une fonction `f` de `RR^2` dans `RR^2`, on cherche des racines c'est à dire des vecteur `vecx in RR^2` tel que `f(vecx)= vec0`.

 

 

---- 23 mars 2021 ----

 


Dominique Mabboux-Stromberg