Lab 1: Solving problems using computational methods¶
Welcome to the first lab session for Computational Modelling! In this first lab we will get Python running using Jupyter lab, and start using two of the most important libraries for scientific computing, NumPy and Matplotlib.
Copyright By PowCoder代写 加微信 powcoder
In the first problem we see one of the simplest ways to combine modelling and computation, as we use a given functional dependence to calculate something of interest. To do this we will need to use NumPy and Matplotlib. In the second problem we will first need to develop a simple model, which we will then solve computationally.
Problem 1: To the top of Mount Everest¶
A useful model for the Earth’s atmosphere is the constant temperature hydrostatic model, in which the pressure $P$ decays with height $z$ like
$$P = P_0 \exp(-z/H),$$
$$H = \frac{k_B T}{mg}$$
is the density scale height. This model follows from the model for an ideal gas. To determine the atmospheric pressure at different altitudes we just need to measure the pressure at sea level, $P_0$, and set a temperature $T$. The other parameters are well-known constants: $k_B$ is the Boltzmann constant, $g$ is the acceleration due to gravity and $m$ is the average molecular mass of air.
We would like to see what this model predicts for atmospheric pressure as a function of altitude. In particular, there are two things we would like to determine: the air pressure at the top of Mount Everest, and the maximum altitude humans can survive at, if the minimum pressure at which we can live is $0.47 P_0$. We shall then compare the predictions of our simple model with the experimental data.
To begin with we need to import the NumPy and Matplotlib.pyplot libraries. Use the import command and use the aliases np and plt for the two libraries.
Given that $k_B = 1.38 \times 10^{-23}$ kgm$^2$s$^{-2}$K$^{-1}$, $g \approx 9.8$ ms$^{-2}$, and $m\approx 29 \times 1.67 \times 10^{-27}$ kg, determine $H$:
#Code for H here
Using the value $P_0 = 101.3$ kPa, determine the pressure at the summit of Mount Everest (which has a height of 8848 m) in kPa. To calculate this we need to use the exponential function. NumPy has all the common mathematical functions defined, so we just need to call it using the standard Python format library.method, so in this case np.exp()
#Code for P here
The experimentally measured value of the atmospheric pressure at the top of Mount Everest is $\approx 30$ kPa. Comment on the success of the model in predicting this, and identify any reasons you can see for a discrepancy.
Add your comments here:
Our goal now is to visualise the pressure as a function of altitude. To do this we first need to obtain the appropriate data from the model. NumPy has a number of functions for creating vectors, so we can use these to first set up a height vector, and from this calculate the associated pressures. The functions linspace() and arange() are both useful for our purposes. Use help() to find out how these functions work, and then set up a vector of heights from 0 to 30,000 metres, in steps of 1 metre, and then a vector of associated pressures.
# Use e.g. help(np.linspace)
# Put your vectors here
We now need to visualise the dependence of $P$ on $z$. To do this we will use routines from the Matplotlib library. A typical process for a quick display is to use the routine plot followed by show. In this case we will also add axes labels using xlabel and ylabel before showing the plot.
# Plot Pressure vs height using Matplotlib routines
# Remember to reference all commands using the Matplotlib alias
Replace the plot command with semilogy, which plots the $y$ axis with logarithmic scaling. What do you observe?
#Plot with logarithmic y axis here
Before we leave the model of atmospheric pressure, what is the prediction of maximum altitude at which humans can survive, if humans can tolerate a pressure of $0.47 P_0$? There are few ways we could solve this, but here we will try to answer using an interactive plot. Unfortunately interactive plots only work with jupyter notebook not jupyter lab. If you’re using Jupyter Lab you will need to save and exit, and reopen with Jupyter Notebook. The command %matplotlib notebook starts Matplotlib in interactive mode. We then need to obtain our pressure versus altitude diagram again, and at the same time, plot the line corresponding to the minimum pressure. The NumPy command ones is useful for creating a vector of constant values. We want the vector to be the same length as the other $P$ vector, and for this the Python function len() is useful.
# %matplotlib notebook Command to start matplotlib plots in interactive mode, along with plt.ion()
# Use len() and np.ones() to set up constant data
%matplotlib notebook
The maximum habitable altitude has been observed to be approximately 6000m. Is the model prediction reasonable?
Checkpoint 1
You have reached the first checkpoint. You need to replot your previous figure, label the axes and include a legend. Possible labels for the lines could be $P$ and $0.47 P_0$. There are a few useful Matplotlib commands that you could use here to produce a publishable graph. First you could clear the previous axes with `clf()`. You could include a legend with `legend([‘label1′,’label2’])` and you could save your figure with `savefig(‘name_of_file.png’)`. Note, to go back to inline (i.e. non-interactive plots) we have to use the command `%matplotlib inline`. Submit your saved png file to the Checkpoint1 submission link in Lab 1 on Canvas.
# Output your final figure for Checkpoint 1 here
%matplotlib inline
We now turn to the second problem in the lab. In this problem we will first need to develop an appropriate representation of our problem. Once we have a model, we can then use the computer to rapidly explore the possible outcomes, and in this way find an optimal solution. This is a common computational modelling task.
Problem 2: Lifeguard to the Rescue!¶
A lifeguard on duty at the beach spots a swimmer in distress. The lifeguard is 20m from the water, and the swimmer is 30m out to sea from the beach, however the lifeguard is 50m down the beach from the swimmer. If the lifeguard could swim and run at the same speed, then the fastest path to the swimmer is the straight line between the lifeguard and the swimmer. However, the lifeguard swims at 1.5 m/s, but runs over the sand at 5 m/s. The question is, how far down the beach should the lifeguard enter the water, so as to minimise the time taken to reach the swimmer?
# Import your libraries here (not necessary, if already imported for first problem)
A good first step is to sketch out a diagram of the situation on a piece of paper. A second step is to introduce suitable notation, and then identify the knowns and the unknowns.
Let’s begin by calculating the time it will take the lifeguard to reach the swimmer when following the path of least distance. Using the relationship $\text{distance} = \text{speed} \times \text{time}$, determine this time below.
You may need to use mathematical functions to do this. NumPy has all the common functions defined. For instance, to calculate the arctangent of a value we call arctan(). Similarly the cosine of an angle is calculated by a call to cos(). You may need to use the Python command print() to print out the value of a variable.
#Write code to calculate time for path of least distance here
The key decision the lifeguard needs to make is where to enter the water. If we define the distance down the beach from the lifeguard as $x$, such that $x=0$ corresponds to the lifeguard entering the water at the closest point to the start point, and $x = 50$ the lifeguard entering the water at the closest point to the swimmer, we can determine how the rescue time depends on $x$.
Solve this problem by setting up a vector of $x$ values, from $x=0$ to $x=50$ and then calculating the rescue time for each $x$ value. Plot these rescue times versus $x$ on a graph, and label the axes.
#Write code to plot rescue time versus x here
We can see in the plot that there is a minimum rescue time. Let’s find this minimum using NumPy routines. The command argmin() returns the index (i.e. the vector address) corresponding to the minimum value in the vector. Use this command to find the value of $x$ corresponding to the minimum rescue time.
# Write code to calculate the value of x which minimises the rescue time, and print out the value
In the 19th century, the French scientist Fermat suggested an elegant method for solving problems
in ray optics in the form of his Principle of Least Time, which states “Light travelling between two points in an optical system follows the path which minimizes the travel time”.
From this follows Snell’s law, which relates the angle at which a beam of light bends when it passes from one medium to another. Snell’s Law states that
$$ \frac{\sin\theta_1}{\sin\theta_2} = \frac{v_1}{v_2}.$$
Calculate the two sides of this equation for the case of the lifeguard travelling along the minimum time path. Does the lifeguard follow Snell’s law?
# Implement Snell’s law here
COSC2002 Checkpoint 2
Save your figure of the rescue time versus $x$ (with labelled axes), and include on the plot $\text{ratio} = \sin\theta_1/\sin\theta_2$, where the right-hand side of the equation is the value you calculated.
Save the resulting figure as a png file and submit to Checkpoint 2 on Canvas.
#Code for saving figure here
#Use plt.text(20,40,’Ratio = ‘+str(ratio)) to place text on the plot
COSC2902 Checkpoint 2
Fermat’s principle of least time also applies for the optics of reflection. In our lifeguard
analogy, this corresponds to the lifeguard running from their start point, touching the waterline, and
then running to the point 20 m from the water, but 50 m further down the beach. Note that the lifeguard doesn’t enter the ocean: the water simply serves as a convenient boundary.
Write a program to determine the path of least time as a function of $x$. Plot the travel time versus $x$, labelling the axes, and include text indicating the value of $x$ which minimises the travel time. Save the figure as a png file and submit it to Checkpoint 2 on Canvas. What does the value of $x$ minimising the travel time imply for the ratio of the angle incidence to the angle of reflection?
# Code for saving figure here
# Use plt.text(20,14,’xmin = ‘+str(xmin)) to place text on the plot
程序代写 CS代考 加微信: powcoder QQ: 1823890830 Email: powcoder@163.com