<jupyter_start><jupyter_text>Lab 6 IntroductionThis lab covers some familiar ground then introduces numerical root finding, matrix manipulation, and curve fitting using NumPy and SciPy.Numerical root finding could be useful for estimating a single critical model parameter.Matrix calculations will be critical for solving systems of differential equations. Being able to handle matrices numerically is also a useful life skill. When your first-year linear algebra is a distant memory and your continued employment is contingent on you finding an eigenvalue, it could help you find your way out of a dark place.Curve fitting will be useful if you want to figure out the parameters of a model given some empirial data.We will cover1. Basic Maths1. Plotting Functions1. Finding Roots of Functions1. Finding Eigenvalues1. Curve FittingHopefully the first two sections are old hat to you. They are included here as a summary. SetupAs always, start with some imports.<jupyter_code>from numpy import linspace, exp, sqrt, cos, sin, pi, array, \
                  random, meshgrid, sqrt, arange, meshgrid
from numpy.linalg import eigvals, eig
from numpy.testing import assert_almost_equal
from scipy.optimize import fsolve, curve_fit
from scipy.integrate import odeint
import seaborn as sns
import pandas as pd<jupyter_output><empty_output><jupyter_text>1. Basic Maths As you know, Python works just fine as a calculator.<jupyter_code>2 + 3
7 - 5
34 * 212
1234 / 5678
3**7
sqrt(2)
cos(pi)
sin(pi)<jupyter_output><empty_output><jupyter_text>Note that the last answer is $1.2246\ldots\times 10^{-16}$, which is very very small and should remind you that non-integer calculations always come with a tiny bit of error.  2. Plotting FunctionsPlotting functions is a very important aspect of mathematics. It gives a very quick and visual way of understanding the behaviour of functions. When we plot a function, say $f(x) = \sin x$, we first must define the range of $x$. Let’s say we want to plot the curve $y = \sin x$ in the range $−5 \leq x \leq 5$. We define the range of x by typing the following command.<jupyter_code>x = linspace(-5, 5, 101)<jupyter_output><empty_output><jupyter_text>This creates a vector `x` whose elements are from −5 to 5 in steps of 0.1. The vector `x` has 101 elements in the vector. (Note that `x` is actually a NumPy array. NumPy calls matrices and vectors "arrays".) You can see what is in `x` like so:<jupyter_code>print(x)
y = sin(x)
data = pd.DataFrame({'x': x, 'y': y})
ax = sns.lineplot(data=data, x='x', y='y')<jupyter_output><empty_output><jupyter_text>Note that you can download the plot by right-clicking on it and saving the image.If you would like fine-grain control of how your plot looks, you can do that using `matplotlib` commands. For instance, `plt.figure(figsize=(5, 5))` sets the shape of the display, which can be useful if you want it to be a square.<jupyter_code>import matplotlib.pyplot as plt
fig = plt.figure(figsize=(5, 5))
ax = sns.lineplot(data=data, x='x', y='y')
ax.set_title('A plot of y=sin(x)')
ax.set_xlabel('x-axis label')
ax.set_ylabel('y-axis label');<jupyter_output><empty_output><jupyter_text>Hence for basic plotting you need to:1. define the range;1. define the function you want to plot;1. put your data in a DataFrame; 1. \[optional\] melt your data frame1. plot the data in your DataFrame. Recall that that Step 3A is only required when you want multiple plots on the same axes. If you set up the DataFrame carefully you can save yourself a bunch of hassle with `matplotlib` trying to get your plot to look pretty.Let’s try plotting another curve $y = x2 − x \sin(3x) − x\mathrm{e}^x$ for $−3 \leq x \leq 3$. Note that when you raise $x$ to the power of 2, multiply it by other arrays, or apply NumPy functions like `sin` and `exp`, NumPy assumes that you want to perform _elementwise_ operations.<jupyter_code>x = linspace(-3, 3, 61)
y = x**2 - x*sin(3*x) - x*exp(x)
data = pd.DataFrame({'x': x, 'y': y})
ax = sns.lineplot(data=data, x='x', y='y')<jupyter_output><empty_output><jupyter_text>3. Finding Roots of FunctionsIn many circumstances one would want to find the value of x when the function $f(x) = 0$, i.e., the roots of $f(x)$. Mathematicians, statisticians, and computer scientists have quite a lot to say about root finding but in many cases you can get away with calling a SciPy function.The first thing you need is a “guess” for the root. A good way of getting an appropriate guess is to simply plot the function and see if you can spot where it intersects the $x$-axis.Say for example, we would like to find a root of the function $f(x) = x−\sin x+1$. So we first plot the curve $y = f(x)$ to see “roughly” where the curve “cuts” the $x$-axis. Remember the four lines of code that will plot this curve?<jupyter_code>x = linspace(-3, 3, 61)
y = x - sin(x) + 1
data = pd.DataFrame({'x': x, 'y': y})
ax = sns.lineplot(data=data, x='x', y='y')<jupyter_output><empty_output><jupyter_text>So from our naked eye, we can guess that the value of the root of the function is around $x = −2$. So now we type in the following commands to get the a more accurate answer.<jupyter_code>def func(x):
    return x - sin(x) + 1

fsolve(func, -2)<jupyter_output><empty_output><jupyter_text>If you want to be fancy, you can do it all on one line using a _lambda_ function. Lambda functions exist in several programming languages and are convenient for defining a function you only intend to use once.<jupyter_code>fsolve(lambda x: x - sin(x) + 1, -2)<jupyter_output><empty_output><jupyter_text>4. Finding EigenvaluesEigenvalues provide very important information regarding the behaviour of various systems. In the next few weeks we will be learning how we can use matrix methods and eigenvalues to solve linear systems of differential equations. Once we go beyond $2 \times 2$ matrices, finding eigenvalues by hand is a pain and using Python easy. How do we find the eigenvalues of matrices using Python? It is disgustingly simple! Read on and find out.Say we would like to find the eigenvalues of matrix A whose elements are\begin{align*}\left(\begin{array}{2}1&2\\3&4\end{array}\right)\end{align*}We first need to define matrix `A` as a NumPy array and we do this by typing the following command<jupyter_code>A = array([[1, 2], [3, 4]])<jupyter_output><empty_output><jupyter_text>which gives<jupyter_code>print(A)<jupyter_output><empty_output><jupyter_text>To obtain the eigenvalues we simply type<jupyter_code>eigvals(A)<jupyter_output><empty_output><jupyter_text>So what? Anyone can find the eigenvalues of a $2 \times 2$ matrix, I hear you say!Let’s now find the eigenvalues of a $3 \times 3$ matrix `B` whose elements are\begin{align*}\left(\begin{array}{2}1&2&3\\4&5&6\\7&8&9\end{array}\right)\end{align*}Imagine having to find the eigenvalues of this $3 \times 3$ matrix by hand. In Python, as before, we need to first define the matrix, and then use the `eigvals` command as follows.<jupyter_code>B = array([[1, 2, 3],
           [4, 5, 6],
           [7, 8, 9]])
eigvals(B)<jupyter_output><empty_output><jupyter_text>Note that one of the eigenvalues is zero. How easy was that! You can do this for any $n \times n$ matrix. Why don’t you make up your own $6 \times 6$ matrix and find the eigenvalues of that matrix. It’s the same two commands you need. One to define the matrix, and then the `eigvals` command. If you would like more practice (and information) on using NumPy arrays, please check out [this lesson](https://swcarpentry.github.io/python-novice-inflammation/02-numpy/index.html) from Software Carpentary. 5. Curve FittingThe following material is just the example from the SciPy documentation that is available [here](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html).Start by defining a model that we want to fit.<jupyter_code>def func(x, a, b, c):
    return a * exp(-b * x) + c<jupyter_output><empty_output><jupyter_text>Define the data to be fit with some noise:<jupyter_code>xdata = linspace(0, 4, 50)
y = func(xdata, 2.5, 1.3, 0.5)
random.seed(1729)
y_noise = 0.2 * random.normal(size=xdata.size)
ydata = y + y_noise

data = pd.DataFrame({'x': xdata, 'data': ydata})
ax = sns.lineplot(data=data, x='x', y='data')<jupyter_output><empty_output><jupyter_text>Fit for the parameters a, b, c of the function `func`:<jupyter_code>popt, pcov = curve_fit(func, xdata, ydata)
print('a, b, c = ', popt)
a, b, c = popt
yfitted = func(xdata, a, b, c)

data = pd.DataFrame({'x': xdata, 'data': ydata, 'fitted': yfitted})
data = data.melt(id_vars=['x'])
ax = sns.lineplot(data=data, x='x', y='value', hue='variable')<jupyter_output><empty_output><jupyter_text>How close were the fitted parameters to the actual parameters?  Exercises  Slope field and DE solution plotPlot on the one figure the slopefield for the DE\begin{align*}\frac{\mathrm{d} y}{\mathrm{d} x} = 2.5y (1 − y)\qquad y(0) = 0.5,\end{align*}and the solutions $y(x)$ with $y(0) = 0.2$, $y(0) = 0.5$ and $y(0) = 0.8$.Start by writing down a new definition for `diff_eq` below. Do not change the function's name or inputs.<jupyter_code>def diff_eq(y, x):
    ### your implementation goes here<jupyter_output><empty_output><jupyter_text>If you have implemented `diff_eq` correctly, the following should print "nice job".<jupyter_code>assert_almost_equal(diff_eq(0.4, 0), 0.6)
assert_almost_equal(diff_eq(0.4, 10), 0.6)
assert_almost_equal(diff_eq(1, 0), 0)
print("nice job")<jupyter_output><empty_output><jupyter_text>`create_slope_field` from Lab 03 is included here for your use.<jupyter_code>def create_slope_field(diff_eq, x, y, args=(), initial_values=()):
    """
    Plot a slope field and numerical solutions to initial value problems.
        
    Parameters
    ----------
    diff_eq : callable(y, x)
        Computes the derivative of y at x.
    x : array
        2-dimensional array of x values to plot slope field.
    y : array
        2-dimensional array of y values to plot slope field.
    args : tuple
        Any additional arguments to diff_eq.
    initial_values : tuple
        Initial values from which to numerically solve diff_eq.
        
    Returns
    -------
    fig : matplotlib.figure
        Resulting matplotlib figure.
    ax : matplotlib.axes
        Resulting matplotlib axes.
    """
    S = diff_eq(y, x, *args)
    L = sqrt(1 + S**2)
    fig, ax = plt.subplots(figsize=(5, 5))
    q = ax.quiver(x, y, 1/L, S/L, scale=25, headwidth=0, headlength=0, color='grey')
    x = linspace(x.min(), x.max())
    data = {'x': x}
    for y0 in initial_values:
        y = odeint(diff_eq, y0, x, args)[:,0]
        data[f'y({x[0]}) = {y0}'] = y
    data = pd.DataFrame(data)
    data = data.melt(id_vars=['x'], value_name='y', var_name='initial value')
    sns.lineplot(data=data, x='x', y='y', hue='initial value')
    return fig, ax<jupyter_output><empty_output><jupyter_text>Look at how `create_slope_field` works and use it to create your plots. This function gives good examples of1. the minimum requirements for comments in your functions (ie. the docstring), and1. use of `odeint`, iteration, and `pandas.DataFrame.melt`.Create the plot for the region $0 < x < 1$ and $0 < y < 1$. To help, the `meshgrid` statement is provided.<jupyter_code>x, y = meshgrid(arange(0, 1, 0.05), arange(0, 1, 0.05))
### use create_slope_field here<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10781552543566274, 'fraction_numerical': 0.018130610807956346, 'mean_word_length': 4.507998061076103, 'pattern_counts': {'":': 0, '<': 85, '<?xml version=': 0, '>': 81, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '20', 'dir': 'jupyter-structured-clean-dedup', 'id': '31057', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 3801, 'n_tokens_neox': 3637, 'n_words': 1473, 'path': 'notebooks/lab-06.ipynb', 'repo_name': 'BenKaehler/mm-labs', 'size': '11362'}