For Loops, If Statements and Euler's Method

Introduction

This week, we will introduce a couple of programming structures: for loops and if statements. Loops allow sections of code to be run many times, iterating over a certain variable in a loop. This is called iteration, and is a very useful part of programming to solve physical problems that evolve over time. If statements allow you to test whether a condition is satisfied and proceed accordingly.

We will use for loops to integrate differential equations numerically. We will start with an equation that you already know how to solve, before moving onto some more complex equations. We will make use of Euler's method. It is the most basic method for solving differential equations.

Useful Links

For Loops - The basics of for loops in Python (This page talks about lists because, as previously mentioned, arrays are a feature of NumPy rather than Python itself.)

If statements - The basics of if statements in Python.

Looping Techniques - More detailed examples of using for loops in different situations, very useful.

Syntax Summary

for x in values:
    print(x)
    # indented loop code

if x > 3:
    print("x is greater than 3")
else:
    print("x is less than (or equal to) 3")

Worksheet

For Loops

The for loop is a construct that appears in virtually all programming languages. It allows you to loop over a set of values repeatedly performing a specified calculation or operation.

For Loops

In Python, a for loop iterates over all of the elements of a list, string or array. The standard structure of a for loop is:

values = array([7, 6, 42])
total = 0
for val in values:
    print(f'Current value is {val}')
    total += val

print(f'The total sum is {total}')

[Line 1] is creating an array of items that the loop will iterate over. The loop will be run once for each item in this list. You can use lists here, or lists of words, letters or numbers.

[Line 2] sets up a variable total with value zero. This is not a necessary part of for loops, it is included here to help illustrate the process.

[Line 3] starts the loop. It takes the first item in the list values and calls it val, then runs the loop, then takes the second item and calls it val, until there are no more items.

[Line 4] is inside the loop. Any indented lines are in the loop. This will print out whatever item is currently called "val".

[Line 5] is also inside the loop. The code

total += val

does the same as

total = total + val

but saves you from repeating total. This syntax is quite commonly used, as it makes code shorter and easier to read, but it can take some getting used to.

[Line 6] is blank

[Line 7] is outside the loop because it is not indented. This line will only be run once the loop has gone over each item in the list values, it prints out the total.

We will adopt the convention that every indent will be 4 spaces wide. The Jupyter editor will do this automatically for you: if you type for val in values: and hit return, the cursor will move to the next line and be indented 4 spaces.

This program should the values 7,6,42 and fianlly 55. This is printing the values from the list in each loop, then printing the final value of total which is 55.

In this case, you can achieve the same thing in one line by simply doing

total = sum(values)

The point here is to illustrate the basic idea of a for loop.

Note that you could make values anything (they don't need to be numbers) for example days of the week:

days = ["Monday", "Tuesday", "Wednesday"]
# loop over values
for day in days:
    print(f"today is {day}" )

Example: Fibonacci Numbers

The Fibonacci numbers are a series of numbers that start off with:

0, 1, 1, 2, 3, 5, 8, ...

The rule for generating the next number is that it must be equal to the sum of the previous two numbers. Algebraically, that is

Fn = Fn-1 + Fn-2

where the subscript n means the n^th Fibonacci number. We must also specify F0 = 0, F1 = 1.

It is straightforward to write a code to generate the first 50 Fibonacci numbers, using a for loop to calculate each subsequent number based on the previous values.

# set up empty array of 50 zeros
n_fib = 50
fib = zeros(n_fib, dtype=int)
# set values of first 2 numbers in the "fib" array
fib[0] = 0
fib[1] = 1
# loop to calculate the remaining values
for i in range(2,n_fib):
    fib[i] = fib[i-1] + fib[i-2]

print(fib)

If Statements

It is often useful to be able to test the validity of a condition and then have the code proceed accordingly. We have already seen the syntax for testing whether a statement is true, for example

a = 4
b = 5
print(a > b)
print(a < b)

Rather than just printing out the result, we can use it to decide what to do next, for example:

n = int( input("Please enter a number\n") )

if (n > 0):
    print("You entered a positive number")
elif (n < 0):
    print("You entered a negative number")
else:
    print("You entered zero")

This illustrates the basic layout of an if statement. First we test whether a given condition is true (n > 0) and if it is, proceed accordingly. You can test a second condition using elif (a contraction of else if) and do something different. The else statement is the path that's taken if none of the above are true.

Note: The indentation is the same as in for loops. It is only by indenting that Python knows where the if/elif/else blocks of code start and end.

Breaking out of loops

Sometimes, we will want to use an if statement within a loop and, if the statement is true, stop the loop. For example, if we have a code to calculate the trajectory of a projectile, it would be very natural to stop the evolution once the height became zero (the projectile hits the ground). This is achieved using the break command. For example

numbers = arange(20)

for num in numbers:
    if num > 10:
        print(f"{num} is bigger than 10")
        break

print(num) 

If you run this code, it will tell you that 11 is greater than 10 and then exit the loop. So, the last line will print 11 (rather than 20 as it would if it got to the end of the loop).

Euler's Method

To briefly describe Euler's Method, if we know that dy/dx = f(x,y) and have an initial values x0 and y0, we can find the value of y at x1 = x0 + h as

y(x₁) = y(x₀) + h f(x₀, y₀)

Once the value is known at x₁ we can use the same procedure to calculate it at the next point. This method becomes much more accurate if smaller steps are used, and this can be easily accomplished using for loops in Python.

Step by Step Method

Before we automate the procedure using for loops, let us work out a short example showing what each step of a for loop will do. To simplify the example, we'll use freefall (no air resistance). F = mg, so using Euler's method with a time step of dt, we have v₁ = v₀ + g dt:

# Manually calculate 2 time steps
t = zeros(3)
v = zeros(3)
dt = 1 # time step
g = - 9.81

t[0] = 0 # Assume initial time is zero 
v[0] = 0 # Assume initial velocity is zero
t[1] = t[0] + dt
v[1] = v[0] + dt * g
t[2] = t[1] + dt
v[2] = v[1] + dt * g

This method is OK if you only want to calculate a couple of points, but after that becomes very tedious. The simple solution to this is to create an iterative program which can run the same line as many times as you want, using for loops.

Iterative Method

First, make an array of time for the loop to run over, and an array of zeros to store your velocities. Don't forget to define all variables used in you loop, in this case g:

t_max = 10 # maximum time
dt = 0.1 # time step
g = -9.81
n_steps = int(t_max/ dt)
t = zeros(n_steps)
v = zeros(n_steps)

# set initial conditions
t[0] = 0
v[0] = 0

for i in range(1,n_steps):
    t[i] = t[i-1] + dt
    v[i] = v[i-1] + dt * g

If you plot v against t, you will get a straight line with a slope of 9.81ms^-2. Of course, we didn't really need a for loop to solve this problem.

Integrating to find position

We have calculated an approximation of the velocity at discrete points. It is easy to add one extra line of code to work out the position as a function of time. Using Euler's method, we can approximate

x[i] = x[i-1] +  dt*v[i-1]

So, the loop at the end of the code becomes:

for i in range(1,n_steps):
    t[i] = t[i-1] + dt
    v[i] = v[i-1] + dt * g
    x[i] = x[i-1] + dt * v[i-1]

Free Fall with Air Resistance

Let's move onto a problem that is more difficult to solve with pen and paper. We will include the effect of air resistance into freefall. Air resistance is well modelled by F_air = - kv². It acts against the direction of motion, and the force is proportional to v². The value of the constant k will depend upon the area of the object and its drag coefficient.

You should now be in a position to write a code to calculate the velocity as a function of time for an object in freefall with air resistance. You can simply alter your code for the system without air resistance to add in the extra terms. Remember that Euler's method gives:

v[i] = v[i-1] +   dt * F(v[i-1])/ m  # F/m is acceleration

Excercises