A friend and I were discussing gravity including retarded potentials from the finite speed of light. This is not general relativity but more like special relativity applied to gravity.
We had the idea to simulate some “galaxies” using Newton’s laws, and then again using retarded gravity, and see what differences occurred (statistically).
In retarded gravity, the force on a given particle at x_i comes from the newton’s equations except evaluated at a given time in the past such that the distance between x_i and x_j is precisely the distance needed for light to have just arrived:
F_{ij} = \frac{Gm_im_j}{R^3} (x_i(t) - x_j(t-\tau_{ij}))
where \tau_{ij} is such that |(x_i(t) - x_j(t-\tau_{ij}))| = c \tau_{ij}
notice that there’s a different delay for each pair of particles.
Now I figure the amazing @ChrisRackauckas and friends have probably implemented delay differential equations with variable delay…
Delay Differential Equations · DifferentialEquations.jl clearly gives us this impression. But it’s not clear to me how to specify a system where the delay for each pair of particles is dependent on solving an algebraic equation like the second equation above.
Also, to initialize this system, it seems we need some history for every particle. Our idea was to run it as Newton’s gravitation for a while… then use the position histories for initializing the retarded gravitation.
Wondering if someone could show me how to specify a problem like this with say 3 particles… I could probably take it from there.
Is this a delay differential-algebraic equation? (the \tau values obeying algebraic equations)? Is that even a problem type that DifferentialEquations.jl can handle?
Yes it’s an interesting discussion, I’ll forward to my friend as well. Any thoughts on how to simulate these equations in Julia (I think actually my friend derived more complete form of equations which have additional parts) but in general the key is this \tau_{ij} which is a continuously dynamically changing quantity obeying the algebraic equation, and which makes the force more complicated.
Easy for you to say Is there an example you can point to? or which section of the docs I should read? Happy to give it a whirl but probably need a little more on-ramp.
Will take a look thanks! The DAE part I knew DifferentialEquations.jl could handle… the part where it’s a combined delay and algebraic was where I got confused. Will come back after giving a read and trying it out.
Do you need to use specific Julia packages? That seems simple to implement using point particles and a integrator with finite-differences, it just depends on putting a branch in the function that computes the forces.
The DDE solver is just the ODE solver with extra steps. The paper will come out soon, but it’s just OrdinaryDiffEq wrapped inside of itself with a compiler trick.
I might wind up having to do a specialized integrator for this in order to keep memory size down (for example it’d be nice to do 1_000_000 particles but a trillion element explicit mass matrix is obviously non-starter) but if I can get it going with debugged standardized libraries for at least a few hundred particles it would be a useful baseline.
I haven’t had a chance to look carefully yet, but it looks like I can construct a DDEProblem or an ODEProblem, I think you suggesting I construct a DDEProblem with a Mass Matrix formulation and that’ll just work?
It seems maybe the Mass Matrix will get rather large quickly. If I have 1000 particles in 3D it’ll be something like 3000 coordinates, and then 1000x1000 time delay equations, so the mass matrix will be of the order of 1M x 1M entries? (though extremely sparse I guess)
Unfortunately that’s pretty fundamental to the issue. The idea is that particles feel the gravity of the distant particle at the position they “see” the particle with finite speed of light. Since the local and distant particles can accelerate chaotically, it’s impossible to eliminate, \tau_{ij} are all time dependent and this is a fundamental part of the problem.
Treating c as infinite gives you Newtons law basically.
But while I’ve got your attention… I’m considering doing a continuum version of this model to try to make it more tractable. It could be implemented a variety of ways, but conceptually the forces at a grid point come from summing the forces between the grid point and all other grid points. This could or could not be truncated at a certain distance, and if that works, it might make things much more tractable… But you now need to keep the history of the density field for some amount of time in the past.
I think I can turn this into a delay differential equation by using a Hilbert Space type approach (ie. represent the density field as a basis expansion, and work with the vector of basis coefficients, making the whole thing be a delay differential equation in those coordinates). Something like a Fourier or a Radial Basis Function expansion…
This approach doesn’t need to keep track of the individual Tau values, because you’d be doing integration over say grid points, and the delay between any two grid points is a fixed easily calculated value.
Is there an example you can point to? or which section of the docs I should read? Happy to give it a whirl but probably need a little more on-ramp.