Calculus 1, Lab 3: Plotting and Approximating!

By Raheem Mir

As the title suggests, in this notebook we are going to be exploring how we can plot functions with SymPy and NumPy.

We’ll also delve into some numerical approximations when we use and implement Newton’s Method in Python!

Getting Started with SymPy

Before we get started, just like with previous notebooks, we first have to “import” (or load in) the SymPy library, to access its functionality:

→ Run the code cell(s) by selecting and hitting shift + enter or by using the play button in the toolbar at the top.

# Make sure to run this cell!
import sympy as sy
sy.init_printing() 

Recap on Importing Python Libraries:

The line: import sympy as sy loads in the SymPy library for us to use. We give it the shorthand name sy, so every time we call something from SymPy, we prefix it with sy.. This is a good way to keep things organized and to prevent confusion, especially when we are working with multiple Python libraries.

The next line: sy.init_printing() calls SymPy’s init_printing() function, which gets the output to be displayed as nicely formatted mathematics.

Now is a good time to tell Sympy that we want to use \(x\) as a variable. This is done by using the symbols() function.

→ Run the code cell below!

x = sy.symbols('x')    # this makes x a symbol

Plotting with SymPy

The way we plot functions using SymPy is actually quite similar to how you would graph something with a TI-84/TI-83 calculator.

We simply use the plot() function, like so:

sy.plot(f, domain, ...)

f represents the function (or functions) we are visualizing, which can be written explicitly or passed in as a variable.
domain refers to the interval we are the plotting the function(s) over, the default is \([-10, 10]\).
... alludes to the many optional arguments we can use to change the plot’s appearance and behaviour.
An example would be the argument show=False, which stops the graph from being displayed until we use the command, .show().

Assigning our plot(s) to a variable is also a good idea, as it’s more efficient when we are manipulating/adjusting an existing plot, or wanting more control over customization.

Now let’s look at a simple example:

Consider \(f(x) = x^2\).

Using SymPy, we can visualize \(f(x)\) over the interval \([-5, 5]\)

→ Run the code cell below!

p = sy.plot(x**2, (x, -5, 5))   # notice how we input the domain

You might notice that we don’t need to type p on the next line under p = sy.plot(x**2, (x, -5, 5)), for the plot to output, like we had to when representing functions with SymPy in the last notebook.

Let’s try another one!

Consider \(f(x) = \cos(2x)\) on the interval \([0, 2\pi]\)

And while we’re at it, let’s change the color of the line!

→ Run the code cell below!

p1 = sy.plot(sy.cos(2*x), (x, 0, sy.pi*2), line_color='purple')   # Remember to use SymPy's built in functions for trig

On to the next curve!

Consider \(f(x) = \ln(x^2)\) on the interval \([-10, 10]\).

Just like with the trigonometric functions, we’ll be using SymPy’s built in function for natural logarithms, sy.ln()

→ Run the code cell below!

p2 = sy.plot(sy.ln(x**2), line_color="green")   # [-10, 10] is the default domain so we don't need to write it

Here’s an interesting one!

Consider \(f(x) = x\sin(3x)\) on the interval \([-4\pi, 4\pi ]\).

→ Run the code cell below!

f = x*sy.sin(3*x)  
p3 = sy.plot(f, (x, -sy.pi*4, sy.pi*4), line_color="red")

Now it’s your turn!

Complete the following exercises! Take a look at the previous examples if you get stuck.

Exercise 1: Using SymPy, plot the function: \(f(x) = \sin(2x) - 2\sin(x)\) on the interval \([-4\pi,4\pi]\).

# Write your code here!

Exercise 2: Using SymPy, plot the function: \(f(x) = x^3 - 3x^2 + 6x + 8\) on the interval \([-2,4]\).

# Write your code here!

Plotting Multiple Functions

SymPy’s plotting functionality also gives us the ability to plot multiple functions together, essentially combining different plots.

Let’s illustrate this with an example, by plotting the following quadratic equations:

\(y = (x+2)^2\) , \(y = (x-2)^2\), and \(y = x^2\) over the interval \([-12, 12]\).

Plotting multiple functions together is straightforward, as all we need to do is add additional arguments in our call to the plot() function.

→ Run the code cell below!

multi_plot = sy.plot((x-2)**2, (x+2)**2, x**2, (x, -12, 12), legend=True)  # adding a legend to our plot

Visualizing a Function and its Derivative

Now that we know how to plot multiple functions together, one exciting application we can explore is visualizaing a function alongside its derivative:

Consider the function \(f(x) = e^{-x^2}\) on the interval \([-3,3]\).

Let’s plot \(f(x)\) and \(f'(x)\) together! We can utilize SymPy’s diff() function (seen in the last notebook) to compute our derivative.

→ Run the code cell below!

f = sy.exp(-x**2)  # representing our function in SymPy
f_prime = f.diff(x) # taking the derivative with respect to x

p4 = sy.plot(f, f_prime, (x, -3, 3), legend=True) 

Now it’s your turn!

Complete the following exercises! Take a look at the previous examples if you get stuck.

Exercise 3: Using SymPy, plot the function: \(f(x) = x^3 - 3x^2\) and \(f'(x)\) on the interval \([0,3]\).

Hint: Calculate \(f'(x)\) using the .diff() function

# Write your code here!

Exercise 4: Using SymPy, plot the function: \(f(x) = \sin(x)\) and \(f'(x)\) on the interval \([2\pi,2\pi]\).

Hint: Calculate \(f'(x)\) using the .diff() function

# Write your code here!

Plotting with NumPy

Remember in our first notebook when we were first introduced to the Python library, NumPy?

So far in this notebook (and the previous one), we have been using the SymPy library, which is great for symbolic manipulations like representing functions, and taking derivatives. NumPy on the other hand is the library of choice when it comes to doing numerical computations. It implements many mathmatical functions, like trigonmetric functions, logarithms, and more.

When NumPy is used in conjunction with Matplotlib, the popular Python plotting library that SymPy’s plotting functionality is based on, we have a powerful framework to create more flexible and customizable plots, that wouldn’t neccessarily be done with SymPy.

Getting Started with NumPy

Just as before, to access the functionality of the libraries NumPy and Matplotlib, we have to “import” or load them in:

→ Run the code cell below!

import numpy as np
import matplotlib.pyplot as plt

Similar to how we gave SymPy the shorthand name sy, when we call functions from NumPy, we use the prefix np, and the shorthand plt when we use functions from Matplotlib.

Plotting Walkthrough

Let’s plot the graph \(y = \tan(x)\) on the interval \([-2\pi, 2\pi]\).

Given its periodic nature and vertical asymptotes at values such as \(x = \dfrac{\pi}{2}\) and \(x = \dfrac{3\pi}{2}\), this will be an interesting curve to explore. Plotting functions with asymptotes can be a challenge for many plotting tools, so let’s see how a numerical plotting approach handles things…

We can use the linspace() function to produce a set of evenly spaced \(x\)-values. In the cell below we are setting \(1000\) points in between \(-2\pi\) and \(2\pi\).

→ Run the code cell below!

x = np.linspace(-2*np.pi, 2*np.pi, 1000)

Now we define our \(y\) by setting it equal to the tan() function from NumPy.
→ Run the code cell below!

y = np.tan(x)

We can use the ylim() function, to set the limits of the \(y\)-axis, just like setting the windows on a TI-84 graphing calculator!
In our second line, we finally get to plotting, calling Matplotlib’s plot() function and passing in our \(x\) and \(y\) values.

→ Run the code cell below to plot \(\tan(x)\)

plt.ylim(-6, 6)  # think of this line as adjusting our windows on a graphing calculator
plt.plot(x, y)

Hmm, is the output what you were expecting? If not, what’s wrong with it? Why might this be the case?

The vertical lines we see in the plot are a result of \(\tan(x)\) having very large positive values to the left of the asymptotes, and very large negative values present immediately to the right of the same asymptotes. Since we are plotting numerically, values near asymptotes aren’t taken into consideration, leading to these peculiar joining lines.

To address this, we can apply a common solution using Numpy, that adjusts for \(y\)-values around the asymptotes, and then plot our updated values.

→ Run the code cell below!

y[:-1][np.diff(y) < 0] = np.nan    # this line helps adjust for values near the asymptotes 
plt.ylim(-6, 6)
plt.plot(x, y)  

Exploring Newton’s Method

There are many occasions where we cannot solve equations by using algebraic methods. In these instances, we can utilize numerical methods to approximate solutions. In the last notebook, we saw the Bisection method for root finding, this time around, we will be working with a different technique, known as Newton’s Method.

Newton’s Method is based on tangent lines. The method starts with an initial guess of where the root of a function \(f(x)\) is, which we’ll call \(x_{0}\). 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)}.\]

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.)

We then repeat the process and see where the tangent line at the point \((x_{1}, f(x))\) meets the \(x\)-axis, which can be called \(x_{2}\). This proccess continues until a certain tolerance level (i.e. accuracy to a specified number of decimal places) or maximum number of iterations is reached.

The formula for Newton’s Method is as follows: $\( x_{n + 1} = x_{n} - \dfrac{f(x_{n})}{f'(x_{n})}, \)\( where \)x_{n + 1}$ is our next iteration.

Now, we could repeatedly do the neccessary computations by hand to approximate a solution for a given equation… or we can use Python to do it all for us!

Let’s translate the formula for Newton’s Method into a Python function!

Implementing Newton’s Method in Python

Defined in the code cell below is a Python function that performs Newton’s Method for approximating solutions, it implements the entire process, from the formula above, to all of the neccessary iterations.

Since it has already been created, we can call it like we would any other function, and get right into using it to solve equations!

We can call it like so:

NewtonsMethod(f, fprime x0, tolerance)

f represents the function we are inputting and trying to solve. Importantly, this is a traditional Python function, made with the def keyword and defined on an input/output basis, not like the symbolic representations with SymPy we’ve been using up to this point.

fprime represents the derivative of the previous function f. It too is a Python function, that we have to create using the def keyword.

x0 represents the initial estimate of the root, our starting point.

tolerance represents the level of accuracy we want (usually a small decimal value like 0.0001), it also partially determines how many iterations our function will perform.

→ Run the code cell below!

def NewtonsMethod(f, fprime, x0, tolerance):
    x = x0
    for iteration in range(0, 100):
        oldx = x
        x = x - f(x) / fprime(x)
        print(f"x{iteration + 1} = {x}")
        if abs(x - oldx) < tolerance:
            return x
        
    return "Did not converge after 100 iterations"

Now let’s do an example!

Let’s approximate the root of \(x^3 - x^2 - 1 = 0\) using Newton’s Method, accurate to three decimal places, with an initial guess of \(x_{0} = 1\).

First, we have to define the two Python functions we’ll need to input into NewtonsMethod():

• f(x) for the equation \(x^3 - x^2 - 1\)

• fprime(x) for the previous equation’s dervative, \(3x^2 - 2x\)

→ Run the code cells below! Notice how we used the def keyword and the format (or syntax) used to define functions.

def f(x):
    return x**3 - x**2 - 1

def fprime(x):
    return 3*x**2 - 2*x 

Now that we’ve defined our inputs, let’s call the NewtonsMethod() function!

→ Run the code cell below!

NewtonsMethod(f, fprime, 1, 0.0001)

Great! NewtonsMethod() took in our inputs, and printed out each iteration of applying Newton’s Method along with the final result!

Let’s try another example!

Let’s approximate the root of \(f(x) = \cos(x)\) using Newton’s Method, accurate to 4 decimal places, with an initial approximation of \(x_{0} = 1.5\).

Just as before, let’s define the Python functions we’ll be passing into NewtonsMethod():

• f2(x) for the function \(f(x) = \cos(x)\)

• f2prime(x) for the previous function’s dervative, \(f'(x) = -sin(x)\)

You’ll notice we’ve added a “2” at the end of the function’s names to prevent confusion and to ensure everything runs smoothly.

→ Run the code cells below!

def f2(x):
    return np.cos(x) # using NumPy's built in cosine function

def f2prime(x):
    return -np.sin(x) # using NumPy's built in sine function

Now that we’ve defined everything we need, let’s call NewtonsMethod() along with our other inputs!

→ Run the code cell below!

NewtonsMethod(f2, f2prime, 1.5, 0.00001)

Since the root of \(\cos(x) = 0\) is known, namely \(\dfrac{\pi}{2}\), how does it compare to the result of our approximation using Newton’s Method?

→ Write your answer in the markdown cell below! (Double-click to edit)

Write Your Answer Here!

Double-click on this cell to edit it and type your answer!

Now it’s your turn!

Complete the following exercises! Take a look at the previous examples if you get stuck.

Exercise 5:
Using Newton’s Method, approximate the solution to \(f(x) = x^5 + x + 1\) to three decimal places, with an initial approximation of \(x_{0} = -1\).

→ Write your answer in the code cells below! The functions you’ll need to pass into NewtonsMethod() are already done for you.

def f3(x):
    return x**5 + x + 1

def f3prime(x):
    return 5*x**4 + 1

# Write the rest of the code here!

Exercise 6:
Using Newton’s Method, approximate the solution to \(f(x) = x^2 -2\) to four decimal places, with an initial approximation of \(x_{0} = 1.5\).

→ Write your answer in the code cells below! The functions you’ll need to pass into NewtonsMethod() are already done for you.

def f4(x):
    return x**2 - 2

def f4prime(x):
    return 2*x 

# Write the rest of the code here!