Content modified under Creative Commons Attribution license CC-BY 4.0, code under BSD 3-Clause License © 2020 R.C. Cooper

Hide code cell source
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('fivethirtyeight')

Homework#

Problems Part 1#

  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)

  1. Use the variables given below, str1 and str2, and check the following

    a. str1<str2

    b. str1==str2

    c. str1>str2

    d. How could you force (b) to be true? Hint or Hint

str1 = 'Python'
str2 = 'python'
  1. 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')
  1. 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#

  1. Create a function called sincos(x) that returns two arrays, sinx and cosx 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\)

  1. 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

  1. 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

  1. 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#

  1. 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?

  1. 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 to np.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 RUN exptaylor(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