Calculus 2, Lab 0: Numpy and Sympy for Calculus

Calculus 2, Lab 0: Numpy and Sympy for Calculus#

Two useful packages for doing calculus in Python are Numpy (which allows for numerical computations) and Sympy (which supports symbolic manipulation). These can be used without any prior experience programming in Python (or at all), especially in the Jupyter notebook format, where code snippets can be provided for you to edit and run as needed.

A Jupyter notebook (such as this one) has two types of content: code cells, which contain Python code, and Markdown cells, which contain text that is formatted using Markdown.

To edit any text in a Jupyter notebook, double-click on it. To create a new cell, click the \(+\) button in the toolbar above. By default, a new cell is a code cell. To change it to Markdown, use the drop-down menu in the toolbar.

To save (and render) a Markdown cell, or to execute a code cell, either hit Shift+Enter on the keyboard, or click the Run button in the toolbar. For help with Markdown formatting, you may find this markdown cheat-sheet useful. You can type mathematics in a Markdown cell using LaTeX syntax. To enter mathematics, wrap it in $ delimiters. For example, $\sin(4\pi/3)=-\frac{\sqrt{3}}{2}$ produces \(\sin(4\pi/3)=-\frac{\sqrt{3}}{2}\). If we use $$ delimiters, we get “display math” on its own line:

\[ \sin(4\pi/3)=-\frac{\sqrt{3}}{2} \]

A great resource for using Python in mathematics is the website of Patrick Walls at UBC. In particular, you will find resources there on how to do root finding, numerical integration, and differential equations.

Numpy#

Often Numpy is used together with Matplotlib: we load these into Python as follows:

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

The import as creates a prefix that allows us to call functions specific to each library. This can be important when multiple Python libraries are in use, because each might have a different definition of a particular function.

Numpy defines all the common mathematical functions you’ll need, such as the sine function. We would call the function \(\sin(x)\) using np.sin (note the np prefix). To calculate the sine of the angle \(\pi/3\), try entering np.sin(pi/3) in the code cell below, and then running the cell.

Did that not work? It’s because pi is not a core object in Python! It’s defined in Numpy though. Try replacing pi with np.pi and re-running the cell above.

Now let’s try plotting. First we have to define a function to plot.

def f(x):
    return np.sin(x)

A faster (but maybe more confusing) way to define our function is to use Python’s lambda syntax. We also could have done

f = lambda x : np.sin(x)

Plotting takes more than just defining the function – we also have to define the input and output variables!

The linspace command in Numpy produces a set of evenly spaced values. In other words, np.linspace(a,b,n) creates a uniform partition of the interval \([a, b]\) into \(n\) subintervals. The result is a NumPy array that stores the partition points. Here’s how we set 100 points between \(-\pi\) and \(\pi\):

x = np.linspace(-np.pi,np.pi,100)

Next, we can define \(y\) as a function of \(x\), and plot.

y = f(x)
plt.plot(x,y)

Numpy has a lot of powerful numerical tools for doing linear algebra, but since Math 1410 is not officially a prerequisite for Math 2565, we won’t describe those here.

Exercise 1#

Plot \(\tan(x)\) for \(x \in [-\pi, \pi]\).

Does your graph look incorrect? It should! Remember that NumPy is a numerical library. It’s creating the graph by plotting the points \((x_i,y_i)\), where \(x_i\) is a point in the linspace partition, and \(y_i=f(x_i)\). It has no way of determining the existence of a vertical asymptote! If \(x_i\) is close to the asymptote, you’ll see a spike, due to the corresponding large value of \(y_i\). But if you happened to end up with \(x_i\) on the asymptote, you’d get a division by zero error!

Newton’s Method#

Recall Newton’s method for finding roots (zeros) of functions:

Initial data:

  • a function \(f(x)\) on an interval \([a, b]\)

  • an initial guess \(x_0\in [a, b]\)

Given this data, we construct the tangent line to \(y=f(x)\) at the point \((x_0,f(x_0))\):

\[y=f(x_0)+f'(x_0)(x-x_0)\]

We want the point where this line meets the \(x\) axis. Setting \(y=0\), we find

\[x = x_0-\frac{f(x_0)}{f'(x_0)}.\]

The first step in Newton's method: the tangent at (x_0,y_0) meets the x axis at (x_1,0).

Note that it is very important that \(x_0\) is not a critical point! (Why?)

We call the value we just found \(x_1\), and then repeat. In most cases, Newton’s Method converges quickly to a root of our function. (In some cases, a poor choice of \(x_0\) can lead to results getting further and further away.)

Exercise 2#

We are going to develop an algorithm to implement Newton’s method. This will work with core Python, even if you do not import the NumPy library.

For now, we will keep things simple, and use a pre-determined number of steps. Recall that to run a for loop in Python with n steps, we use the syntax for i in range(n):. We will also hard-code the function and its derivative. Later, we will learn how to use SymPy to compute derivatives (and, if necessary, convert that SymPy derivative into a NumPy function).

To avoid conflict with the work we did above, we will call our function g(x). If you want to use f(x), you can, but beware that if you re-run any of the earlier cells, they will use the new definition of f(x).

Problem: find a root of \(g(x)=x^3+x-5\) on the interval \([0,2]\), using initial guess \(x_0=1\).

Aside: why do you know there must be a root on this interval?

Here is your initial setup:

def g(x):
    return x**3+x-5
def gprime(x):
    return 3*x**2+1
xval=[1]

Note that values is an array, containing a single element, which is our initial guess, \(x_0\):

print(values[0])

The reason we use an array is that in a for loop with index i, Python will not understand a variable like x_i. But if we make x into a list, then we can reference the list entries: x[i] is understood.

In the cell below, create a for loop that will create a list containing the next 6 values determined by Newton’s method. Some hints:

  • begin with the line for i in range(6): (and keep the indentation that will be automatically inserted)

  • recall that you can add the next entry in the xval array using the syntax xval.append(), where the parentheses should contain a formula used to determine the next value

  • the formula you want is the Newton’s method formula given above!

  • your last line should be print(values)

  • if you have more than 3 lines, you’ve done too much work

Exercise 3#

This time, we’ll get a little bit fancier. Generally, when using Newton’s method, you want to continue until you’ve found a value for your root to a specified degree of accuracy. The way you know that you have, for example, the first 5 decimal places correct is if the first 5 decimal places do not change when you increment from \(x_i\) to \(x_{i+1}\).

Since we generally do not know how many steps this will take, we will use a while loop instead of a for loop.

If we want 5 decimal places correct, then we want \(|x_i - x_{i+1}|<0.000001\) for some \(i\). (Probably \(0.000005\) would do, but let’s play it safe.) We will use NumPy to provide the absolute value: the function is np.abs().

First, some initial setup:

import numpy as np
xval=[0,1]
i=1
tolerance=0.000001

Since we need to compute a difference of two terms, we need to initialize our array with two terms rather than one. The i=1 line initializes our index, and we’ve set the tolerance to 5 decimal places.

The syntax of a while loop looks like:

while <some condition>:
    <some step, or sequence of steps>
    i+=1

The last line adds one to the index; otherwise, your while loop will run forever!

  • The condition you’re going to want will be an inequality comparing the absolute value of \(x_i-x_{i+1}\) to the tolerance

  • The step you’re going to want to execute is the same as it was in Exercise 2: append the next value to the xval array

  • below the while loop, you’ll again want to print the xval array, or possibly the last term in the xval array, which you can access using the syntax xval[-1:]

Problem: construct a while loop that gives the value of the root from Exercise 2, accurate to 7 decimal places.

Note: you will need to adjust the value of tolerance above! Your output will also contain more than 7 decimal places, so adjust your answer accordingly.

Challenge: Python provides a round() function for rounding to a certain number of decimal places. But the output from your while loop will be considered a list, not a floating point number.

Can you find a way to convert your list value to a floating point value so that you can round to the desired number of decimals?

SymPy#

SymPy allows for symbolic manipulation. It’s also really nice to use in a Jupyter notebook, because you can get the output to display as nicely formatted mathematics, using the init_printing() function.

import sympy as sy
sy.init_printing()

A common first step is to tell SymPy that you’re going to use x and y as variables:

x, y = sy.symbols('x y')

When using mathematical functions, we need to make sure we’re using the SymPy version. We can define functions of one or more variables:

f = (x**2)*sy.sin(2*x)
g = x*y*sy.exp(4*x*y)
f,g

Note the use of * for any multiplication, and ** for exponents. We can take derivatives:

sy.diff(f,x)

sy.diff(f,x,x)

sy.diff(g,x)

sy.diff(g,y)

sy.diff(g,x,y)-sy.diff(g,y,x)

Note that for the function \(g(x,y) = xye^{4xy}\) those are partial derivatives, and the fact that that last line is zero is a theorem!

To do indefinite integrals, we just use the integrate command. We need to specify a function and the integration variable.

sy.integrate(f,x)

(Yes, that’s some integration by parts you’re seeing there!)

For definite integrals, we just have to add in the limits of integration. But note that we need to enter the ordered triple (x.0,sy.pi/2). This time, \(\pi\) comes from SymPy, not Numpy:

sy.integrate(f,(x,0,sy.pi/2))

Just for fun: try changing sp.pi to np.pi. What happens to your output?

We can also plot using SymPy. The basic command is quite simple: plot(f).

sy.plot(f)

You can modify the plot in many ways. For example, to specify a domain, we can write something like plot(f,(x,-2,2)). For more advanced features, such as adding labels and colours, see the SymPy documentation.

To plot more than one graph on the same axes, we need to store each graph, and then combine them. When we first enter the plots, we don’t want them to show, so we use the show=False option. The first way to combine the plots is with the append command. The [0] is there because the object plot(f,show=False) contains two entries in a list: f, and show=False, and we want the first of these.

p1 = sy.plot(f,show=False)
p2 = sy.plot(x**2-2*x,show=False)
p1.append(p2[0])
p1

That likely didn’t do anything, because we forgot that we want to “show” the plot. Instead, we do:

p1.show()

Another option in the documentation is the extend command. I’m not sure what the difference is!

p1 = sy.plot(f,show=False)
p2 = sy.plot(x**2-2*x,show=False)
p1.extend(p2)
p1.show()

Exercise: find the critical points of \(f(x) = \dfrac{x^2+3}{x-1}\). (You may find that Sympy’s factor command is useful!)

Exercise: plot the following on a single set of axes.

  1. The function \(f(x) = x^2\ln(x)\)

  2. The derivative \(f'(x)\)

  3. The antiderivative \(F(x) = \int_1^x f(t)\,dt\)

  4. The tangent line to \(f(x)\) when \(x=1\)

Contents