Simulation #2: Ball motion in gravitational field (part 2 - program)

Previous simulation #1 was visualisation rather than real simulation, because it does not show for user any calculated physical properties. This simulation is different, because:

• it calculates physical properties that describe physical phenomenon (ball motion in gravitational field) with units;
• it compares results with theory, to assess calculation error.

Numerical simulations are used for solving complex equations like differential equations for which finding analytical solutions is very difficult or the solution even does not exists. Ball motion in gravitational field, on the other hand is described by well known simple physics formulas. Therefore numerical methods has no practical usage in that case, because even simple calculator can do the work. However, I introduce numerical methods for this case for educational purpose. I am going to show that we can calculate different physical properties describing ball's motion without using final formulas.

In details, the purpose of program that realizes this simulation is:

• to make possible for user to enter easily initial values of physical properties that describe ball in $t=0$ moment (initial velocity $\vec{V}_0=[V_{x0}, V_{y0}]$ of ball and its initial height $h_0$);
• to make possible for user to enter size of time step $dt$ used in simulation;
• to calculate physical properties that describe ball motion with units (max height $h_{max}$, final velocity $V_{end}$, max distance $d$, time of flight $t_{end}$, time to reach max height $t_{max}$);
• to visualize ball motion (without alignment of time in simulation with clock time, but with alignment with calculation progress);
• to calculate delta between simulated results and theoretical results (based on formulas from part 1);
• to make possible for user to store results from at least one simulation run to be able to compare them with results from other run.

Program/simulation presentation

Below picture (Fig 1) presents UI of the program. It contains 3 radio boxes to select 1 from 3 motions. Each one has related set of input fields in which we can set entry conditions.

Fig 1.

When simulation is finished, results values that describe ball motion are calculated (which is indicated by red outline around resize button). Results can be saved in the second table. It is very useful if we want to compare results between two different simulations. Fig. 2 shows presentation of calculated results:

Fig 2.

Below picture (Fig 3) shows how time step $dt$ size increases/decrease error of calculations. When the time step decreases, the accuracy of calculations increases and vice versa. Using second table we can compare results easily.

Fig 3.

It can be seen that ball does not hit the ground on above pictures. It is because, simulation calculates and visualizes the last step when the ball height is still above the ground: $y\geq 0$. Case $y=0$ for last step happens only for enough small time step value $dt$ or accidentally for larger $dt$ (in lucky cases). The smaller $dt$ is - the closer ball to the ground will be. $dt$ with order of magnitude $10^{-4}$ is enough to reach $y=0$ for all types of motions (all entry conditions). Below picture (Fig 4.) shows results of calculation with $dt=0.0001$:

Fig 4.

It can be seen that calculations are very accurate. Difference between theoretical and simulated values are exactly the same with the 2nd decimal place. Excellent view!

The following figure (Figure 5) shows the possibility of adjusting the angle and distance of the camera setting from the object on the scene in VPython. Change can be done by hold and drag right or both buttons of the mouse:

Fig 5.

It looks interesting and can be useful for example for more visually complex simulations (with many moving objects) like for example model of solar system.

Implementation

To write this program, I used combination of wxPython and VPython libraries. wxPython is library for making UI in Python and VPython is scientific library (for Physics mainly) that display 3D objects and animations. Objects in that meaning are for example: vectors, spheres, cubes.. They can look like 2D, if we do not move them along z-axis, however in fact they are always 3D (it can be seen on Fig. 5). I do not use tkinter library for UI like in Simulation#1, because VPython uses wxPython. VPython usually opens animation in new window and I want to present everything in only one window, so to be able to do this we need to use wxPython.

I) Program structure:

Program consists of two classes: Ball and Simulation. Beyond that, it contains couple of functions whose are related to UI. The following tables contains list of variable and methods with short description about each one:

Class: Ball
Variables	Description
t	time
dt	time step
g	gravitational acceleration
x0,y0	initial coordinates of the ball (for t=0)
vx0,vy0	initial velocity of the ball (for t=0)
x,y	current ball coordinates
vx,vy	current ball velocity
nextY	auxiliary variable, it stores value of y for next calculation's interation
mysphere	Vpython object (sphere) which represents the ball
velocityVector	Vpython object (arrow) which represents velocity vector of the ball
Methods	Description
nextStep	calculates next ball position and velocity after time t=dt

Class: Simulation
Variables	Description
g	gravitational acceleration
ifSimulationWorks	auxiliary variable of boolean type, is True if simulation is already working
simulationCalculatedParameters	dictionary, where keys are physical properties that describe ball motion (hmax, Vend, d, hmin, tend, vxend, vyend) and values are their results
initialData	dictionary, stores: initial position, initial velocity and time step values
Methods	Description
runSimulation	runs simulation
stopSimulation	stops simulation
initEntryData	auxiliary method, reads data from all input fields and store them in initialData dictionary
calculateEndParameters	method calculates physical properties that describe ball motion: hmax, Vend, d, hmin, tmax, tend, vxend, vyend
getTheoreticalResults	method calculates theoretical values for physical properties that describe ball motion: hmax, Vend, d, tmax, tend

Table related to class Simulation does not contain all variables and methods used in implementation. Some of them are UI related and I skipped them.

You can find project files on github: https://github.com/sim-num/BallInGravityField.


II) Steps to program simulation:

1. Set initial ball's position and velocity
2. Calculate position and velocity for next simulation step
3. If calculated position is above the ground accept results, if not - stop simulation
4. Update calculated ball's physical properties
5. Go to step 2

III) Implementation of the steps:

• Ad 1. (Set initial ball's position and velocity)
  
  Initial position and velocity is set by user in proper fields in UI. Then, those values are read by method initEntryData from class Simulation and stored in dictionary initialData:

  
      def initEntryData(self):
          if radio1.GetValue():  #if radio button nr 1 is selected (is true)
              self.initialData={'x0':0,'y0':int(inputH0_1.GetValue()), 'vx0':int(inputVx0_1.GetValue()), 'vy0':int(inputVy0_1.GetValue()), 'dt':float(inputDt.GetValue())}
          if radio2.GetValue():
              self.initialData={'x0':0,'y0':int(inputH0_2.GetValue()), 'vx0':int(inputVx0_2.GetValue()), 'vy0':int(inputVy0_2.GetValue()), 'dt':float(inputDt.GetValue())}
          if radio3.GetValue():
              self.initialData={'x0':0,'y0':int(inputH0_3.GetValue()), 'vx0':int(inputVx0_3.GetValue()), 'vy0':int(inputVy0_3.GetValue()), 'dt':float(inputDt.GetValue())}
  

initEntryData method is called at the beginning of simulation inside runSimulation method from class Simulation (see code in Ad 4.). Then initialData dict is passed to contructor of class Ball to initialize new Ball instance:

class Ball(object):
    ...
    def __init__(self, initialData, dtt):
        self.vx0 = initialData['vx0']
        self.vy0 = initialData['vy0']
        self.x0 = initialData['x0']
        self.y0 = initialData['y0']

NOTICE: "$\ldots$" in code represents some skipped part of the code irrelevant in described context.


• Ad 2. (Calculate position and velocity for next simulation step)
  
  This step is implemented by method nextStep from class Ball:

  
      def nextStep(self):
          self.nextY = self.y + self.vy * self.dt
          if self.nextY > 0:
           self.t+=self.dt
           self.vy = self.vy-self.g*self.dt
           self.x = self.x + self.vx * self.dt
           self.y =self.y + self.vy * self.dt
           ...
  
  

  $V_x$ component of ball's velocity vector is constant, therefore it is not calculated during simulation. What we need to calculate are components of position vector: $x$, $y$ and velocity component $V_y$. They can be denoted by $x(t)$, $y(t)$, $V_y(t)$, because they are time dependent. All of them are calculated numerically for each time step. It is done using above algorithm. This part of algorithm is so important that requires another part of this topic for exmplanation.
  
• Ad 3. (If calculated position is above the ground accept results, if not - stop simulation)
  
  In method nextStep, we check if after calculation of next step, ball position is still above the ground. If yes - then method returns True, of not - then returns False, which is the indicator to stop the simulation.

    def nextStep(self):
        self.nextY = self.y + self.vy * self.dt
        if self.nextY > 0:
         self.t+=self.dt
         self.vy = self.vy-self.g * self.dt
         self.x = self.x + self.vx * self.dt
         self.y = self.y + self.vy * self.dt
         ...
         return True
        else:
         return False

• Ad 4. (Update calculated ball's physical properties)
  
  nextStep method is called inside method runSimulation from class Simulation. After each simulation step, for which ball is still above the ground, as I mentioned previously - nextStep method returns value True. Next, this value is assigned to variable status. Its value is condition for while loop in which we invoke nextStep method and calculateEndParameters method. This last method is responsible for calculation of physical properties that are stored in simulationCalculatedParameters dictionary. The code below shows hierarchy of methods and while loop:

  
      def runSimulation(self, event=None):
          ...
          self.simulationCalculatedParameters={'hmax':-1,'vxend':-1,'vyend':-1,'Vend':-1,'d':-1,'tend':-1,'hmin':1000000}
          ...
          self.initEntryData()
          self.ourBall=gravityFieldBall.Ball(self.initialData, float(inputDt.GetValue()))
  
  
          if self.ifSimulationWorks == False:
              self.ifSimulationWorks = True
              status=True
              ...
              i=0
              ...
              while self.ifSimulationWorks and status:
                  self.calculateEndParameters()
                  ...
                  status=self.ourBall.nextStep()
                  ...
                  i += 1
          ...
          self.ifSimulationWorks=False
          return
  
  

  Code of calculateEndParameters method is presented below:

  
      def calculateEndParameters(self):
          if self.ourBall.y>self.simulationCalculatedParameters['hmax']:
              self.simulationCalculatedParameters['hmax']=self.ourBall.y
              self.simulationCalculatedParameters['tmax']=self.ourBall.t
  
          if self.ourBall.y>0 and self.ourBall.y<self.simulationCalculatedParameters['hmin'] and self.ourBall.vy<0:
              self.simulationCalculatedParameters['hmin']=self.ourBall.y
              self.simulationCalculatedParameters['d']=self.ourBall.x
              self.simulationCalculatedParameters['vxend']=self.ourBall.vx
              self.simulationCalculatedParameters['vyend']=self.ourBall.vy
              self.simulationCalculatedParameters['Vend']=math.sqrt(self.ourBall.vx**2+self.ourBall.vy**2)
              self.simulationCalculatedParameters['tend']=self.ourBall.t
  

  First if condition checks if current y is the largest value found, if yes - then it is stored in simulationCalculatedParameters dict as hmax value with corresponding tmax value to it. The second if condition checks if ball is still above the ground, if current y position component is at the same time the smallest value found . The last part of 2nd condition checks if ball is falling down (Vy velocity component is negative). We need to consider only falling part of the ball motion for searching physical properties like: $V_{end}$, $t_{end}$, $d$. hmin is auxiliary variable, used to figure out if current y position is the lowest position of the ball above the ground (during falling down part).
• Ad 5. (Go to step 2)
  
  This step is realized by while loop inside runSimulation method.

IV) Simulation visualisation

This time I used VPython library for making simulation's visualisation. Because I used it for the first time here, I will list every object I used with quick description:

• display - to combine VPython window with wxPython window, we define display object. It creates VPython scene. To create it, we put window type object from wxPython libary as constructor parameter:

  
  w = window(width=L, height=W,
             menus=False, title='Ball motion in gravity field - 2D simulation',
             style=wx.SYSTEM_MENU | wx.CAPTION | wx.CLOSE_BOX)
  
  #create scene in VPython inside window:
  scene = display(window=w, x=10, y=105, width=600, height=500, center=(40,28,0), background=(0,0,0), autoscale=True)
  

• label - object is used for creating label (for example white box with "h= V=" values):

  
  hV_label = label(pos=(45, 55, 0), text='h=\nV=', xoffset=1, line=0, box=True, opacity=0)
  

  Another example is axis description:

  
  labelX = text(text='x [m]', depth=0.4, color=color.white, height=1.5, pos=[ymax, -3 * tic_h], font='serif')
  

• arrow - object is used for creating axes in XY coordinate system:

  
  xmin = 0.
  xmax = 50.
  ymin = 0.
  ymax = 55.
  
  # tick marks
  tic_dx = 5
  tic_h = 1
  self.xaxis = arrow(pos=(xmin, 0, 0), axis=(xmax+tic_dx, 0, 0), shaftwidth=0.2)        # axes
  
  self.yaxis = arrow(pos=(0, ymin, 0), axis=(0, ymax+tic_dx, 0), shaftwidth=0.2)
  

  The same object can be used to create $\vec{g}$ vector:

  
  self.gravityVector = arrow(pos=(70, 55, 0), axis=(0, -5, 0), shaftwidth=0.5, color=color.blue)
  

• curve - I use it to create markers on the axes:

  
  for i in arange(xmin,xmax+tic_dx,tic_dx):
   tic = curve(pos=[(i,-0.5*tic_h),(i,0.5*tic_h)])
  

• sphere - object is used for creating ball:

  
  self.mysphere = sphere(pos=vector(self.x0, self.y0, 0), radius=self.R, color=color.red)
  

Each object has some attributes like for example visible. I use visible attribute to show or hide vector or sphere:

         if showVelcheck.GetValue()==True:
            self.ourBall.velocityVector.visible=True
         else:
            self.ourBall.velocityVector.visible=False

Summary

Simulation is done according to the plan. It shows that simple formulas can be calculated in other way - numerically. Simulation can be extended by adding wind factor or air resistance. I will add them in next parts. Algorithm used for calculations can be also improved. Its description and improvement methods will be shown in part 3. During this simulation I learnt that VPython is powerful framework for making visualisation. Therefore, I am going to use it for the future simulations.