{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"ChEn-3170: Computational Methods in Chemical Engineering Fall 2018 UMass Lowell; Prof. V. F. de Almeida **13Nov2018**\n",
"\n",
"# 12. Linear Least-Squares Fitting w/ Fourier Basis Functions\n",
"$ \n",
" \\newcommand{\\Amtrx}{\\boldsymbol{\\mathsf{A}}}\n",
" \\newcommand{\\Bmtrx}{\\boldsymbol{\\mathsf{B}}}\n",
" \\newcommand{\\Mmtrx}{\\boldsymbol{\\mathsf{M}}}\n",
" \\newcommand{\\Imtrx}{\\boldsymbol{\\mathsf{I}}}\n",
" \\newcommand{\\Pmtrx}{\\boldsymbol{\\mathsf{P}}}\n",
" \\newcommand{\\Lmtrx}{\\boldsymbol{\\mathsf{L}}}\n",
" \\newcommand{\\Umtrx}{\\boldsymbol{\\mathsf{U}}}\n",
" \\newcommand{\\Smtrx}{\\boldsymbol{\\mathsf{S}}}\n",
" \\newcommand{\\xvec}{\\boldsymbol{\\mathsf{x}}}\n",
" \\newcommand{\\avec}{\\boldsymbol{\\mathsf{a}}}\n",
" \\newcommand{\\bvec}{\\boldsymbol{\\mathsf{b}}}\n",
" \\newcommand{\\cvec}{\\boldsymbol{\\mathsf{c}}}\n",
" \\newcommand{\\rvec}{\\boldsymbol{\\mathsf{r}}}\n",
" \\newcommand{\\mvec}{\\boldsymbol{\\mathsf{m}}}\n",
" \\newcommand{\\gvec}{\\boldsymbol{\\mathsf{g}}}\n",
" \\newcommand{\\zerovec}{\\boldsymbol{\\mathsf{0}}}\n",
" \\newcommand{\\norm}[1]{\\bigl\\lVert{#1}\\bigr\\rVert}\n",
" \\newcommand{\\transpose}[1]{{#1}^\\top}\n",
" \\DeclareMathOperator{\\rank}{rank}\n",
"$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"---\n",
"## Table of Contents\n",
"* [Introduction](#intro)\n",
"* [Fourier Expansion](#fexp)\n",
"* [Experimental Data](#data)\n",
"* [Linear System](#ls)\n",
"* [LS Data Fitting](#lsdf)\n",
"* [Modal Analysis](#modal)\n",
"* [Power Spectrum](#spectrum)\n",
"* [Interactive Results](#ires)\n",
"---"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Introduction\n",
"The least-squares method with Fourier basis functions is a powerful computational tool for data fitting and data analysis, the course notes OneNote [ChEn-3170-gen-lsq](https://studentuml-my.sharepoint.com/:o:/g/personal/valmor_dealmeida_uml_edu/Ep9qLSMssi1MjWvl2wYlTHkBSd5aUUoo1fIoe5pswIV0vA?e=Mhkd6f) collect elements of the theory. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Fourier Expansion\n",
"\n",
"The Fourier expansion for approximating a function is\n",
"\n",
"\\begin{equation*}\n",
"b(t) = \\sum\\limits_{k=0}^N \\alpha_k\\,\\cos(k\\,\\mu\\,t) + \\beta_k\\,\\sin(k\\,\\mu\\,t)\n",
"\\end{equation*}\n",
"\n",
"it is a particular form of the generic linear combination of functions ${f_j}(t)$\n",
"\n",
"\\begin{equation*}\n",
"b(t) = \\sum\\limits_{j=1}^{2N+1} x_j\\,f_j(t) .\n",
"\\end{equation*}\n",
"\n",
"If we have a set of values of the independent variable $t_i, i=1,\\ldots,m$, the above Fourier expression when applied to every $t_i$ gives\n",
"\n",
"\\begin{equation*}\n",
"\\bvec = \\Amtrx\\,\\xvec\n",
"\\end{equation*}\n",
"\n",
"where $\\Amtrx = \\begin{pmatrix}\n",
"1 & \\cos(\\mu\\,t_1) & \\sin(\\mu\\,t_1) & \\cos(2\\mu\\,t_1) & \\sin(2\\mu\\,t_1) & \\ldots & \\cos(N\\mu\\,t_1) & \\sin(N\\mu\\,t_1) \\\\\n",
"1 & \\cos(\\mu\\,t_2) & \\sin(\\mu\\,t_2) & \\cos(2\\mu\\,t_2) & \\sin(2\\mu\\,t_2) & \\ldots & \\cos(N\\mu\\,t_2) & \\sin(N\\mu\\,t_2) \\\\\n",
"\\vdots & \\vdots & \\vdots & \\vdots & \\vdots & \\vdots & \\vdots & \\vdots \\\\\n",
"1 & \\cos(\\mu\\,t_m) & \\sin(\\mu\\,t_m) & \\cos(2\\mu\\,t_m) & \\sin(2\\mu\\,t_m) & \\ldots & \\cos(N\\mu\\,t_m) & \\sin(N\\mu\\,t_m) \\\\\n",
" \\end{pmatrix}$, \n",
" $\\xvec = \\begin{pmatrix}\n",
" \\alpha_0 \\\\ \n",
" \\alpha_1 \\\\\n",
" \\beta_1 \\\\\n",
" \\vdots \\\\\n",
" \\alpha_N \\\\\n",
" \\beta_N \\\\ \n",
" \\end{pmatrix}$, \n",
"and \n",
"$\\bvec = \\begin{pmatrix}\n",
" b_1 \\\\ \n",
" b_2 \\\\ \n",
" \\vdots \\\\ \n",
" b_m \\\\ \n",
"\\end{pmatrix} $."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Experimental Data\n",
"Data for this lecture is found in the `data/` directory of the course [repository](https://github.com/dpploy/chen-3170/tree/master/notebooks/data). The data is organized in two columns of NO$_2$ mass concentration ($\\mu$g/m$^3$) versus time (h); see `data/urban-no2.dat`."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"'''Function: read experimental data'''\n",
"\n",
"def read_experimental_data( filename ):\n",
" import io # import io module\n",
" finput = open(filename, 'rt') # create file object\n",
"\n",
" import numpy as np\n",
"\n",
" for line in finput:\n",
" \n",
" line = line.strip()\n",
" \n",
" if line[0] == '#': # skip comments in the file\n",
" continue\n",
" \n",
" var_line = line.split(' = ')\n",
" \n",
" if var_line[0] == 'n_pts':\n",
" n_pts = int(var_line[1])\n",
" time_expt = np.zeros(n_pts)\n",
" no2_mass_cc_expt = np.zeros(n_pts)\n",
" idx = 0 # counter\n",
" else:\n",
" data = line.split(' ')\n",
" time_expt[idx] = float(data[0])\n",
" no2_mass_cc_expt[idx] = float(data[1])\n",
" idx += 1\n",
" \n",
" return (n_pts, time_expt, no2_mass_cc_expt)"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"m = 25\n",
"time [h] = [ 0. 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. 16. 17.\n",
" 18. 19. 20. 21. 22. 23. 24.]\n",
"NO2 [ug/m^3] = [110.49 73.72 23.39 17.11 20.31 29.37 74.74 117.02 298.04 348.13\n",
" 294.75 253.78 250.48 239.48 236.52 245.04 286.74 304.78 288.76 247.11\n",
" 216.73 185.78 171.19 171.73 164.05]\n"
]
}
],
"source": [
"'''Read experimental data'''\n",
"\n",
"import numpy as np\n",
"(n_pts, time_expt, no2_mass_cc_expt) = read_experimental_data('data/urban-no2.dat')\n",
" \n",
"print('m = ',n_pts)\n",
"np.set_printoptions(precision=2)\n",
"print('time [h] =',time_expt)\n",
"print('NO2 [ug/m^3] =', no2_mass_cc_expt)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"'''Function: plot experimental data'''\n",
"\n",
"def plot_experimental_data( time_expt, no_mass_cc_vec ):\n",
" \n",
" import matplotlib.pyplot as plt\n",
"\n",
" plt.figure(1, figsize=(7, 7))\n",
"\n",
" plt.plot(time_expt, no_mass_cc_vec,'r*',label='experimental')\n",
" \n",
" plt.xlabel(r'Time [h]',fontsize=18)\n",
" plt.ylabel(r'NO$_2$ [$\\mu$g/m$^3$]',fontsize=18)\n",
" plt.title('Release of NO$_2$ in Urban Area',fontsize=20)\n",
" plt.legend(loc='best',fontsize=12)\n",
" plt.xticks(fontsize=16)\n",
" plt.yticks(fontsize=16)\n",
" plt.grid(True)\n",
" plt.show()\n",
" print('')\n",
" \n",
" return"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"