Fundamentals of Iterative Prediction with Varying Force

It is rare to have a force which is perfectly constant, and iterative analysis of more realistic varying-force systems is substantially more complicated. A toy model demonstrates how programs may be written to analyze these systems.

Main Idea

The physics of iterative prediction with varying force is the same as for prediction with constant force, but it is necessary to generalize the mathematical expressions, which adds complexity to the code.

A Mathematical Model

To begin with, consider a one dimensional force, which may vary with both as a function of time and/or dependent variables such as position and velocity. Then we write this as [math]\displaystyle{ F(t,x,v) }[/math]. Now, using the momentum principle, we know that [math]\displaystyle{ F = \frac{\text{d}p}{\text{d}t} }[/math], which in discrete terms is [math]\displaystyle{ \Delta p = F\Delta t }[/math].

Just as with a constant force, this lets us write out for some iteration at [math]\displaystyle{ (t_0,x_0,v_0) }[/math] that

[math]\displaystyle{ p_{final} = p_{initial} + F(t_0,x_0,v_0)\Delta t }[/math]

Which we combine with kinematics to produce a new set of variables [math]\displaystyle{ (t_1,x_1,v_1) }[/math]. The difference we now have is that whereas before [math]\displaystyle{ F(t_0,x_0,v_0) = F(t_1,x_1,v_1) }[/math], we must now recalculate [math]\displaystyle{ F(t_1,x_1,v_1) }[/math] using the relevant formula. This will take the form of an extra step in each iteration. It is important to note that although we write the force as a function of all of these variables, in most cases it will only depend upon one of them. In a spring, for example, we will see that [math]\displaystyle{ F(t,x,v) = F(x) }[/math], meaning that only the position is necessary to compute the force.

A Computational Implementation

Code for an implementation of varying force iterative prediciton can be found in this Google colaboratory jupyter notebook. The following code snippet covers the physics calculation portions:

def calcmotion(self): #this is the method of the class which does the actual computation
   #First, we want to know how many steps we will calculate, and round this to an integer
   #Ideally, the user will make dur%tstep = 0, but we can't be sure, so we need to sanitize the input
   numsteps = int(self.dur/self.tstep)
   
   #Now we make an array to hold all of our data. For more information on numpy arrays, see 
   #https://docs.scipy.org/doc/numpy/user/basics.types.html
   #It will have 5 columns and numsteps number of rows
   #The columns are, in order, time, position, momentum, velocity, and force
   #np.zeros takes a tuple for the shape of the array, and the optional argument dtype for the data type
   #the plus 1 makes it so that we will end at the duration, instead of the duration - 1 timestep 
   self.data = np.zeros((numsteps+1,5),dtype=float)
   #Now we perform the actual work of iterating:
   for i in range(numsteps+1): 
     #Here we compute the time stamp of each step, but considering the iteration number and the timestep width
     self.data[i,0] = i*self.tstep
     #If this is the first step, we set our initial position and momentum
     if i ==0:
       self.data[i,1] = self.xinit
       self.data[i,2] = self.pinit
     #Otherwise, we use kinematics to compute a new position, and the momentum principle to compute a new momentum
     else:
       #xfinal = xinitial + v*delta t
       self.data[i,1] = self.data[i-1,1] + self.data[i-1,3]*self.tstep
       #pfinal = pinitial + F*delta t
       self.data[i,2] = self.data[i-1,2] + self.data[i-1,4]*self.tstep
     #We compute velocity using v = p/m
     self.data[i,3] = self.data[i,2]/self.mass  
       
     #Finally, we compute the force at our new position
     self.data[i,4] = self.forcefunc(self.data[i,0],self.data[i,1],self.data[i,2]) 
   return self.data

First, time data is computed by simply iterating over each time step, since we are using even steps. Next, if this is the first iteration, initial position and momentum data is input. If it is not the first iteration, then the momentum and position of the previous step are read to determine the new position, and the force and momentum of the previous step are read to determine the new force. Finally, the velocity column is filled out using the definition of momentum, and the new force is computed using the user defined force function.

The rest of the notebook consists of the code used to define these functions, and to view the results, both as plots and as animations.

Examples

Simple

Imagine a force on an object is defined by

[math]\displaystyle{ F=kx }[/math]

Where [math]\displaystyle{ k = 3 N/m }[/math]. Now, assume the object has an initial position of [math]\displaystyle{ x_0 = 5 m }[/math], and an initial momentum of [math]\displaystyle{ p = 0 }[/math]. After a time step of [math]\displaystyle{ 0.2 s }[/math], what will be the object's new momentum according to our model?

Solution

We have the formula [math]\displaystyle{ \Delta p = F \Delta t }[/math]. We may compute [math]\displaystyle{ F = (3 N/m)(5m) = 15 N }[/math], so plugging this in gives [math]\displaystyle{ \Delta p = (15 N)(0.2 s) = 3 N\cdot s = 3 \frac{kg\cdot m}{s} }[/math]

Moderately Difficult

Now, using the a varying force iterative motion script - either the one given above or your own - let's compute a motion prediction. Take the force defined above, with [math]\displaystyle{ k = 3 N/m }[/math] as above. Choose the initial position this time to be [math]\displaystyle{ 0.1 m }[/math], and let the initial momentum still be zero. Taking time steps of [math]\displaystyle{ 0.2 s }[/math] as above, what will be its position after 5 seconds? Include a plot of position versus time. What standard function does this plot resemble?

Solution

Plugging the problem into the above program produces a final position of [math]\displaystyle{ 84.8 m }[/math] (if this seems a bit small to you, make sure to do the difficult portion of the problem), and a plot that looks like this:

Error creating thumbnail: sh: /usr/bin/convert: No such file or directory Error code: 127

caption

This plot should be recognizable as an exponential function. In fact any force defined with the generic formula [math]\displaystyle{ F = kx }[/math] will produce an exponential position curve, so long as [math]\displaystyle{ k \gt 0 }[/math]. As it happens, the case [math]\displaystyle{ k \lt 0 }[/math] happens much more frequently in nature, and will be dealt with in the page (insert page on harmonic motion here).

Difficult

Using the same situation as above, let's now turn down the time step to [math]\displaystyle{ 0.01 s }[/math]. What is the final position now? Keep decreasing the time step until this converges to a final value. What is this value? If you can, try writing this in terms of fundamental constants and the parameters of the problem (hint: think about the discussion of the what the function looks like given above). If you haven't had any experience with differential equation this will be difficult, but the answer is very instructive. Think about why changing the time step caused such a radical difference, and give a brief explanation of your thoughts.

Solution

At a time step of [math]\displaystyle{ 0.01 s }[/math], the answer will come out to be [math]\displaystyle{ 267.84 m }[/math], which is dramatically larger than the result above. However, as we continue reducing the time step, we see that it converges to a value of [math]\displaystyle{ 288.45 m }[/math]. This may be expressed exactly as [math]\displaystyle{ 0.05 e^{\sqrt{3}t} }[/math]. To arrive at this value requires an understanding of differential equations, so if this is not familiar to you don't worry. One straightforward way of arriving at the solution is to use a language such as Mathematica/WolframAlpha. The below gives Mathematica code (if you are a Georgia Tech student, you have a freely available Mathematica license, and it can be a very useful program), in which one is able to specify the initial condition:

DSolve[{x"[t] == 3*x[t], x[0] == 0.1, x'[0] == 0}, x, t]

This is the simplified version required to make WolframAlpha run:

DSolve[x"[t] == 3*x[t], x,t]

This version will leave constants which you must solve. For a description of how this is done, see the solution below. Either way, you may notice that the result that these give (with constants filled in) is actually [math]\displaystyle{ 0.05 e^(-\sqrt{3} t) (1 + 1. e^(2 \sqrt{3} t)) }[/math]

This plot should be recognizable as an exponential function. In fact any force defined with the generic formula [math]\displaystyle{ F = kx }[/math] will produce an exponential position curve, so long as [math]\displaystyle{ k \gt 0 }[/math]. As it happens, the case [math]\displaystyle{ k \lt 0 }[/math] happens much more frequently in nature, and will be dealt with in the page (insert page on harmonic motion here).