{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Homework 7\n",
"\n",
"Student's Name: PLEASE INSERT YOUR NAME HERE"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"__Directions:__ Add work to this notebook to solve the problems below. Some other cells have been included which are important for testing your work and should give you some feeback.\n",
"\n",
"Check your work. I give some tests that should be passed if you get it correct.\n",
"\n",
"**Your notebook should run without printing errors and without user input. If this is not the case, points will be deducted from your grade.**\n",
"\n",
"Problem Sources:\n",
"* __LL__: *Programming for Computations - Python* by Svein Linge and Hans Petter Langtangen, 2nd edition.\n",
"* __L__: *A Primer on Scientific Programming with Python* by Hans Petter Langtangen, 2nd edition.\n",
"* __TAK__: *Applied Scientific Computing With Python* by Peter R. Turner, Thomas Arildsen, and Kathleen Kavanagh."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Standard imports\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"import math as m\n",
"from mpmath import mp, iv"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Data Source:**\n",
"\n",
"Monthly $CO_2$ measurements made at the Mauna Lao Observatory are available from the [co2.earth website](https://www.co2.earth/monthly-co2). I have included the important information from the table of that data for years begining with 2010 and placed it in the following dictionary. (I extracted this data from [this text file](ftp://aftp.cmdl.noaa.gov/products/trends/co2/co2_mm_mlo.txt).)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"data = { 2010.042: 388.7,\n",
" 2010.125: 390.2,\n",
" 2010.208: 391.17,\n",
" 2010.292: 392.46,\n",
" 2010.375: 393.0,\n",
" 2010.458: 392.15,\n",
" 2010.542: 390.2,\n",
" 2010.625: 388.35,\n",
" 2010.708: 386.85,\n",
" 2010.792: 387.24,\n",
" 2010.875: 388.67,\n",
" 2010.958: 389.79,\n",
" 2011.042: 391.33,\n",
" 2011.125: 391.86,\n",
" 2011.208: 392.6,\n",
" 2011.292: 393.25,\n",
" 2011.375: 394.19,\n",
" 2011.458: 393.74,\n",
" 2011.542: 392.51,\n",
" 2011.625: 390.13,\n",
" 2011.708: 389.08,\n",
" 2011.792: 389.0,\n",
" 2011.875: 390.28,\n",
" 2011.958: 391.86,\n",
" 2012.042: 393.12,\n",
" 2012.125: 393.86,\n",
" 2012.208: 394.4,\n",
" 2012.292: 396.18,\n",
" 2012.375: 396.74,\n",
" 2012.458: 395.71,\n",
" 2012.542: 394.36,\n",
" 2012.625: 392.39,\n",
" 2012.708: 391.11,\n",
" 2012.792: 391.05,\n",
" 2012.875: 392.98,\n",
" 2012.958: 394.34,\n",
" 2013.042: 395.55,\n",
" 2013.125: 396.8,\n",
" 2013.208: 397.43,\n",
" 2013.292: 398.41,\n",
" 2013.375: 399.78,\n",
" 2013.458: 398.61,\n",
" 2013.542: 397.32,\n",
" 2013.625: 395.2,\n",
" 2013.708: 393.45,\n",
" 2013.792: 393.7,\n",
" 2013.875: 395.16,\n",
" 2013.958: 396.84,\n",
" 2014.042: 397.85,\n",
" 2014.125: 398.01,\n",
" 2014.208: 399.77,\n",
" 2014.292: 401.38,\n",
" 2014.375: 401.78,\n",
" 2014.458: 401.25,\n",
" 2014.542: 399.1,\n",
" 2014.625: 397.03,\n",
" 2014.708: 395.38,\n",
" 2014.792: 396.03,\n",
" 2014.875: 397.28,\n",
" 2014.958: 398.91,\n",
" 2015.042: 399.98,\n",
" 2015.125: 400.28,\n",
" 2015.208: 401.54,\n",
" 2015.292: 403.28,\n",
" 2015.375: 403.96,\n",
" 2015.458: 402.8,\n",
" 2015.542: 401.31,\n",
" 2015.625: 398.93,\n",
" 2015.708: 397.63,\n",
" 2015.792: 398.29,\n",
" 2015.875: 400.16,\n",
" 2015.958: 401.85,\n",
" 2016.042: 402.56,\n",
" 2016.125: 404.12,\n",
" 2016.208: 404.87,\n",
" 2016.292: 407.45,\n",
" 2016.375: 407.72,\n",
" 2016.458: 406.83,\n",
" 2016.542: 404.41,\n",
" 2016.625: 402.27,\n",
" 2016.708: 401.05,\n",
" 2016.792: 401.59,\n",
" 2016.875: 403.55,\n",
" 2016.958: 404.45,\n",
" 2017.042: 406.17,\n",
" 2017.125: 406.46,\n",
" 2017.208: 407.22,\n",
" 2017.292: 409.04,\n",
" 2017.375: 409.69,\n",
" 2017.458: 408.88,\n",
" 2017.542: 407.12,\n",
" 2017.625: 405.13,\n",
" 2017.708: 403.37,\n",
" 2017.792: 403.63,\n",
" 2017.875: 405.12,\n",
" 2017.958: 406.81,\n",
" 2018.042: 407.96,\n",
" 2018.125: 408.32,\n",
" 2018.208: 409.41,\n",
" 2018.292: 410.24,\n",
" 2018.375: 411.24,\n",
" 2018.458: 410.79,\n",
" 2018.542: 408.71,\n",
" 2018.625: 406.99,\n",
" 2018.708: 405.51,\n",
" 2018.792: 406.0,\n",
" 2018.875: 408.02,\n",
" 2018.958: 409.07,\n",
" 2019.042: 410.83,\n",
" 2019.125: 411.75,\n",
" 2019.208: 411.97,\n",
" 2019.292: 413.32,\n",
" 2019.375: 414.66,\n",
" 2019.458: 413.92,\n",
" 2019.542: 411.77,\n",
" 2019.625: 409.94,\n",
" 2019.708: 408.54,\n",
" 2019.792: 408.53,\n",
" 2019.875: 410.27,\n",
" 2019.958: 411.76,\n",
" 2020.042: 413.4,\n",
" 2020.125: 414.11,\n",
" 2020.208: 414.5}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"A dictionary represents a map. Here our map sends a month to the average ppm of $CO_2$ measured that month. Each month is represented by the point in time, the middle of the month measured in years. For example, since `2020.042` is the first entry larger than `2020`, this represents January 2020 and the average of $CO_2$ measurements that month was `413.4 ppm`."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We haven't discussed dictionaries much. But, you can learn about them in L § 6.2. Here we convert the data to numpy arrays:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"months = []\n",
"measurements = []\n",
"for month in data:\n",
" months.append(month)\n",
" measurements.append(data[month])\n",
"months = np.array(months)\n",
"measurements = np.array(measurements)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Below we plot the data."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"plt.plot(months,measurements,\".\")\n",
"plt.xlabel(\"time in years\")\n",
"plt.ylabel(\"average $C0_2$ measurement in ppm\")\n",
"plt.title(\"Monthly C02 measurements\")\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The goal of this homework assignment is to try to fit a curve to this data. Looking at the graph, we see two dominant features:\n",
"* The measurements seem to oscillate with a period of one year. (E.g., there is one local maximum per year.)\n",
"* Each year, there seems to be a linear increase.\n",
"\n",
"This suggests we might be able to fit a curve to this data of the form\n",
"$$f(t) = c_0 + c_1 t + c_2 \\sin(2 \\pi t) + c_3 \\cos(2 \\pi t),$$\n",
"where $c_0$, $c_1$, $c_2$, and $c_3$ are real constants. The collection of all functions of this form is a vector space, we will denote by $W$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 1. Converting coefficients to a function:\n",
"\n",
"The space $W$ has a natural coordinate system in ${\\mathbb R}^4$ given by coefficients. That is, $c \\in {\\mathbb R}^4$ corresponds to the function \n",
"$$f(t) = c_0 + c_1 t + c_2 \\sin(2 \\pi t) + c_3 \\cos(2 \\pi t).$$\n",
"\n",
"Write a function `to_function(c)` which takes as input a numpy array `c` consisting of four coefficients and returns the function $f$ above. Use currying to return a function. (Currying was discussed in the last example in the [Plotting Graphs notebook](http://wphooper.com/teaching/2020-spring-computation/nb/Plotting%20Graphs.html).) Be sure to use `cos` and `sin` from `numpy` so that the functions work with `numpy` arrays."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Test 1: This should produce a line.\n",
"x = np.linspace(2010,2021,1000)\n",
"f = to_function(np.array([1,1,0,0]))\n",
"plt.plot(x,f(x))\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Test 2: This should produce a wave with a period of 1 year.\n",
"x = np.linspace(2010,2021,1000)\n",
"f = to_function(np.array([0,0,1,-1]))\n",
"plt.plot(x,f(x))\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 2. Finding an orthogonal basis"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"A standard basis for $W$ given as coefficients consists of the following four vectors."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"coefs_0 = np.array([1,0,0,0])\n",
"coefs_1 = np.array([0,1,0,0])\n",
"coefs_2 = np.array([0,0,1,0])\n",
"coefs_3 = np.array([0,0,0,1])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Each coefficient vector represents a function, which can be evaluated on times in our data set. This gives alternate coordinates for $W$ in ${\\mathbb R}^{123}$. We can \"pull back\" the dot product on ${\\mathbb R}^{123}$ to an inner product on $W$. The following is this pull back inner product:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def ip(coefs_A, coefs_B):\n",
" return to_function(coefs_A)(months) @ to_function(coefs_B)(months)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# This is another test of the to_function from problem 1.\n",
"assert abs(ip(coefs_0, coefs_1)-247860.375)<10**-8, \"There is something wrong with to_function.\""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Use the Gram-Schmidt process to find an othogonal basis for $W$ with respect to the inner product given in the `ip` function above. Store the four elements of the orthogonal basis in the variables `orth_0`, `orth_1`, `orth_2`, and `orth_3`."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"orth_0 = "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"orth_1 = "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"orth_2 = "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"orth_3 = "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Test: The vectors should should each be coefficient vectors, with 4 entries:\n",
"assert len(orth_0)==4, \"The length of orth_0 should be 4.\"\n",
"assert len(orth_1)==4, \"The length of orth_1 should be 4.\"\n",
"assert len(orth_2)==4, \"The length of orth_2 should be 4.\"\n",
"assert len(orth_3)==4, \"The length of orth_3 should be 4.\""
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Test: The vectors should be orthogonal. Here we compute the square of the \n",
"# cosine of the angle between two vectors.\n",
"assert ip(orth_0,orth_1)**2 / ip(orth_0,orth_0) / ip(orth_1,orth_1) < 10**-15, \\\n",
" \"Vectors orth_0 and orth_1 are not orthogonal.\"\n",
"assert ip(orth_0,orth_2)**2 / ip(orth_0,orth_0) / ip(orth_2,orth_2) < 10**-15, \\\n",
" \"Vectors orth_0 and orth_2 are not orthogonal.\"\n",
"assert ip(orth_0,orth_3)**2 / ip(orth_0,orth_0) / ip(orth_3,orth_3) < 10**-15, \\\n",
" \"Vectors orth_0 and orth_3 are not orthogonal.\"\n",
"assert ip(orth_1,orth_2)**2 / ip(orth_1,orth_1) / ip(orth_2,orth_2) < 10**-15, \\\n",
" \"Vectors orth_1 and orth_2 are not orthogonal.\"\n",
"assert ip(orth_1,orth_3)**2 / ip(orth_1,orth_1) / ip(orth_3,orth_3) < 10**-15, \\\n",
" \"Vectors orth_1 and orth_3 are not orthogonal.\"\n",
"assert ip(orth_2,orth_3)**2 / ip(orth_2,orth_2) / ip(orth_3,orth_3) < 10**-15, \\\n",
" \"Vectors orth_2 and orth_3 are not orthogonal.\""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 3. Project the data to the space W\n",
"\n",
"Recall ${\\mathbb R}^{123}$ represents the space of all possible measurements that could have been taken at the times in our data set. Our data set `measurements` of $CO_2$ measurements represents a single point in this space. Find the orthogonal projection of `measurements` to $W$ as a list of four coefficients. Store the coefficients as a numpy array in the variable `proj`."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We remark that the elements of ${\\mathbb R}^{123}$ associated to your orthogonal basis can be computed with:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"orth_img_0 = to_function(orth_0)(months)\n",
"orth_img_1 = to_function(orth_1)(months)\n",
"orth_img_2 = to_function(orth_2)(months)\n",
"orth_img_3 = to_function(orth_3)(months)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Also, in ${\\mathbb R}^{123}$, we are using the usual dot product. Compute `proj` below."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"proj = "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Instead of a test, you can visually check if your coefficients track the data. Here is a plot that draws both the data points in blue and the function corresponding to `proj` in red. This gives a least squares approximation."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"proj_f = to_function(proj)\n",
"t = np.linspace(2010,2021,1000)\n",
"plt.plot(t, proj_f(t), \"r\")\n",
"plt.plot(months, measurements, \".\")\n",
"plt.xlabel(\"time in years\")\n",
"plt.ylabel(\"average $C0_2$ measurement in ppm\")\n",
"plt.title(\"Monthly C02 measurements\")\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As a remark, the linear growth term is the most troubling. Your computation suggests that the $CO_2$ is increasing at the following rate (measured in ppm/year):"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"proj[1]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Also, you might guess that we will hit the `milestone` of 500 ppm of $CO_2$ in the atmosphere around the following year:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(500 - proj[0]) / proj[1]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If you are interested in learning more, a more serious discussion can be read here:\n",
"[https://e360.yale.edu/features/how-the-world-passed-a-carbon-threshold-400ppm-and-why-it-matters](https://e360.yale.edu/features/how-the-world-passed-a-carbon-threshold-400ppm-and-why-it-matters)."
]
}
],
"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.7.3"
}
},
"nbformat": 4,
"nbformat_minor": 2
}