Sci Programing 2

Follow-up on Scientific programming as of Wednesday (or last week) Victor Miclovich November 11, 2009

1 y(t) = v0 t − gt2 2 The above expression that says a lot about the trajectory of a ball. It actually computes the height of a ball in vertical motion. What if we throw the ball with an initial velocity having an angle θ with the horizontal ”plane”? This kind of problem has been encountered by many of you... at least if you went to a high school physics class. The ball will follow a trajectory y = f (x) through the air, where f (x) = x tan θ −

1 gx2 + y0 2vo2 cos2 θ

In this expression, is a horizontal coordinate... the distance in the direction of the x axis, g is the acceleration of gravity, v0 is the size1 of the initial velocity which makes an angle θ with the x axis, and (0, y0 ) is the initial position of the ball. Our programming goal is to make a program evaluating equation above. The program should write out the value of all the involved variables and what their units are. Solution: We use the SI System and assume that v0 has been given in km/h; g = 9.81m/s2 ; x,y, and y0 are measure in meters; and θ in degrees. The program has got four parts; that are natural. • Initialization of input data • Importation of functions and π from math • conversion of v0 and θ to m/s and radians, respectively. 1

I use size or magnitude here, velocity has speed and direction, don’t forget simple Physics!

1

• And, finally, the evaluation of the right hand-side of the expression. The complete program is shown below2 : program 1 Python solving for y g = 9.81 #m/s**2 v0 = 15 # km/h theta = 60 # degrees x = 0.5 # m y0 = 1 # m print v0 theta y0 """ %

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

from math import pi, tan, cos # convert v0 to m/s # and theta to radians: v- = v0/3.6 theta = theta*pi/180 y = (x*tan(theta) - 1/(2*v0**2)*g*x**2 / ((cos(theta))**2) + y0) print ’y = %.1f m’ % y

Explanation of code The backslash in the triple-quoted multi-line string makes the string continue on the next line without a newline. This means that removing the backslash results in a blank line above the v0 line and a blank line between x and y lines in t he output on the screen. Another point to note is the expression 1/(2*v0**2), which might seem as a candidate for unintended integer division. However, the conversion of v0 to m/s involves a division by 3.6, which results in v0 being float, and therefore 2*v0**2 being float. We shall execute this scripts using Spyder3 . 2

We choose to write out all numerical values with one decimal I hope to configure all these machines with PyLab (plotting), SciPy(signal and image processing), Matplotlib (2D and 3D graphing) and Numpy (linear Algebra) packages... so you won’t have any excuses about heavy course works when you can do things with Python. 3

2

Some pointers to take note of: • Understand the problem you want to solve • Understand what the program is supposed to do • Understand how to translate the problem description into a series of Python statements • Another equally important step is: Testing! A complete solution to programming exercises therefore consists of two parts: the program text and a demonstration that the program works correctly.

Some Python words to know • object: anything that a variable (name) can refer to4 (number, string, function, module,...) • variable: name of an object • statement: an instruction to the computer, usually written on a line in a Python program (multiple statements on a line must be separated by semicolons) • expression: a combination of numbers, text, variables, and operators that results in a new object, when being evaluated • assignment: a statement binding an evaluated expression (object) to a variable (name) • algorithm: detailed recipe for how to solve a problem by programming • code: program text (or synonym for program) • implementation: same as code • executable: the file we run to start a program • verification: providing evidence that the program works correctly • debugging: locating and correcting errors in a program.

4

But objects can exist without being bound to a name:print ’Hello!’ first makes a string object of the text in quotes and then the contents of this string object, without a name, is printed

3

How to cook the perfect egg As an egg cooks, the proteins first denature and then coagulate. When the temperature exceeds a critical point, reactions begin and proceed faster as the temperature increases. In the egg white the proteins start to coagulate for temperatures above 70◦ C. For a hard boiled egg, the center of the yolk should be allowed to reach 70◦ C. The following formula expresses the time t it takes (in seconds) for the center of the yolk to reach the temperature Ty (in Celsius degrees): 2

1



T0 − Tw t= 2 ln 0.76 Ty − Tw 3 Kπ 2 ( 4π 3 ) M 3 cρ 3



Here, M , ρ, c, and K are properties of the egg: M is the mass, ρ is the density, c is the specific heat capacity, and K is thermal conductivity. Relevant values are M = 47 g for a small egg and M = 67 g for a large egg, ρ = 1.038g cm−3 , c = 3.7Jg −1 K −1 , and K = 5.4 × 10−3 W cm−1 K −1 . Furthermore, Tw is the temperature (in C degrees) of the boiling water, and T0 is the original temperature (in C degrees) of the egg before being put in the water. Implement the formula in a program, set Tw = 100 C and Ty = 70 C, and compute t for the large egg taken from the fridge (T0 = 4 C) and from room temperature (T0 = 20 C)5 . Name of program should be: egg.py

Other curious stuff ! The piecewise constant function  0 < t < T2 ,  1, f (x) = 0, t = T2 ,  −1, T2 < t < T can be approximated by the sum n 2(2i − 1)πt 4X 1 S(t; n) = sin π 2i − 1 T i=1

4 = π



2πt 1 6πt 1 10πt sin + sin + sin + ... T 3 T 5 T



It can be shown that S(t; n) → f (t) as n → ∞. Write a program that prints out the value of S(αT ; n) for α = 0.01, T = 2π, and n = 1, 2, 3, 4. 5

Notice that I mention two different ”original” temperatures, your code should answer t for different values of T0

4

Let S(=S(t; n)), t, alpha, and T be variables in the program. A new S, corresponding to a new n, should be computed by adding one term to the previous value of S, i.e., by a statement like S=S + term. Run the program also for α = 14 . Does the approximation S(αT ; 4) seem to be better for α = 1/4 than for α = 0.01? Name the program file: compare fun sum.py6

What next? Over the next weeks, months or years, we shall be looking at various engineering functions, interesting mathematics and the like... just be prepared to use this tool in your work too!

6

Email me your answers: [email protected]

5