If we are to use python to write powerful mathematical programs, we need those programs to be able to make decisions and perform repetitive operations. All programming languages can do those things. The commands that do that are called flow control statements. In particular, every programming language possesses at least two statements that allow us to control the flow of execution. Decisions are made using "if" statements, while repetitive operations can be done using "for" statements.
The "if" statement in python has the following structure.
if condition1: statements to execute elif condition2: more statements else: default statementsRemember that in python, indentations are important! Only the first two lines are really required. The other lines can be added if more options are needed. The execution would go like this: python looks at condition1 (e.g. x>0) and if it is true, it performs the statements to execute. If condition1 is false, then python checks condition2 (e.g. x>-1), and if that is true, it executes the more statements. If condition1 and condition2 are both false, python executes default statements. In this construction, "if", "elif", and "else" are keywords, while the rest is content you will create yourself. Note that "elif" stands for "else if", i.e., if the previous condition was not true, try this. Let's look at an example. We want to write a function that will perform the composite trapezoidal rule to approximate the integral of a function designated here by fct, i.e. it approximates
def trapezoid(fct,a,b,n): if n<=0: print "Error: n must be positive" return False
We first define the function, called "trapezoid". It takes four arguments: fct, a, b, and n. The last number n denotes the number of subintervals used for the composite trapezoidal rule. It seems important that the number n be positive. Otherwise the scheme will fail. Thus, we test n to see if it is less than or equal to zero. If so, the the program does the indented code, which prints an error message and quits the function. When we run this code we get the following results.
>>> import pythonfcts as mf >>> def f(x): ... return x*x ... >>> mf.trapezoid(f,0,1,-4) Error: n must be positive FalsePython understands several basic comparison operators, as detailed in Table 1. It can also combine these using logical "and" and "or" operators. For example, we could have written the "if" statement equivalently as if n<0 or n==0.
|==||...is equal to|
|!=||...is not equal to|
|<||...is less than|
|>||...is greater than|
|<=||...is less than or equal to|
|>=||...is greater than or equal to|
The "for" keyword allows us to perform some set of commands repetitively, with a counting variable that changes in some prescribed way. The basic syntax looks like this.
for i in list: some commands that might depend on iIn this, the word "for" must appear exactly as shown. The variable that follows is denoted "i" here, but could have any name. The keyword "in" must come next, and then some list (frequently created as a range or arange) from which the values of i appears. The colon at the end of the line is not optional. After that we use some uniform indentation to put as many commands as we want to be repeated. Here is an example.
>>> for i in [-1,1.5,7]: ... j = i**2 ... print j ... 1 2.25 49
By the way, the ** operator raises the argument to its left to the power given on its right - it is an exponentiation operator. Thus, i**2 = pow(i,j). In these commands, our counting variable is called i, and we tell it to take on the values -1, 1.5, and 7, in turn. Each time it takes one of those values, we compute a variable j as the square of i, and then print j. We could do the list differently, using an arange, viz.
>>> from numpy import * >>> for i in arange(-1,3): ... j = i**2 ... print j ... 1 0 1 4
This time we need numpy. The point is that "for" is a plain python construct - we only needed numpy when we wanted to use an arange.
In our trapezoidal rule example, recall that the trapezoidal rule can be written as when and . We can type this into our function file to make our trapezoid function look like the following.
def trapezoid(fct,a,b,n): if n<0 or n==0: print "Error: n must be positive" return False h = (double(b)-double(a))/double(n) total = 0.0 for i in arange(n+1): total += fct(a+i*h) total *= 2.0 total -= fct(a)+fct(b) return total*h/2.0
Have we ever mentioned that indentations are very important in python? The indentations here indicate the code that runs if the condition is true, and the code that runs repeatedly in the "for" construction. After the "if" part we discussed earlier, we initialize a variable called total to zero. This is where we accumulate the sum from the equation. We then use a "for" construction to run through indices i starting with 0 and ending with n, using these to construct points at which to evaluate the functions. We multiply these function values by two as in the equation, and then subtract off the numbers on the ends. Notice the odd notation total += 2.0*fct(a+i*h). That += means "is equal to itself plus the stuff on the right." Or again, we could say that as "compute the quantity total+2.0*fct(a+i*h) and then put that new value in the variable total. This kind of notation has become common in most modern programming languages - thank you Dennis Ritchie.
The next line is not inside the "for" construction, because it is not indented. Thus, it is only executed once. It takes the total and multiplies it by the 2 from the equation. We could have done this inside the "for" construction, but then that would have been n multiplications. This way we have only one multiplication, and we get to illustrate the use of the *= operator. Finally we just subtract fct(a) and fct(b) from the total. Finally, we multiply the whole total by the quantity outside the brackets in the equation, and send it back to the python session as the answer.
Now we can run this function in our python session.
>>> reload(mf) <module 'myfunctions' from 'myfunctions.pyc'> >>> mf.trapezoid(f,0,1,100) 0.33335000000000004 >>> mf.trapezoid(f,0,1,1000) 0.33333349999999995
We used the function f that we defined earlier. The correct answer would be 1/3, so we see that our approximation is correct to four decimal digits with 100 subintervals, and six digits with 1000 subintervals. It seems like our function works.