Refroidissement de la Lune

Intérêt

La détermination d'un profil théorique est motivée par le fait que le profil observé dans le cas sans source semble auto-similaire , c'est à dire que la courbe semble identique et varie seulement avec un facteur lié au temps. C'est pourquoi on a tenté une approche avec une séparation des variables pour la température.

Méthode

On part du cas de l'équation de la chaleur adimensionnée sans termes sources définie par: $\frac{\partial T}{\partial t}=\frac{1}{\rho^2}\frac{\partial}{\partial\rho}\left({\rho^2}\frac{\partial T}{\partial\rho}\right)$ On fait donc l'hypothèse de la séparation des variables, d'où: $T(\rho,t)=X(\rho)Y(t)$ En insérant ce résultat dans l'équation de la chaleur, on obtient une équation de cette forme pour la fonction radiale: ${\rho^2}{\frac{{\partial^2}X}{\partial\rho^2}}+2\rho{\frac{\partial X}{\partial\rho}}+{\omega^2}{\rho^2}X=0$ le oméga provenant d'une constante que l'on retrouve en développant les équations et qui est arbitraire. En posant $X(\rho)={\rho^\nu}{z(\rho)}$ on obtient l'équation suivante: $\rho{\frac{\partial z}{\partial \rho}+{\rho^2}{\frac{{\partial^2}z}{\partial{\rho^2}}}+({\omega^2}{\rho^2}-{\nu^2})z=0$ si on impose $\nu=-\frac{1}{2}$ on retrouve cette forme pour les équations. Les solutions

sont les fonctions de Bessel donc

peut s'écrire: $X(\rho)=\sqrt{\rho}\left({a_0}{J_{-\frac{1}{2}}}(\omega\rho)+{a_1}{Y_{-\frac{1}{2}}(\omega\rho)\right)$ Les conditions aux limites de gradient de la température nul au centre et de variation de la température au bord proportionnelle (avec $\beta$ )à une température extérieure (voir cas précédent), impose comme condition pour le paramètre $\omega$ : $\frac{\beta}{\omega}=\tan{\omega}$ qui possède plusieurs solutions. Dans le programme, ces valeurs sont déterminées par la méthode de Newton jusqu'à un certain ordre car la fonction tangente est périodique et donc l'équation possède plusieurs solutions. Par cette méthode on recherche les zéros de la fonction correspondante à la différence des deux membres. On voit donc que la température est une somme de l'ensemble des solutions précédente avec la valeur de $\omega$ calculée par la méthode de Newton correspondante et en faisant le même travail pour la partie temporelle de la séparation des variables, on obtient finalement pour la température: $T(\rho,t)={\sum_{n=1}^{\infty}}{a_{n0}}{e^{-{\omega_n}^2 t}}{\sqrt{\frac{2}{{\omega_n}{\pi}}}}{\cos{{\omega_n}{\rho}}}$ En imposant les conditions initiales pour la température égale à 1 pour tout $\rho$ à t=0, on se retrouve avec un système linéaire à résoudre pour déterminer les coefficients ${a_{n0}}$ quand on tronque le développement à un certain ordre N. L'inversion de ce système par la méthode de la décomposition LU et de l'inversion qui en suit décrite plus haut permet de déterminer ces coefficients. Une fois tous les paramètres obtenus, il suffit de calculer la valeur de la température en chaque point et à un instant t dans le programmme par la formule précédente de la température théorique obtenue et de l'afficher à l'écran. Dans le cas théorique, le paramètre du flux $\beta$ est limité à certaines valeurs car la sommation effectuée dans le programme est réalisée des valeurs les plus faibles supposée des coefficients aux plus grands et quand ce paramètre est trop grand les calculs sont faussés. La courbe théorique sur le graphe est tracée en rouge et est proche du calcul numérique à certains moments et éloignée à d'autres, car notre hypothèse de séparation des variables n'est qu'une approximation résultante d'une constatation faite à partir du profil calculé par notre programme.

Ici, nous souhaitons considérer le cas où la Lune possède une source radiative interne, dont les particules radioactives vont sédimenter vers le centre de l'astre, refroidissant la surface et entraînant la formation d'une croûte. Nous nous placerons dans un cas isotrope où les différentes grandeurs physiques ne dépendront donc que du temps et du rayon

(coordonnées sphériques) par rapport au centre de l'astre.

Mise en équation du problème de sédimentation

En coordonnées cartésiennes :

Considérons les particules de la source radioactive. Deux phénomènes physiques sont alors à prendre en compte : d'une part la diffusion des particules, d'autre part la sédimentation des particules les plus massives.

Figure 3 :Elément de fluide considéré pour la détermination des équations du mouvement

Supposons un flux de particules dans la direction des $x$

positifs (figure 3). Supposons que notre fluide se comporte comme un gaz parfait de particules se déplaçant toutes avec une vitesse moyenne égale à la vitesse thermique des particules, $V$

. Appliquons le modèle du " $\frac{1}{6}$ '' (c'est-à-dire, considérons que les particules se déplacent toutes selon les axes). Dans ce cas, le flux de particules à travers la section $S$

pendant

vaut : $\frac{dN}{dt}=\frac{1}{dt}\left(\frac{1}{6}n(x,t)\cdot\left(V-v(x,t)\right)Sdt-\frac{1}{6}n(x+dx,t)\cdot\left(V+v(x+dx,t)\right)Sdt\right)$
$\frac{dN}{dt}=\frac{1}{6}\left(\left(n(x,t)-n(x+dx,t)\right)\cdot V-n(x,t)v(x,t)-n(x+dx,t)v(x+dx,t)\right)S$
où $dN$

est le nombre élémentaire de particules traversant $S$

pendant l'intervalle de temps $dt$

est la densité de particules en $x$

à l'instant $t$

, et

la vitesse des particules sédimentant.

Si l'on suppose que la distance typique séparant les particules en

(autrement dit

) est de l'ordre du libre parcours moyen $l_{pm}$ des particules, alors le flux de particules F(x)

par unité de temps à travers la surface

sera :
$F(x)=\frac{1}{6}\left(-Vl_{pm}\frac{\partial n}{\partial x}\left(x,t\right)-n(x,t)v(x,t)-n(x+dx,t)v(x+dx,t)\right)S$
que l'on choisira de réécrire sous la forme :
$F(x)=-D(x){\frac{\partial n}{\partial x}}S-{\alpha}n(x,t)v(x,t)S-{\alpha}n(x+dx,t)v(x+dx,t)S$
où D(x)

est le coefficient de diffusion des particules et peut dépendre de la position

, et ${\alpha}$ un coefficient que l'on peut assimiler au pourcentage de particules qui vont effectivement sédimenter dans la direction

avec la vitesse -v(x,t)

(en supposant, pour simplifier, que ${\alpha}$ ne dépend ni de la position, ni du temps).

Si de plus, on émet l'hypothèse que le nombre de particules sédimentant (ainsi que leur vitesse de sédimentation) dans la direction $-x$

est le même en $x$

, alors l'expression du flux devient :
$F(x)=-D(x){\frac{\partial n}{\partial x}}S-2{\alpha}n(x,t)v(x,t)S$
que l'on réécrira (en posant $\beta=2\alpha$ sans dimension) :
$F(x)=-D(x){\frac{\partial n}{\partial x}}S-{\beta}n(x,t)v(x,t)S$ (11)
Avec (11), on vérifie bien :
- d'une part, que le flux diffusif se fait bien des densités élevées vers les densités faibles ( $\frac{\partial n}{\partial x}<0$ et donc $-D(x){\frac{\partial n}{\partial x}}S>0$ lorsque les densités élevées sont au centre)
- d'autre part, que le flux de particules sédimentant est bien dirigé vers le centre, des densités faibles vers les densités élevées $-{\beta}n(x,t)v(x,t)S<0$ pour $\beta>0$ )
Note : Ceci permet donc de faire une simple vérification une fois la programmation faite, puisqu'on s'attendra alors, pour des valeurs extrêmes de $\beta$ , à obtenir soit de la sédimentation si $\beta>0$ , soit un phénomène de pseudodiffusion si $\beta<0$ ).

Figure 4 : Elément de fluide pour la conservation du nombre de particules dans le volume

A l'aide de la figure 4, nous allons établir l'équation différentielle pour n(x,t)

. Le nombre élémentaire de particules traversant une surface

pendant un intervalle de temps

est :

$d^2N={\frac{\partial n}{\partial t}}dtdV={\frac{\partial n}{\partial t}}dtSdx$ (12)
Le flux élémentaire

entre

est :
$dF=F(x+dx)-F(x)={\frac{\partial F}{\partial x}}dx={\frac{\partial}{\partial x}\left[-D(x){\frac{\partial n}{\partial x}}-{\beta}n(x,t)v(x,t)\right]}Sdx$
c'est-à-dire :
$dF=-\left[{\frac{\partial}{\partial x}}\left(D(x){\frac{\partial n}{\partial x}}\right)+{\beta}\frac{\partial\left(nv\right)}{\partial x}\right]Sdx$ (13)
La variation du nombre de particules dans le volume V=Sdx

entre

vaut :

n(x,t+dt)Sdx=n(x,t)Sdx+rentrant-sortant=n(x,t)Sdx+F(x)dt-F(x+dx)dt

(14)
En utilisant (12), (13), et (14), on obtient finalement l'équation du problème :
$\frac{\partial n}{\partial t}={\frac{\partial}{\partial x}}\left(D(x){\frac{\partial n}{\partial x}}\right)+{\beta}\frac{\partial\left(nv\right)}{\partial x}$

En coordonnées sphériques :

On considère un cas isotrope de sorte que $n\equiv n(r,t)$ et $v\equiv v(r,t)$ .

Figure 5 : Schéma descriptif du problème en coordonnées sphériques, avec symétrie sphérique

En reprenant le problème par analogie au cas de la section précédente (et en s'aidant de la figure 5), on trouve d'une part, que le flux à travers une surface S(r)

par unité de temps vaut :
$F(r)=-D(r){\frac{\partial n}{\partial r}}S(r)-{\beta}n(r,t)v(r,t)S(r)$
et que :
$d^2N=\frac{\partial n}{\partial t}S(r)dt=F(r)S(r)dt-F(r+dr)S(r+dr)dt=entrant-sortant$
où S(r)

en sphérique vaut : $S(r)=4\pi r^2$ .
Note : Rigoureusement, il faudrait écrire dans les équations précédentes dS(r)

et non

. Dans le cas général, $dS=dS(r,\Omega)=r^2d\Omega$ . Cependant, en considérant un cas isotrope, aucune grandeur ne dépend de l'angle solide, et par conséquent, on peut intégrer sur tout l'angle solide, ce qui revient à n'intégrer que

selon $\Omega$ ce qui donne $S(r)=4\pi r^2$ .

Finalement, après avoir remplacé S(r)

par son expression puis simplifié par $4\pi$ , on obtient :
$\frac{\partial n}{\partial t}=-\frac{1}{r\lyxmathsym{\texttwosuperior}}\frac{\partial\left(r^2F\right)}{\partial r}$
Ce qui conduit à l'équation suivante :
$\frac{\partial n}{\partial t}= \frac{1}{r^2}\frac{\partial}{\partial r}\left(r^2D(r)\frac{\partial n}{\partial r}\right)+\frac{\beta}{r^2}\frac{\partial (r^2nv)}{\partial r}$ (15)

A présent, nous allons dédimensionner l'équation (15) avant de lui appliquer l'algorithme de Crank-Nicholson. Tout d'abord, considérons deux constantes $D_{0}$ et $v_{0}$ , grandeurs dimensionnées caractéristiques du coefficient de diffusion D(r)

et de la vitesse de sédimentation des particules v(r,t)

, telles que
$D(r)={D_{0}}\times f(r)$
et
$v(r,t)={v_{0}}\times g(r,t)$
où f(r)

sont sans dimension.
Puis, appliquons les changements de variables suivants :
r=Rr'

$t={\tau}t'$
où $R$

et $\tau$ sont des grandeurs caractéristiques du problème (ici,

correspondrait au rayon de la source, et $\tau$ au temps caractéristique d'évolution de la densité). L'équation (15) devient donc :
$\frac{\partial n}{\partial t'}= \frac{1}{r'^2}\frac{\tau D_0}{R^2}\frac{\partial}{\partial r'}\left( r'^2f\frac{\partial n}{\partial r'}\right)+\frac{\beta v_0\tau}{Rr'^2}\frac{\partial (r'^2ng)}{\partial r'}$
avec $f\equiv f(r')$ et $g\equiv g(r',t')$ .

Enfin, choisissons $D_{0}$ ,

, et $\tau$ telles que $\frac{\tau D_{0}}{R^2}=1$ . Posons $\gamma=\frac{\beta v_{0}\tau}{R}$ (notons que, comme $D_{0}$ est en m^2/s

et que $\beta$ est sans dimension, l'homogénéïté pour les deux égalités est bien respectée). L'équation adimensionnée est donc :
$\frac{\partial n}{\partial t'}=\frac{1}{r'^2}\frac{\partial}{\partial r'}\left(r'^2f\frac{\partial n}{\partial r'}\right)+\frac{\gamma}{r'^2}\frac{\partial (r'^2ng)}{\partial r'}$ (16)
Par commodité, on notera pour la suite r=r'

dans les équations (

sans dimension).

Algorithme de Crank-Nicholson pour la sédimentation

On considère que le coefficient de diffusion et la vitesse moyenne de sédimentation des particules sont constants.

Ecriture de l'équation sous forme algorithmique

Pour ce cas, cela signifie que l'on prendra f(r)=1

, de sorte que l'équation (16) s'écrive :
$\frac{\partial n}{\partial t}={\frac{1}{r²}}{\frac{\partial}{\partial r}}\left(r²{\frac{\partial n}{\partial r}}\right)+\frac{\gamma}{r²}\frac{\partial (r²n)}{\partial r}$ (17)
Pour appliquer Crank-Nicholson, nous fixons un pas spatial

et un pas temporel

. Ensuite, il va s'agir d'établir l'expression discrétisée de l'équation précédente en $(i,j+\frac{1}{2})$ (indice spatial, indice temporel), $i$

étant des entiers. On a alors :
${\frac{\partial n}{\partial t}}^{i,j+\frac{1}{2}}=\frac{n^{i,j+1}-n^{i,j}}{k}$ (18)
De la même façon, nous aurons pour une dérivée première spatiale : ${\frac{\partial n}{\partial r}}^{i+\frac{1}{2},j}=\frac{n^{i+1,j}-n^{i,j}}{h}$
Note : Si jamais l'indice $i+\frac{1}{2}$ (respectivement $j+\frac{1}{2}$ ) apparaît pour une variable, étant donné que l'on a discrétisé en nombres entiers de

et de

, alors on calculera la valeur en ce point comme étant la valeur moyenne du point suivant i+1

(respectivement j+1

) et du point précédent

(respectivement

).
Exemple :
$n^{i+\frac{1}{2},j}=\frac{1}{2}\left(n^{i+1,j}+n^{i,j}\right)$
ou encore
$n^{i,j+\frac{1}{2}}=\frac{1}{2}\left(n^{i,j+1}+n^{i,j}\right)$

de même pour toute dérivée partielle d'ordre

par rapport à une variable

, on aura :
$\frac{\partial^k n^{i+\frac{1}{2},j}}{\partial x^k}=\frac{1}{2\times p}\left({\frac{\partial^k n^{i+1,j}}{\partial x^k}}+{\frac{\partial^k n^{i,j}}{\partial x^k}}\right)$
où

est le pas relatif à la variable

(dans notre cas, p=k

, ou

Notons que dans l'équation (17) , le premier terme du membre de droite de l'égalité est en fait le laplacien en coordonnées sphériques. Pour des raisons de commodité, nous allons utiliser l'autre forme du laplacien sphérique d'une fonction f(r)

qui est :
$\Delta (f) \equiv {\frac{1}{r}}{\frac{\partial^2 }{\partial r²}\left(rf\right)}$
Dès lors, on obtient que :
$\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2\frac{\partial n}{\partial r}\right)^{i,\,j}=\frac{1}{r}\frac{\partial ^2}{\partial r^2}(rn)^{i,\,j}=\frac{r^{i+1}n^{i+1,\,j}-2r^in^{i,\,j}+r^{i-1}n^{i-1,\,j}}{r^ih^2}$

Il s'ensuit que le laplacien de

pris en $(i,j+\frac{1}{2})$ s'écrit :
${\frac{1}{r}}{\frac{\partial^2 }{\partial r²}\left(rn\right)}^{i,j+\frac{1}{2}}={\frac{1}{2}}\left({\frac{1}{r}}{\frac{\partial^2 }{\partial r²}\left(rn\right)}^{i,j+1}+{\frac{1}{r}}{\frac{\partial^2 }{\partial r²}\left(rn\right)}^{i,j}\right)$ $=\frac{r^{i+1} n^{i+1,j+1} -2r^i n^{i,j+1} +r^{i-1} n^{i-1,j+1} + r^{i+1} n^{i+1,j} -2r^i n^{i,j} +r^{i-1} n^{i-1,j}}{r^i h^2}$ (19)

De même, on trouve que :
$\left(\frac{\gamma}{r^2} {\frac{\partial (r^2n)}{\partial r}}\right)^{i,j+\frac{1}{2}}=\frac{\gamma}{(r^2)^{i}}\frac{(r^2n)^{i+1,j+1} -(r^2n)^{i-1,j+1} +(r^2n)^{i+1,j} -(r^2n)^{i-1,j}}{4h}$ (20)
A présent, nous allons poser $r^{i}=ih$ et $\rho=\frac{k}{2h^2}$ .
Note : $(r^2n)^{i}=(r^{i})^2$ ainsi, on voit facilement que $\frac{(r^2n)^{i+1,j+1}}{{r^{i}}^2}=\left(\frac{i+1}{i}\right)^{2}$
En prenant (18), (19) et (20) que l'on injecte dans (17), on obtient alors l'équation suivante :
$-{\frac{i-1}{i}\left(1-\frac{i-1}{2i}\gamma h\right)\rho n^{i-1,j+1}+\left(1+2\rho \right)n^{i,j+1}-\frac{i+1}{i}\left(1+\frac{i+1}{2i}\gamma h\right)\rho n^{i+1,j+1}}$ $={\frac{i-1}{i}\left(1-\frac{i-1}{2i}\gamma h\right)\rho n^{i-1,j}+\left(1-2\rho \right)n^{i,j}+\frac{i+1}{i}\left(1+\frac{i+1}{2i}\gamma h\right)\rho n^{i+1,j}}$ (21)

Détermination des conditions aux limites :

Pour résoudre complètement l'équation de départ (17), il ne nous manque plus que les conditions aux limites.

Condition limite en :
Pour obtenir une condition en , nous allons procéder de la même façon que dans la cas sans source, c'est-à-dire en intégrant l'équation (17) sur une sphère infinitésimale de centre 0 et de rayon $\epsilon$ . Nous considèrerons qu'en , le gradient de densité est nul pour n'importe quel temps donc quelque soit ( ${\frac{\partial n}{\partial r}}^{0,j}=0$ ).
$\int{\frac{\partial n}{\partial t}dV}=\int{{\frac{1}{r^2}}{\frac{\partial}{\partial r}}\left(r^2{\frac{\partial n}{\partial r}}\right)}dV+\int{\frac{\gamma}{r^2}\frac{\partial (r^2n)}{\partial r}}dV$
où $dV=r² sin(\theta)drd\theta d\phi$ , ce qui donne, lorsque l'on considère un problème à symétrie sphérique (après simplification par $4\pi$ ): $\int_{r=0}^{r=\epsilon} {\frac{\partial n}{\partial t}r^2dr}=\int_{r=0}^{r=\epsilon} {{\frac{\partial}{\partial r}}\left(r^2{\frac{\partial n}{\partial r}}\right)}dr+\int_{r=0}^{r=\epsilon} {{\gamma}\frac{\partial (r^2n)}{\partial r}}dr$    (22)
Notons que la seule différence avec le cas de diffusion seule (analogie avec l'équation de chaleur mais en remplaçant par ) est la présence d'un second terme à intégrer dans le membre de droite de l'équation ci-dessus. Nous allons donc ne nous occuper que de l'intégration de ce terme et de son écriture via Crank-Nicholson, l'intégration des deux autres termes ayant déjà été faite dans la cas de l'équation de la chaleur avec diffusion seule.
On a :
$\int_{0}^{\epsilon} {\gamma \frac{\partial (r^2n)}{\partial r}}dr=\gamma \int_{0}^{\epsilon} {\partial (r^2n)}=\gamma {\left[r^2n\right]}_{0}^{\epsilon}=\gamma \epsilon ^2 n(\epsilon,t)$    (23)
En remplaçant (23) dans (22) et en utilisant le résultat pour la diffusion, on trouve alors :
$\left({{\frac{\partial n}{\partial t}}\vline}_{r=0}+{{\frac{\partial n}{\partial t}}\vline}_{r=h}\right){\frac{\epsilon ^3}{6}}=\epsilon ^2 {{\frac{\partial n}{\partial r}}\vline}_{r=\epsilon}+\gamma \epsilon ^2 n(\epsilon,t)$
En reprenant la discrétisation précédente, et en choisissant $\epsilon=\frac{h}{2}$ (milieu de la première maille, on obtient après une petite séance de gymnastique mathématique, la condition au centre :
$\left(1+12\rho -\frac{\gamma h}{2} \rho \right)n^{0,j+1} +\left(1-12\rho-\frac{\gamma h}{2} \rho \right)n^{1,j+1}$ $=\left(1-12\rho +\frac{\gamma h}{2} \rho \right)n^{0,j} +\left(1+12\rho +\frac{\gamma h}{2} \rho \right)n^{1,j}$    (24)
(où $\rho$ vaut toujours $\rho=\frac{k}{2h^2}$ ).
Seconde condition aux limites :
Comme deuxième condition, nous prendrons la conservation de la masse des particules radioactives de la source, :
$\int_V \rho (\overrightarrow{r})dV=M$
Comme nous considérons un cas isotrope, nous avons $\rho(\overrightarrow{r})=\rho(r)$ et donc en intégrant sur les deux angles $\theta$ et $\phi$ :
$\int_{r=0}^{r=R} \rho (r)r^2dr=\frac{M}{4\pi}$ ou encore, avec le changement de variable $r\equiv Rr$ :
$\int_{r=0}^{r=1} \rho (r)r^2dr=\frac{M}{4\pi R^3}$
Ce qui, discrétisé de la même façon que précédemment et en utilisant la méthode des trapèzes peut se réécrire : $\sum_{i=0}^{n-1} \left(\frac{\rho^{i}(r)r_{i}^{2}+\rho^{i+1}(r)r_{i+1}^{2}}{2}\right)\left(r_{i+1}-r_i\right)=\frac{M}{4\pi R^3}$
A partir de cette expression, on choisira de mettre la deuxième condition sous la forme suivante, ( $n(r)=\frac{\rho(r)}{Am_{p}}$ , où correspond au nombre de nucléons de l'élément radioactif et $m_{p}$ à la masse d'un proton ; rappelons que $r_{i}=ih$ ) donnant une condition en :
$n(R,t)=n^{n,j}=\frac{M/Am_p}{2\pi n² h^3 R^3}-\sum_{i=0}^{n-2} \frac{i² n^{i,j}+(i+1)² n^{i+1,j}}{n²}-\frac{\left(n-1\right)²}{n²}n^{n-1,j}$ (23)
On remarque que dans la somme, les termes en $i^2n^{i,j}$ et en $(i+1)^2n^{i+1,j}$ se ressemblent, reprenons-les et essayons d'arranger la somme différemment : $\sum_{i=0}^{n-2} \frac{i^2 n^{i,j}+(i+1)^2 n^{i+1,j}}{n^2}=\sum_{i=0}^{n-2} \frac{i^2 n^{i,j}}{n^2}+\sum_{i=0}^{n-2} \frac{(i+1)^2 n^{i+1,j}}{n^2}$
Prenons le membre en et appliquons-lui le changement de variable suivant : . Ce terme se réécrit donc : $\sum_{i=0}^{n-2} \frac{i^2 n^{i,j}}{n^2}=\frac{0^2n^{0,j}}{n^2}+\sum_{i=1}^{n-2} \frac{i^2 n^{i,j}}{n^2}=\sum_{i=1}^{n-2} \frac{i^2 n^{i,j}}{n²}$
$\sum_{i=1}^{n-2} \frac{i^2 n^{i,j}}{n^2}=\sum_{i'=0}^{n-3} \frac{(i'+1)^2 n^{i'+1,j}}{n^2}$
ce qui donne en utilisant le changement : $\sum_{i'=0}^{n-3} \frac{(i'+1)^2 n^{i'+1,j}}{n^2}=\sum_{i=0}^{n-3} \frac{(i+1)^2 n^{i+1,j}}{n^2}$
En remarquant que : $\frac{(n-1)^2}{n^2}n^{n-1,j}=\frac{\left((n-2)+1\right)^2}{n^2}n^{(n-2)+1,j}$
on obtient : $\frac{(n-1)^2}{n^2}n^{n-1,j}+\sum_{i=0}^{n-3} \frac{(i+1)^2 n^{i+1,j}}{n^2}=\sum_{i=0}^{n-2} \frac{(i+1)^2 n^{i+1,j}}{n^2}$
D'où finalement, on en déduit l'expression simplifiée de (23) :
$n(R,t)=n^{n,j}=\frac{M/Am_p}{2\pi n^2 h^3 R^3}-\frac{2}{n^2}\sum_{i=0}^{n-2} {(i+1)^2 n^{i+1,j}}$
que l'on pourra réécrire () :
$n(R,t)=n^{n,j}=\frac{M/Am_p}{2\pi n^2 h^3 R^3}-\frac{2}{n²}\sum_{i=0}^{n-1} {i^2 n^{i,j}}$
c'est-à-dire :
$n(R,t)=n^{n,j}=\frac{M/Am_p}{2\pi n^2 h^3 R^3}-\frac{2}{n^2}\sum_{i=1}^{n-1} {i^2 n^{i,j}}$ (25)
puisque $0*n^{0,j}=0$ .
En utlisant (21), (24) et (25), on peut alors résoudre l'agorithme numériquement.

Note : Un petit calcul rapide nous permet d'avoir une idée de la densité d'élément radioactif que l'on peut espérer trouver pour une masse donnée de cet élément. Partons d'un état initial uniforme, dans ce cas, dans ce cas, on a directement (en utilisant le fait que $\sum_{0}^{n-1}{i²}=(n-1)n(2n-1)/6\sim n^3/3$ , suffisament grand) :
$n_{(/m^3)}^{i,0}\sim\frac{9M}{4\pi R^3h^3n^3Am_p}=\frac{9M}{4\pi R^3Am_p}\sim3.4\cdot10^{5}M_{(kg)}$
lorsque l'on considère de l'uranium () avec $R_{lune}=1737.4km$ , étant la masse totale d'uranium à l'intérieur de la Lune que l'utilisateur pourra choisir.