Flying A Rocket, Numerically.
A brief introduction to simulations and converting discretized Ordinary Differential Equations into a python code to model a physical phenomenon.
The Universe is continuous or discrete is up for a debate (Greek philosopher Zeno’s glance over two giants: General theory of relativity & Quantum theory). Nonetheless, a physical phenomenon such as flow of an inviscid fluid even though seems continuous from a distance, we still need to discretize them in order to quantitatively ‘measure’ any of its property. And when we, humans ( we will keep an Earthling’s perspective in this article ), measure something we always lose some accuracy in the process. Isn’t it surprising? We shall get back to this a little later.
Studying a physical phenomena with the help of conventional physical experiments can be extremely difficult (High energy plasma and electric discharges) and sometimes impossible (Black holes). Simulations enable us to build a model of a system and allows us to do virtual experiments.
So, what does this very process of discretization, actually mean? one may ask. Before we go further, we need to put the whole narrative in a more technical language to better understand the underlying principle. The continuous flow, we have seen earlier is described by a continuous function, say, velocity of a particle is converted into a discrete form that can be represented in the finite memory of the computer.
The next and most important question is how do we do that ? This is where numerical methods come into action. Surely, there are other methods (Spectral Element Element, Lattice Boltzmann Method) as well but for the sake of brevity, we will stick to the simplest numerical method to solve ordinary differential equations (ODEs). So, let’s begin out journey through this beautiful world of simulations where we will get to know a little bit about numerical integration method and the process of discretization. Then, we will simulate a rocket flight plotting its purely vertical trajectory to observe the power of simulation with a few lines of python code.
Rocket Science
The motion of a rocket in a purely vertical flight is governed by a set of first-order differential equations, as given below.
Where v is the velocity of the rocket, h is the altitude of the rocket.
The mass of the remaining propellant is given by:
One might have noticed that since our third equation is also time-dependent, we actually have three equations to solve now. We have three variables to solve for, i.e velocity (v), height (h) and mass of the propellant (Mp).
Notice that the variables are present in other equations as well. This is called coupling of equations, since the variables are also dependent on other variables at each time step.
Time Integration
We want to find out the trajectory of the rocket. This can be done by integrating the differential equations with time. The temporal solution for these differential equations can be obtained in two steps:
- Finding the derivatives of all our variables at a point in time.
- Integrating these solutions with time.
Note that there exists an analytical solution for these set of differential equations . However, since we are not verifying the accuracy of our numerical solution, we will just focus on the numerical analysis aspect of the equations.
Let’s rewrite the time derivative with prime.
Another way to look at these first-order differential equations is to use vector notation. We can make a vector with our three variables.
Then differential system can be written as
If the Right-hand side is termed as f(u), we can write the full vector system as u′(t)=f(u). Notice, that we have dropped the vector sign just for the sake of clarity. We shall, henceforth, use this notation often to describe the full vector system of equations.
We can define a function in python to write this right-hand side.
Initial Value Problem
Since, we move in temporal dimension we need a starting point which is at u(t=0)=u_0. We need to have certain initial values for all the variables in our equations. This is called Initial Value Problem.
Discretization
We have used differential equations with a time derivative. But what exactly is a derivative ? For the curve of u = u(t), a derivative represents the slope of the tangent at a point on this curve. The definition of the derivative u′ is given as
If the step Δt is very small, we can approximate the derivative by dropping the limit,
What does this approximation implies ? It means we lose a certain accuracy while determining the actual value. Rings a bell ? I talked about losing accuracy in the very beginning ! This is one of those instances where we lose accuracy depending on the order of the derivative. Here, we are using first-order derivative.
What does the above equation really means? We can find the new value of u(t+Δt) by using an initial value of u at time t and u′ which we already know and stepping by Δt. Similarly, we can find u(t+2Δt) and so on. This is actually the second step i.e integration of solutions with time.
We have actually discretized the time domain into small time steps Δt. We find the numerical solution of variables for a range of values of time each separated by Δt. These time values form our time-grid. The first point of this time-grid is the initial time t=0.
Euler’s Method
Based on the approximate equation to determine u(t+Δt) above, the numerical solution of the equation consists of computing a sequence of approximate solutions using the following equation:
Numerical Solution
Now that we have all theory we need, we can solve these differential equations using python and its efficient libraries.
Firstly, we shall import some important libraries.
We define a certain time interval which can be arbitrary. We create an array of u for holding a vector of all variables. For the initial step, we fill the u with initial values of all variables defined above.
Next, we plot the trajectory of the rocket which is u[:,0] as a function of t.
As we can observe, the rocket goes up to a maximum height of around 1300 m in 15.7 seconds. At this moment, the fuel inside the rocket is empty and the rocket starts to descend.
Conclusions
This was a brief introduction to simulating a physical phenomenon using python which can be done on a laptop. The beauty of simulations lies in the approximations methods with which we can study a lot of complex phenomena. However, as we saw above, we lose some accuracy when we approximate the derivatives.
We can further calculate the rocket’s velocity and the time taken by the rocket to hit ground after it runs out of fuel. Lastly, we can compare the results with the analytical solution to check the accuracy of our numerical solution. The complete code can be found in the link here.
