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.