Gravitational Force Interaction in MATLAB 

Introduction

The Gravitational Force is one of the oldest forces that has been discovered by Physicists and is still baffling on some levels. Sir Isaac Newton was the first to define the force mathematically. He presented an analytical relation between the magnitude of the force and the masses interacting and the distance of separation between them. In his first law of motion, Newton said that force is the amount of change in momentum producing it. (i.e., a force applied to an object has an equal and opposite reaction) Hence, Gravitational force results from mass and momentum interaction according to gravitation equations as stated below.

Originally, gravitational force was defined by Newton as a result of two masses. It was discovered that such a situation typically did not occur in nature. Still, the realization didn't stop scientists from exploring further into this force which would eventually become known as gravitational interaction. The main problem with this theory is its inability to predict experimental results accurately when another body was brought into play. Eventually, it was discovered that an increase in the number of objects (m1, m2) involved in the equation would prove to be more accurate.

So far, there hasn't been any experimentation conducted to test the validity of this theory or any conservation of energy equations in general. How can one be able to confirm a conservation equation if they haven't conducted any experiments? The closest thing to an experiment would be observations. We can confirm that gravitational interaction exists between two bodies, which means we have confirmation about the conservation of energy for sure. However, keeping these problems for the core Physicists, we will just look upon the basic form of the equation and use it to describe motions using MATLAB.

Simulating Gravitational Force in free space

The Gravitation Equation is an approximation that has been derived based on Newton's Law of Universal Gravitation. In Physics, Newton's Law states that each point mass attracts every other point mass by a directly proportional force to the product of their masses and inversely proportional to the square of the distance between them. The proportionality constant, the gravitational constant, in this case, has a very precise value of devised by very sophisticated experimentations.

While two bodies in free space is a very trivial and straightforward model of Gravitational force interaction, it is fundamental in describing every other complex motion. So, we will begin our simulation journey by having a mass of 10kg fixed in spatial dimensions at a point and a free body interacting in its Gravitational field. A free-body model is one where no external forces act on it andare simply subject to its mass and gravity.

We have our object of mass 10kg and an empty space where another object of unspecified mass will be our free body. Since this is a Newtonian simulation, with every tick of time, the force applied to the mass will be equal to the Gravitational constant * mass of Earth * distance between them. Since we are not computing acceleration but using the time taken to reach an equilibrium state, we will need a function that calculates how long it takes a given mass at rest to come to rest (i.e., distance reduction).

While the formula is descriptive enough for calculations, it cannot be used directly for simulating the scenario. The formula must be converted into an ordinary differential equation (ODE) format describing the conditions or forces it depends on. So, the ODE analog of the equation is as shown below.

where g also depends upon the masses concerned and the distance between the objects, as already mentioned before. Here is the respective MATLAB code for using the ODE for simulating the scenario described above.

Runge-Kutta Method

The solver technique you just saw in the MATLAB code must have raised many questions in your mind, like what exactly is the ode45 function, which was able to solve the equation. But for that, we need to go to the basics for a while.

Questions concerning ordinary differential equations (ODEs) may always be simplified to the study of differential equations of first order. For example, the following second-order equation

can be rewritten as two first-order differential equations

where z is a new variable. Following the same trend, we can reformulate the former differential equation concerning the gravitational force into two subsequent ODEs.

The ode45 function calculates the next point in the resultant sequence with the var_r function linked to it, presenting the two ODEs after calculating the variable distance from the fixed point. Let us see how the function operates internally.

Runge-Kutta methods propagate a solution over an interval by combining the information from several Euler-style steps (each involving one evaluation of the right-hand function) and then using the information obtained to match a Taylor series expansion up to some higher order. The formula for the Euler method is

which advances a solution from to .

The formula is unsymmetrical: It advances through an interval h but uses derivative information only at the beginning of that interval.

There are many versions of the Runge Kutta method. Still, the most used version is the classical fourth-order Runge-Kutta formula, as shown below:

The fourth-order Runge-Kutta method requires four evaluations of the right-hand side per step as shown in the following figure.

The Runge-Kutta method treats every step in a sequence of steps identically. Prior behavior of a solution is not used in its propagation. This is mathematically proper since any point along the trajectory of an ordinary differential equation can serve as an initial point. All steps are treated identically and make it easy to incorporate Runge-Kutta into relatively simple "driver" schemes. The function 'ode45' is the same implementation of the described algorithm calculating the points in numerical order for our purpose as well.

Simulating a free fall

The former case was a typical example of a second-order ordinary differential equation into action. The same procedure can be implemented for any other ODE and replace the initial conditions on the way. Now, we move to a procedure solving symbolic ODEs using the 'dsolve' function in MATLAB. Simulating a freefall on Earth is not as complex as the previous one. In this scenario, the change in the value of g (acceleration due to gravity can be taken as a constant throughout the process, simplifying the process. With the same ODE in hand, we simulate a condition where a ball is set to freefall from a height of 100m. Upon each impact, the ball loses 20% of the impact velocity due to inelastic collision. We try to render the animation till the max height reduces to only 5m. Let's look at the MATLAB code now.

This case didn't use the numerical ODE solver and instead uses the symbolic solver where only two variables are considered. The height is progressing over time, and that's what the ODE solver tries to describe with the help of an equation.

Conclusion & Discussion

As Physicists boast, Mathematics is the language of the Universe, and if they are correct, solving differential equations is a powerful tool for understanding the Universe. We went through two different approaches while solving the ordinary differential equation. Each has its own domain of operations application. It is now the person's discretion to use it to find which one is more suitable for the specific purpose.