{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Introduction to programming for Geoscientists (through Python)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Solutions to Lecture 1 Exercises: Computing with formulas\n", "## Gerard J. Gorman (g.gorman@imperial.ac.uk) http://www.imperial.ac.uk/people/g.gorman" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* **Compute the growth of money in a bank**
\n", "Let *p* be a bank's interest rate in percent per year. An initial amount *A* has then grown to $$A\\left(1+\\frac{p}{100}\\right)^n$$ after *n* years. Write a program for computing how much money 1000 euros have grown to after three years with a 5% interest rate." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The amount of money in the account after 3 years is: 1157.63 euros\n" ] } ], "source": [ "# Solution\n", "\n", "p = 5 # Interest rate in percent per year\n", "A = 1000.0 # Initial amount of money in euros\n", "n = 3 # Number of years since initial deposit\n", "\n", "amount = A*(1 + p/100.0)**n\n", "\n", "print(\"The amount of money in the account after %d years is: %.2f euros\" % (n, amount))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* **Find error(s) in a program**
\n", "Suppose somebody has written the simple one-line program below for computing $\\sin(1)$. Try to run this program. What is the problem? Try to fix it. (Reminder - the math function needs to be imported)." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "ename": "SyntaxError", "evalue": "invalid syntax (, line 2)", "output_type": "error", "traceback": [ "\u001b[0;36m File \u001b[0;32m\"\"\u001b[0;36m, line \u001b[0;32m2\u001b[0m\n\u001b[0;31m x=1; print 'sin(%g)=%g' % (x, sin(x))\u001b[0m\n\u001b[0m ^\u001b[0m\n\u001b[0;31mSyntaxError\u001b[0m\u001b[0;31m:\u001b[0m invalid syntax\n" ] } ], "source": [ "# Broken program\n", "x=1; print 'sin(%g)=%g' % (x, sin(x))" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "sin(1)=0.841471\n" ] } ], "source": [ "# Fixed program (Solution). You must import the 'sin' function from the 'math' module.\n", "from math import sin\n", "x=1; print('sin(%g)=%g' % (x, sin(x)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* **Evaluate a Gaussian function**
\n", "The bell-shaped Gaussian function,\n", "$$f(x)=\\frac{1}{\\sqrt{2\\pi}s}\\exp\\left(-\\frac{1}{2} \\left(\\frac{x-m}{s}\\right)^2\\right)$$\n", "is one of the most widely used functions in science and technology. The parameters $m$ and $s$ are real numbers, where $s$ must be greater than zero. Write a program for evaluating this function when $m = 0$, $s = 2$, and $x = 1$. Verify the program's result by comparing with hand calculations on a calculator." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The value of the Gaussian function with x = 1.00 is 0.1760\n" ] } ], "source": [ "# Solution\n", "\n", "from math import exp, pi, sqrt\n", "\n", "# Parameters.\n", "m = 0.0 # The mean\n", "s = 2.0 # The standard deviation\n", "x = 1.0\n", "coefficient = 1.0/(sqrt(2*pi)*s)\n", "\n", "result = coefficient*exp(-0.5*((x - m)/s)**2)\n", "\n", "print(\"The value of the Gaussian function with x = %.2f is %.4f\" % (x, result))" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The value of the Gaussian function with x = 1.00 is 0.1760\n" ] } ], "source": [ "# Alternative Solution. This uses two concepts that will be covered\n", "# in later lectures: 1) Python functions, and 2) casting.\n", "\n", "from math import exp, pi, sqrt\n", "\n", "def f(x):\n", " # Parameters. These could be given to the function as arguments.\n", " m = 0 # The mean\n", " s = 2 # The standard deviation\n", " coefficient = 1/(sqrt(2*pi)*s)\n", " result = coefficient*exp(-0.5*(float(x - m)/s)**2)\n", " return result\n", "\n", "x = 1\n", "print(\"The value of the Gaussian function with x = %.2f is %.4f\" % (x, f(x)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* **Compute the air resistance on a football**
\n", "The drag force, due to air resistance, on an object can be expressed as\n", "$$F_d = \\frac{1}{2}C_D\\rho AV^2$$\n", "where $\\rho$ is the density of the air, $V$ is the velocity of the object, $A$ is the cross-sectional area (normal to the velocity direction), and $C_D$ is the drag coefficient, which depends heavily on the shape of the object and the roughness of the surface.

\n", "The gravity force on an object with mass $m$ is $F_g = mg$, where $g = 9.81ms^{−2}$.

\n", "Write a program that computes the drag force and the gravity force on an object. Write out the forces with one decimal in units of Newton ($N = kgm/s^2$). Also print the ratio of the drag force and the gravity force. Define $C_D$, $\\rho$, $A$, $V$, $m$, $g$, $F_d$, and $F_g$ as variables, and put a comment with the corresponding unit.

\n", "As a computational example, you can initialize all variables with values relevant for a football kick. The density of air is $\\rho = 1.2 kg m^{−3}$. For any ball, we have obviously that $A = \\pi a^2$, where $a$ is the radius of the ball, which can be taken as $11cm$ for a football. The mass of the ball is $0.43kg$. $C_D$ can be taken as $0.2$.

\n", "Use the program to calculate the forces on the ball for a hard kick, $V = 120km/h$ and for a soft kick, $V = 10km/h$ (it is easy to make the mistake of mixing inconsistent units, so make sure you compute with V expressed in m/s)." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The force due to gravity is 4.2 N\n", "The drag force acting on the ball after a hard kick is 5.1 N\n", "The ratio of the drag force and the gravitational force is 1.2\n", "The drag force acting on the ball after a soft kick is 0.0 N\n", "The ratio of the drag force and the gravitational force is 0.0\n" ] } ], "source": [ "# Solution\n", "\n", "from math import pi\n", "\n", "# The density of air, in kg/(m**3)\n", "density = 1.2\n", "\n", "# Mass of the ball, in kg\n", "mass = 0.43\n", "\n", "# Radius of the ball, in m\n", "a = 0.11\n", "\n", "# Cross-sectional area (m**2)\n", "A = pi*(a**2)\n", "\n", "# Non-dimensional\n", "drag_coefficient = 0.2\n", "\n", "# Acceleration due to gravity, in m/(s**2)\n", "g = 9.81 \n", "\n", "# Force due to gravity - independent of the velocity\n", "F_g = mass*g\n", "print(\"The force due to gravity is %.1f N\" % F_g)\n", "\n", "# Velocity of the ball after a hard kick, converted from km/h to m/s\n", "V_hard = 33.333333333 \n", "\n", "# Drag force acting on the ball\n", "F_d = 0.5*drag_coefficient*density*A*(V_hard**2) \n", "print(\"The drag force acting on the ball after a hard kick is %.1f N\" % F_d)\n", "print(\"The ratio of the drag force and the gravitational force is %.1f\" % (F_d/F_g))\n", "\n", "# Velocity of the ball after a hard kick, converted from km/h to m/s\n", "V_soft = 2.777777778 \n", "\n", "# Drag force acting on the ball\n", "F_d = 0.5*drag_coefficient*density*A*(V_soft**2)\n", "\n", "# The drag force will be so small that this prints out as \"0.0 N\"\n", "print(\"The drag force acting on the ball after a soft kick is %.1f N\" % F_d)\n", "print(\"The ratio of the drag force and the gravitational force is %.1f\" % (F_d/F_g))\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* **Find errors in the coding of a formula**
\n", "Given a quadratic equation,\n", "$$ax^2 + bx + c = 0,$$\n", "$$x1 = −b+\\frac{\\sqrt{b^2 −4ac}}{2a},$$ and\n", "$$x2 = −b−\\frac{\\sqrt{b^2 −4ac}}{2a}.$$\n", "Why does the following program not work correctly?" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "ename": "ValueError", "evalue": "math domain error", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mValueError\u001b[0m Traceback (most recent call last)", "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m()\u001b[0m\n\u001b[1;32m 2\u001b[0m \u001b[0ma\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m2\u001b[0m\u001b[0;34m;\u001b[0m \u001b[0mb\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m;\u001b[0m \u001b[0mc\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m2\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 3\u001b[0m \u001b[0;32mfrom\u001b[0m \u001b[0mmath\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0msqrt\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 4\u001b[0;31m \u001b[0mq\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msqrt\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mb\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mb\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0;36m4\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mc\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 5\u001b[0m \u001b[0mx1\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0mb\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mq\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 6\u001b[0m \u001b[0mx2\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0mb\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mq\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;31mValueError\u001b[0m: math domain error" ] } ], "source": [ "# Original program\n", "a = 2; b = 1; c = 2\n", "from math import sqrt\n", "q = sqrt(b*b - 4*a*c)\n", "x1 = (-b + q)/2*a\n", "x2 = (-b - q)/2*a\n", "print(x1, x2)" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "-0.5 -2.0\n" ] } ], "source": [ "# Solution\n", "a = 2; b = 5; c = 2\n", "from math import sqrt\n", "\n", "# The error above was caused by taking the square root of a negative number. Try changing b to 5.\n", "q = sqrt(b*b - 4*a*c)\n", "\n", "# You need to make sure that -b is divided by 2*a, and also include brackets in the denominator...\n", "x1 = (-b + q)/(2*a)\n", "\n", "# ...and again here. This is because of operator precedence.\n", "x2 = (-b - q)/(2*a)\n", "print(x1, x2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* **Convert from meters to British length units**
\n", "Make a program where you set a length given in meters and then compute and write out the corresponding length measured in inches, in feet, in yards, and in miles. Use the fact that one inch is 2.54 cm, one foot is 12 inches, one yard is 3 feet, and one British mile is 1760 yards. As a verification, a length of 640 meters corresponds to 25196.85 inches, 2099.74 feet, 699.91 yards, or 0.3977 miles." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "640.0000 meters = 25196.8504 inches\n", "640.0000 meters = 2099.7375 feet\n", "640.0000 meters = 699.9125 yards\n", "640.0000 meters = 0.3977 miles\n" ] } ], "source": [ "# Solution\n", "\n", "meters = 640\n", "\n", "# 1 inch = 2.54 cm. Remember to convert from 2.54 cm to 0.0254 m here.\n", "inches = meters/0.0254 \n", "print(\"%.4f meters = %.4f inches\" % (meters, inches))\n", "\n", "feet = inches/12.0 # 1 foot = 12 inches\n", "print(\"%.4f meters = %.4f feet\" % (meters, feet))\n", "\n", "yards = feet/3.0 # 1 yard = 3 feet\n", "print(\"%.4f meters = %.4f yards\" % (meters, yards))\n", "\n", "miles = yards/1760.0 # 1 yard = 1760 miles\n", "print(\"%.4f meters = %.4f miles\" % (meters, miles))" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.2" } }, "nbformat": 4, "nbformat_minor": 1 }