# Euler algorithm for the computation of the evolution of planetary rings

Hi,

I am currently working on a program to compute the evolution of planetary rings over large time scale. For the moment I’m trying to replicate the results of an algorithm described p6 of this article :

https://hal.science/hal-00676215/document

While the article is using a Runge-Kutta method I choose to use an explicit Euler method.
The problems I encounter are :

• the evolution is too slow
• As the time increases, my Gaussian function shifts on the left while its center should not move.

Here is the code :

``````import numpy as np
import matplotlib.pyplot as plt

G = 6.67430e-11   #gravitational constant
Ms = 5.6834e26    #Saturn's mass in kg
Mm = 3.84e19       #Mima's mass en kg

rhoP = 900.0                      #particle density
rp = 1.0                              #particle radius in m
mp = rhoP*4/3*np.pi*rp**3 #particle mass in kg
nu = 0.1                             #constant viscosity in  m2.s-1
dispR = np.sqrt(G*mp/rp)  #radial dispertion en m.s-1

Rmin, Rmax = 60e6, 140e6
dr = 100e3                                                                     #element size in m
Rmid = np.arange(Rmin, Rmax+dr, dr)                         #position of elements mid points
RInf = np.arange(Rmin-dr/2, Rmax, dr)                 #position of elements lower limit
RSup = np.arange(Rmin+dr/2, Rmax+1+dr/2, dr) #position of elements upper limit
Nr = int((Rmax-Rmin)/dr)          #number of elements for the finite element method

Tmin, Tmax = 0, 1e5*365.25*24*3600
dt = 0.5*(dr**2)                                 #time step
T = np.arange(Tmin, Tmax, dt)
Nt = int((Tmax-Tmin)/dt)                  #number of time steps
#omega = np.sqrt(G*Ms/(R**3))      #orbital frequency in Hz

#Initial state
def repartSigma(centre, width):
return 6*np.exp(-(Rdemi-centre)**2 / (2*(largeur**2)))

#evolution
def f(nu, sigma) :
return nu*sigma*np.sqrt(Rdemi)

def calcDerivF(nu, sigma) :
fct = f(nu, sigma)
derivF = np.zeros_like(Rdemi)
for i in range (Nr) :
derivF[i] = (fct[i+1]-fct[i])/(dr)
return derivF

def calcFlux(R, nu, sigma) :
derivF = calcDerivF(nu, sigma)
flux = 6*np.pi*np.sqrt(R)*derivF
return flux

def calcDerivSigma(nu, sigma):
fluxInf = calcFlux(RInf, nu, sigma)
fluxSup = calcFlux(RSup, nu, sigma)
fluxInf[0] = fluxSup[0]
fluxSup[-1] = fluxInf[-1]
deriveeSigma = (fluxSup - fluxInf)/(2*np.pi*Rmid*dr)
return derivSigma

def euler(nu, centre, witdh) :
sigma = repartSigma(centre, width)
for i in range (Nt) :
derivSigma = calcDerivSigma(nu, sigma)
sigma += derivSigma*dt
print (i/Nt)    #progress check
return sigma

def calcMass(R, sigma) :
mass = 0
for i in range (Nr) :
mass += sigma[i]*2*np.pi*R[i]*dr
return masse

sigmaIni = repartSigma(110e6, 3e5)
sigmaFin= euler(0.1, 110e6, 3e5)

plt.figure()
plt.plot(Rdemi, sigmaIni, "r", label='0 years')
plt.plot(Rdemi, sigmaFin, "b", label=f'{Tmax/(365.25*24*3600)} ')
plt.title(f'{Tmax/(365.25*24*3600)} years')
plt.xlabel('Distance à Saturn (100000 km)')
plt.ylabel(r'\$\Sigma\$  (\${10}^{4}kg.m^{-2}\$)')
plt.xlim(100e6, 120e6)
plt.ylim(0, 8)
plt.legend()
plt.show()
``````

You seem to be asking a question that need mathamatics knowledge rather then python knowledge to help you?

Can you frame your problem as a python question?

Is there some place where `Rdemi`, `largeur` get defined?

Yes, this is what happens if you use “explicit Euler”. It isn’t stable. I would imagine that is part of why the article uses Runge-Kutta. You could use a variation, called “implicit Euler”, which should avoid this effect.