$$ \newcommand{\tp}{\thinspace .} $$

Ch.1: Computing with formulas

Hans Petter Langtangen [1, 2]

[1] Simula Research Laboratory
[2] University of Oslo, Dept. of Informatics

Aug 21, 2016










Why program?

Everybody in this country should learn how to program a computer... because it teaches you how to think. Steve Jobs, 1955-2011.









The teaching strategy is example-based









The learning strategy is about doing exercises









Chapter 1 is about evaluating formulas

Why?









Formulas and arithmetics are fundamental...

A physicist, a biologist and a mathematician were at a cafe when across the street two people entered a house. Moments later three people came out. The physicist said, "Hmm, that must be a measurement error." The biologist wondered, "It must be reproduction!" And the mathematician said, "If someone goes into the house, it will be empty again."









Evaluating a mathematical formula

Height of a ball in vertical motion.

$$ y(t) = v_0t- \frac{1}{2}gt^2 $$ where

Task: given \( v_0 \), \( g \) and \( t \), compute \( y \).









Use a calculator? A program is much more powerful!

What is a program?

A sequence of instructions to the computer, written in a programming language, somewhat like English, but very much simpler - and very much stricter.

This course teaches the Python language.

Our first example program:

Evaluate \( y(t) = v_0t- \frac{1}{2}gt^2 \) for \( v_0=5 \), \( g=9.81 \) and \( t=0.6 \): $$ y = 5\cdot 0.6 - \frac{1}{2}\cdot 9.81 \cdot 0.6^2$$

The complete Python program:

print 5*0.6 - 0.5*9.81*0.6**2









How to write and run the program

Step 1. Write the program in a text editor, here the line

print 5*0.6 - 0.5*9.81*0.6**2

Step 2. Save the program to a file (say) ball1.py. (.py denotes Python.)

Step 3. Move to a terminal window and go to the folder containing the program file.

Step 4. Run the program:

Terminal> python ball1.py

The program prints out 1.2342 in the terminal window.









In this course we probably use computers differently from what you are used to









A short program can calculate any integral

You cannot calculate this integral by hand: $$ \int_{-\infty}^1 e^{-x^2}dx\tp$$

A little program can compute this and "all" other integrals:

from numpy import *

def integrate(f, a, b, n=100):
    """
    Integrate f from a to b,
    using the Trapezoidal rule with n intervals.
    """
    x = linspace(a, b, n+1)    # Coordinates of the intervals
    h = x[1] - x[0]            # Interval spacing
    I = h*(sum(f(x)) - 0.5*(f(a) + f(b)))
    return I

# Define my special integrand
def my_function(x):
    return exp(-x**2)

minus_infinity = -20  # Approximation of minus infinity
I = integrate(my_function, minus_infinity, 1, n=1000)
print 'Value of integral:', I

The program computes an approximation with error \( 10^{-12} \) within 0.1 s (\( n=10^6 \))!









Computers are very picky about grammar rules and typos

Look at

print 5*0.6 - 0.5*9.81*0.6**2
write 5*0,6 - 0,5*9,81*0,6^2

Would you consider these two lines to be equal?









Programming opens up a new life

People only become computer programmers if they're obsessive about details, crave power over machines, and can bear to be told day after day exactly how stupid they are. G. J. E. Rawlins









Store numbers in variables to make a program more readable

From mathematics you are used to variables, e.g., $$ v_0=5,\quad g=9.81,\quad t=0.6,\quad y = v_0t -\frac{1}{2}gt^2$$

We can use variables in a program too, and this makes the last program easier to read and understand:

v0 = 5
g = 9.81
t = 0.6
y = v0*t - 0.5*g*t**2
print y

This program spans several lines of text and use variables, otherwise the program performs the same calculations and gives the same output as the previous program









There is great flexibility in choosing variable names

initial_velocity = 5
accel_of_gravity = 9.81
TIME = 0.6
VerticalPositionOfBall = initial_velocity*TIME - \ 
                         0.5*accel_of_gravity*TIME**2
print VerticalPositionOfBall

(Note: the backslash allows an instruction to be continued on the next line)

Good variable names make a program easier to understand!









Some words are reserved in Python

Certain words have a special meaning in Python and cannot be used as variable names. These are: and, as, assert, break, class, continue, def, del, elif, else, except, exec, finally, for, from, global, if, import, in, is, lambda, not, or, pass, print, raise, return, try, with, while, and yield.









Comments are useful to explain how you think in programs

Program with comments:

# program for computing the height of a ball
# in vertical motion
v0 = 5    # initial velocity
g = 9.81  # acceleration of gravity
t = 0.6   # time
y = v0*t - 0.5*g*t**2  # vertical position
print y

Note:









Comments are not always ignored....

Normal rule: Python programs, including comments, can only contain characters from the English alphabet.

hilsen = 'Kjære Åsmund!'  # er æ og Å lov i en streng?
print hilsen

leads to an error:

SyntaxError: Non-ASCII character ...

Remedy: put this special comment line as the first line in your program

# -*- coding: utf-8 -*-

Or stick to English everywhere in a program









The printf syntax gives great flexibility in formatting text with numbers

Output from calculations often contain text and numbers, e.g.,

At t=0.6 s, y is 1.23 m.

We want to control the formatting of numbers: no of decimals, style: 0.6 vs 6E-01 or 6.0e-01. So-called printf formatting is useful for this purpose:

t = 0.6;  y = 1.2342
print 'At t=%g s, y is %.2f m.' % (t, y)

The printf format has "slots" where the variables listed at the end are put: %g \( \leftarrow \) t, %.2f \( \leftarrow \) y









Examples on different printf formats

%g         most compact formatting of a real number
%f         decimal notation (-34.674)
%10.3f     decimal notation, 3 decimals, field width 10
%.3f       decimal notation, 3 decimals, minimum width
%e or %E   scientific notation (1.42e-02 or 1.42E-02)
%9.2e      scientific notation, 2 decimals, field width 9
%d         integer
%5d        integer in a field of width 5 characters
%s         string (text)
%-20s      string, field width 20, left-adjusted
%%         the percentage sign % itself

(See the the book for more explanation and overview)









Using printf formatting in our program

Triple-quoted strings (""") can be used for multi-line output, and here we combine such a string with printf formatting:

v0 = 5
g = 9.81
t = 0.6
y = v0*t - 0.5*g*t**2

print """
At t=%f s, a ball with
initial velocity v0=%.3E m/s
is located at the height %.2f m.
""" % (t, v0, y)

Running the program:

Terminal> python ball_print2.py

At t=0.600000 s, a ball with
initial velocity v0=5.000E+00 m/s
is located at the height 1.23 m.









Some frequently used computer science terms

Computer science meaning of terms is often different from the human language meaning









A program consists of statements

a = 1      # 1st statement (assignment statement)
b = 2      # 2nd statement (assignment statement)
c = a + b  # 3rd statement (assignment statement)
print c    # 4th statement (print statement)

Normal rule: one statement per line, but multiple statements per line is possible with a semicolon in between the statements:

a = 1;  b = 2;  c = a + b;  print c









Assignment statements evaluate right-hand side and assign the result to the variable on the left-hand side

myvar = 10
myvar = 3*myvar   # = 30









Syntax is the exact specification of instructions to the computer

Programs must have correct syntax, i.e., correct use of the computer language grammar rules, and no misprints!

This is a program with two syntax errors:

myvar = 5.2
prinnt Myvar

    prinnt Myvar
               ^
SyntaxError: invalid syntax

Only the first encountered error is reported and the program is stopped (correct the error and continue with next error)

Programming demands significantly higher standard of accuracy. Things don't simply have to make sense to another human being, they must make sense to a computer. Donald Knuth, computer scientist, 1938-









Blanks (whitespace) can be used to nicely format the program text

Blanks may or may not be important in Python programs. These statements are equivalent (blanks do not matter):

v0=3
v0  =  3
v0=   3
v0 = 3

Here blanks do matter:

counter = 1
while counter <= 4:
    counter = counter + 1   # correct (4 leading blanks)

while counter <= 4:
counter = counter + 1       # invalid syntax

(more about this in Ch. 2)









A program takes some known input data and computes some output data

v0 = 3;  g = 9.81;  t = 0.6
position = v0*t - 0.5*g*t*t
velocity = v0 - g*t
print 'position:', position, 'velocity:', velocity









An operating system (OS) is a set of programs managing hardware and software resources on a computer









Evaluating a formula for temperature conversion

Given \( C \) as a temperature in Celsius degrees, compute the corresponding Fahrenheit degrees \( F \): $$ F = \frac{9}{5}C + 32 $$

Program:

C = 21
F = (9/5)*C + 32
print F

Execution:

Terminal> python c2f_v1.py
53









We must always check that a new program calculates the right answer

Using a calculator:

9/5 times 21 plus 32 is 69.8, not 53.









The error is caused by (unintended) integer division

Corrected program (with correct output 69.8):

C = 21
F = (9.0/5)*C + 32
print F









Everything in Python is an object

Variables refer to objects:

a = 5       # a refers to an integer (int) object
b = 9       # b refers to an integer (int) object
c = 9.0     # c refers to a real number (float) object
d = b/a     # d refers to an int/int => int object
e = c/a     # e refers to float/int => float object
s = 'b/a=%g' % (b/a)  # s is a string/text (str) object

We can convert between object types:

a = 3              # a is int
b = float(a)       # b is float 3.0
c = 3.9            # c is float
d = int(c)         # d is int 3
d = round(c)       # d is float 4.0
d = int(round(c))  # d is int 4
d = str(c)         # d is str '3.9'
e = '-4.2'         # e is str
f = float(e)       # f is float -4.2









Arithmetic expressions are evaluated as you have learned in mathematics









Standard mathematical functions are found in the math module

Compute \( \sqrt{2} \) using the sqrt function in the math module:

import math
r = math.sqrt(2)
# or
from math import sqrt
r = sqrt(2)
# or
from math import *   # import everything in math
r = sqrt(2)









Another example on computing with functions from math

Evaluate $$ Q = \sin x\cos x + 4\ln x $$ for \( x=1.2 \).

from math import sin, cos, log
x = 1.2
Q = sin(x)*cos(x) + 4*log(x)   # log is ln (base e)
print Q









Computers have inexact arithmetics because of rounding errors

Let us compute \( 1/49\cdot 49 \) and \( 1/51\cdot 51 \):

v1 = 1/49.0*49
v2 = 1/51.0*51
print '%.16f %.16f' % (v1, v2)

Output with 16 decimals becomes

 0.9999999999999999 1.0000000000000000









Test that a calculation is correct

What is printed?

a = 1; b = 2;
computed = a + b
expected = 3
correct = computed == expected
print 'Correct:', correct

Change to a = 0.1 and b = 0.2 (expected = 0.3). What is now printed? Why? How can the comparison be performed?









Answer to exercise on previous slide: use equality test with tolerance!

>>> a = 0.1; b = 0.2; expected = 0.3
>>> a + b == expected
False

>>> print '%.17f\n%.17f\n%.17f\n%.17f' % (0.1, 0.2, 0.1 + 0.2, 0.3)
0.10000000000000001
0.20000000000000001
0.30000000000000004
0.29999999999999999

>>> a = 0.1; b = 0.2; expected = 0.3
>>> computed = a + b
>>> diff = abs(expected - computed)
>>> tol = 1E-15
>>> diff < tol









Another example involving math functions

The \( \sinh x \) function is defined as $$\sinh (x) = \frac{1}{2}\left(e^{x} - e^{-x}\right) $$

We can evaluate this function in three ways:

  1. math.sinh
  2. combination of two math.exp
  3. combination of two powers of math.e

from math import sinh, exp, e, pi
x = 2*pi
r1 = sinh(x)
r2 = 0.5*(exp(x) - exp(-x))
r3 = 0.5*(e**x - e**(-x))
print '%.16f %.16f %.16f' % (r1,r2,r3)

Output: r1 is \( 267.744894041016\underline{4369} \), r2 is \( 267.744894041016\underline{4369} \), r3 is \( 267.744894041016\underline{3232} \) (!)









Python can be used interactively as a calculator and to test statements

Terminal> python
Python 2.7.6 (r25:409, Feb 27 2014, 19:35:40)
...
>>> C = 41
>>> F = (9.0/5)*C + 32
>>> print F
105.8
>>> F
105.8

Previous commands can be recalled and edited









Python has full support for complex numbers

>>> a = -2
>>> b = 0.5
>>> s = complex(a, b)  # make complex from variables
>>> s
(-2+0.5j)
>>> s*w                # complex*complex
(-10.5-3.75j)
>>> s/w                # complex/complex
(-0.25641025641025639+0.28205128205128205j)
>>> s.real
-2.0
>>> s.imag
0.5

See the book for additional info









Python can also do symbolic computing

>>> from sympy import *
>>> t, v0, g = symbols('t v0 g')
>>> y = v0*t - Rational(1,2)*g*t**2
>>> dydt = diff(y, t)                     # 1st derivative
>>> dydt
-g*t + v0
>>> print 'acceleration:', diff(y, t, t)  # 2nd derivative
acceleration: -g
>>> y2 = integrate(dydt, t)
>>> y2
-g*t**2/2 + t*v0









SymPy can do a lot of traditional mathematics

>>> y = v0*t - Rational(1,2)*g*t**2
>>> roots = solve(y, t)    # solve y=0 wrt t
>>> roots
[0, 2*v0/g]

>>> x, y = symbols('x y')
>>> f = -sin(x)*sin(y) + cos(x)*cos(y)
>>> simplify(f)
cos(x + y)
>>> expand(sin(x+y), trig=True)  # requires a trigonometric hint
sin(x)*cos(y) + sin(y)*cos(x)









Summary of Chapter 1 (part 1)









Summary of Chapter 1 (part 2)

Mathematical functions like \( \sin x \) and \( \ln x \) must be imported from the math module:

from math import sin, log
x = 5
r = sin(3*log(10*x))

Use printf syntax for full control of output of text and numbers!

Important terms: object, variable, algorithm, statement, assignment, implementation, verification, debugging









Programming is challenging

Alan Perlis, computer scientist, 1922-1990.









Summarizing example: throwing a ball (problem)

We throw a ball with velocity \( v_0 \), at an angle \( \theta \) with the horizontal, from the point \( (x=0,y=y_0) \). The trajectory of the ball is a parabola (we neglect air resistance): $$ y = x\tan\theta - \frac{1}{2v_0}\frac{gx^2}{\cos^2\theta} + y_0 $$









Summarizing example: throwing a ball (solution)

Program:

g = 9.81    # m/s**2
v0 = 15     # km/h
theta = 60  # degrees
x = 0.5     # m
y0 = 1      # m

print """v0    = %.1f km/h
theta = %d degrees
y0    = %.1f m
x     = %.1f m""" % (v0, theta, y0, x)

# convert v0 to m/s and theta to radians:
v0 = v0/3.6
from math import pi, tan, cos
theta = theta*pi/180

y = x*tan(theta) - 1/(2*v0)*g*x**2/((cos(theta))**2) + y0

print 'y     = %.1f m' % y