Content modified under Creative Commons Attribution license CC-BY 4.0, code under BSD 3-Clause License © 2020 R.C. Cooper
Show code cell source
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('fivethirtyeight')
Homework#
Problems Part 1#
Calculate some properties of a rectangular box that is 12.5”\(\times\)11”\(\times\)14” and weighs 31 lbs
a. What is the volume of the box?
b. What is the average density of the box?
c. What is the result of the following logical operation,
volume>1000
(in inches^3)
Use the variables given below,
str1
andstr2
, and check the followinga.
str1<str2
b.
str1==str2
c.
str1>str2
str1 = 'Python'
str2 = 'python'
The following code has an error, fix the error so that the correct result is returned:
y is 20 and x is less than y
x="1"
y=20
if x<y and y==20:
print('y is 20 and x is less than y')
else:
print('x is not less than y')
There is a commonly-used programming question that asks interviewees to build a fizz-buzz result.
Here, you will build a similar program, but use the numbers from the class, 3255: \(3,~2,~5\rightarrow\) “computational”, “mechanics”, “rocks!”. You should print out a list of numbers, if the number is divisible by 3, replace the 3 with “computational”. If the number is divisible by 2, replace with “mechanics”. If the number is divisible by 5, replace the number with “rocks!”. If the number is divisible by a combination, then add both words e.g. 6 is divisible by 3 and 2, so you would print out “computational mechanics”.
Here are the first 20 outputs your program should print,
index |
printed output |
---|---|
0 |
Computational Mechanics Rocks! |
1 |
1 |
2 |
Mechanics |
3 |
Computational |
4 |
Mechanics |
5 |
Rocks! |
6 |
Computational Mechanics |
7 |
7 |
8 |
Mechanics |
9 |
Computational |
10 |
Mechanics Rocks! |
11 |
11 |
12 |
Computational Mechanics |
13 |
13 |
14 |
Mechanics |
15 |
Computational Rocks! |
16 |
Mechanics |
17 |
17 |
18 |
Computational Mechanics |
19 |
19 |
Problems Part 2#
Create a function called
sincos(x)
that returns two arrays,sinx
andcosx
that return the sine and cosine of the input array,x
.a. Document your function with a help file in
'''help'''
b. Use your function to plot sin(x) and cos(x) for x=\(0..2\pi\)
Use a for-loop to create a variable called
A_99
, where every element is the product of the two indices from 0 to 9 e.g. A_99[3,2]=6 and A_99[4,4]=16.
a. Calculate the mean and standard deviation of A_99
b. time your script that creates A_99
takes its mean and standard deviation using %%timeit
Use the two arrays, X and Y, given below to create A_99 using numpy array math rather than a for-loop.
X, Y = np.meshgrid(np.arange(10), np.arange(10))
b. Calculate the mean and standard deviation of A_99
c. time your script that creates A_99
takes its mean and standard deviation using %%timeit
d. create a filled contour plot of X, Y, A_99 contourf plot documentation
The following linear interpolation function has an error. It is supposed to return y(x) given the the two points \(p_1=[x_1,~y_1]\) and \(p_2=[x_2,~y_2]\). Currently, it just returns and error.
def linInterp(x,p1,p2):
'''linear interplation function
return y(x) given the two endpoints
p1=np.array([x1,y1])
and
p2=np.array([x2,y2])'''
slope = (p2[2]-p1[2])/(p2[1]-p1[1])
return p1[2]+slope*(x - p1[1])
Problems Part 3#
The growth of populations of organisms has many engineering and scientific applications. One of the simplest models assumes that the rate of change of the population p is proportional to the existing population at any time t:
\(\frac{dp}{dt} = k_g p\)
where \(t\) is time in years, and \(k_g\) is growth rate in [1/years].
The world population has been increasing dramatically, let’s make a prediction based upon the following data saved in world_population_1900-2020.csv:
year |
world population |
---|---|
1900 |
1,578,000,000 |
1950 |
2,526,000,000 |
2000 |
6,127,000,000 |
2020 |
7,795,482,000 |
a. Use a growth rate of \(k_g=0.013\) [1/years] and compare the analytical solution (use initial condition p(1900) = 1578000000) to the Euler integration for time steps of 20 years from 1900 to 2020 (Hint: use method (1)- plot the two solutions together with the given data)
b. Discussion question: If you decrease the time steps further and the solution converges, will it converge to the actual world population? Why or why not?
Note: We have used a new function np.loadtxt
here. Use the help
or ?
to learn about what this function does and how the arguments can change the output. In the next module, we will go into more details on how to load data, plot data, and present trends.
import numpy as np
year, pop = np.loadtxt('../data/world_population_1900-2020.csv',skiprows=1,delimiter=',',unpack=True)
print('years=',year)
print('population =', pop)
years= [1900. 1950. 2000. 2020.]
population = [1.578000e+09 2.526000e+09 6.127000e+09 7.795482e+09]
print('average population changes 1900-1950, 1950-2000, 2000-2020')
print((pop[1:] - pop[0:-1])/(year[1:] - year[0:-1]))
print('average growth of 1900 - 2020')
print(np.mean((pop[1:] - pop[0:-1])/(year[1:] - year[0:-1])))
average population changes 1900-1950, 1950-2000, 2000-2020
[18960000. 72020000. 83424100.]
average growth of 1900 - 2020
58134700.0
d. What happens to the numerical model as the number of timesteps increases? Does it look more like the analytical curve, \(e^k(t-1900)\), or the measured population data, pop
?
In the freefall example you used smaller time steps to decrease the truncation error in our Euler approximation. Another way to decrease approximation error is to continue expanding the Taylor series. Consider the function f(x)
\(f(x)=e^x = 1+x+\frac{x^2}{2!}+\frac{x^3}{3!}+\frac{x^4}{4!}+...\)
We can approximate \(e^x\) as \(1+x\) (first order), \(1+x+x^2/2\) (second order), and so on each higher order results in smaller error.
a. Use the given
exptaylor
function to approximate the value of exp(1) with a second-order Taylor series expansion. What is the relative error compared tonp.exp(1)
?b. Time the solution for a second-order Taylor series and a tenth-order Taylor series using
%%timeit
. How long would a 100,000-order series take (approximate this, DON’T RUNexptaylor(1, 1000000)
)c. Plot the relative error as a function of the Taylor series expansion order from first order upwards. (Hint: use method (4) in the comparison methods from the “Truncation and roundoff error accumulation in log-log plot” figure)
from math import factorial
def exptaylor(x,n):
'''Taylor series expansion about x=0 for the function e^x
the full expansion follows the function
e^x = 1+ x + x**2/2! + x**3/3! + x**4/4! + x**5/5! +...'''
if n<1:
print('lowest order expansion is 0 where e^x = 1')
return 1
else:
ex = 1+x # define the first-order taylor series result
for i in range(1,n):
ex+=x**(i+1)/factorial(i+1) # add the nth-order result for each step in loop
return ex