CoCalc -- PlotsAndApplications.ipynb

| Hosted by CoCalc | Download

Kernel: Python 3 (system-wide)

Basic plotting using matplotlib

First we'll discuss basic plotting of functions in python.

Then we'll go on to some applications of MAT 105 to advanced material in the CS major.

How do you do I/O in python?
How do you define a function in python?
How do you plot a function in python?

I/O in python

Input is how you get information into a program. Output is how information comes back into the real world.

One way to get information into a program is simply to hard code it like this:

x = 17

The simplest way to get information out of a program in python is to use the print() command.

Let's warm up by doing these two things together.

In [1]:

x = 17
print(x)

17

Advanced input

It's useful to be able to hard code constants into a program, but there are other kinds of input.

User input
File input
Command line input

Let's get some user input in the next cell.

In [2]:

x = input("Type a number:")
print(x)

Type a number:

5

Advanced output

Just printing numbers to the screen is kind of boring. We can make fancier output by using "formatting".

This lets your program incorporate variable values into output statements.

In [3]:

x = input("Type a number")
y = input("Type an animal")
z = input("Type a location")

print("I want {} {}'s delivered to {} immediately.".format(x,y,z))

Type a number

Type an animal

Type a location

I want 6 zebra's delivered to zanzibar immediately.

Representing functions

Now that we can do a little bit of I/O let's work on printing some functions.

But how do you represent a function in python?

There are lots of ways.

For plotting it's useful to represent a function as

A list of numbers in the domain
A list of corresponding numbers in the range

Remember that the domain of a function is the set of inputs and the range is the set of outputs.

Let's do a quick example.

In [4]:

dom = [0,1,2,3,4,5]
ran = [0,1,4,9,16,25]

The above code describes a little bit of the function $f(x)=x^2$ .

You might notice that it's only defined at a few numbers, whereas $f(x)=x^2$ really has $(-\infty,\infty)$ as its domain.

But computers are finite machines and you always have to be satisfied with some finite description of any mathematical object.

Later we'll see how to add more details.

On to plotting!

The first think to do is to import the plotting library matplotlib.

In [5]:

import matplotlib.pyplot as plt

Most plotting in python is done through the pyplot object.

Here we'll be calling it plt for short.

It can do all kinds of things.

Let's make a simple plot of $f(x) = x^2$ .

In [6]:

plt.plot(dom,ran)
plt.show()

That was easy, right?

We just made a crude graph of $y=x^2$ .

Over the next few cells we'll improve the plot by adding

a title
a legend
enriching the domain and range

First let's add a title.

In [7]:

plt.plot(dom,ran)
plt.title(r"A plot of the function $y=x^2$")
plt.show()

Again, that was very easy. We were even a little fancy because we used $\LaTeX$ inside our title.

That's why we had to use "r" and the dollar signs.

If you don't care about how the plot looks you could leave those parts out.

In [8]:

plt.plot(dom,ran)
plt.title("A plot of the function y=x^2")
plt.show()

Adding a legend

Now let's see how to add a legend to the plot.

A legend isn't very interesting unless there is more than one curve.

So at the same time we add the legend, we'll also add a plot of $y=x^3$ .

In [9]:

sq_dom = [0,1,2,3,4,5]
sq_ran = [0,1,4,9,16,25]

cu_dom = [0,1,2,3,4,5]
cu_ran = [0,1,2**3,3**3,4**3,5**3]

plt.plot(sq_dom,sq_ran,label=r"$y=x^2$")
plt.plot(cu_dom,cu_ran,label=r"$y=x^3$")
plt.title(r"$y=x^2$ vs $y=x^3$" )
plt.legend()

plt.show()

That looks pretty nice!

By the way, if you'd rather plot points than lines, you can easily do that.

Here's how. We just change plot to scatter in the code we had before.

Scatterplots can be very useful for representing data.

In [10]:

sq_dom = [0,1,2,3,4,5]
sq_ran = [0,1,4,9,16,25]

cu_dom = [0,1,2,3,4,5]
cu_ran = [0,1,2**3,3**3,4**3,5**3]

plt.scatter(sq_dom,sq_ran,label=r"$y=x^2$")
plt.scatter(cu_dom,cu_ran,label=r"$y=x^3$")
plt.title(r"$y=x^2$ vs $y=x^3$ as a scatterplot" )
plt.legend()

plt.show()

Smoothing out the curve

That was fun, but the curve is a little rough looking. It's also a lot of work to type out the domain and the range.

Is there a better way to do that?

Yes, definitely. It uses a powerful library called "numpy" which is short for Numerical Python.

In [11]:

import numpy as np

Getting the domain

Like we said before, the domain of any function is finite inside a computer.

But we should be able to still get a lot of resolution.

Numpy let's us do that.

One way is to use the arange command which returns a list of numbers between two given values.

Below we're going to get a list of numbers between 0 and 5 which are each separated by a distance of 0.1.

In [12]:

start = 0
stop = 5
step = 0.1
dom = np.arange(start,stop,step)
print(dom)

[0.  0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.  1.1 1.2 1.3 1.4 1.5 1.6 1.7
 1.8 1.9 2.  2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3.  3.1 3.2 3.3 3.4 3.5
 3.6 3.7 3.8 3.9 4.  4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8 4.9]

This will greatly increase the resolution of our plot. We could smooth things out even more by making step smaller.

This improves the domain, but what about the range? Do we have to manually square all of these numbers?

No!

In [13]:

ran = dom**2
print(ran)

[0.000e+00 1.000e-02 4.000e-02 9.000e-02 1.600e-01 2.500e-01 3.600e-01
900e-01 6.400e-01 8.100e-01 1.000e+00 1.210e+00 1.440e+00 1.690e+00
960e+00 2.250e+00 2.560e+00 2.890e+00 3.240e+00 3.610e+00 4.000e+00
410e+00 4.840e+00 5.290e+00 5.760e+00 6.250e+00 6.760e+00 7.290e+00
840e+00 8.410e+00 9.000e+00 9.610e+00 1.024e+01 1.089e+01 1.156e+01
225e+01 1.296e+01 1.369e+01 1.444e+01 1.521e+01 1.600e+01 1.681e+01
764e+01 1.849e+01 1.936e+01 2.025e+01 2.116e+01 2.209e+01 2.304e+01
401e+01]

To define the range, we used a trick called broadcasting.

This means you apparently do something to a list, but it actually gets done to each member of the list.

By "squaring" the domain, we actually produced a list of each domain element, squared.

We could have cubed the elements the same way.

Let's revisit our plot.

In [14]:

sq_dom = np.arange(start,stop,step)
sq_ran = sq_dom**2

cu_dom = np.arange(start,stop,step)
cu_ran = cu_dom**3

plt.plot(sq_dom,sq_ran,label=r"$y=x^2$")
plt.plot(cu_dom,cu_ran,label=r"$y=x^3$")
plt.title(r"$y=x^2$ vs $y=x^3$" )
plt.legend()

plt.show()

That's starting to look pretty professional!

Polynomials

What if you want to plot something more complex, like a quadratic polynomial?

Broadcasting works with all kinds of arithmetic operations.

Let's plot a function of the form $y=ax^2 + bx +c$ .

This time instead of writing "dom" I'm going to use the name x for the same thing.

Python doesn't care what you name your variables -- I'm just doing it for convenience.

Similarly I'll now use y for the range.

In [15]:

a,b,c = 2,2,3
x = np.arange(-10,10,0.1)
y = a*x**2 + b*x + c

plt.plot(x,y)
plt.title(r"$y={}x^2+{}x+{}$".format(a,b,c))
plt.show()

Notice that the line

y = a*x**2 + b*x + c

seems to be defining a single value y, but it's actually working with lists through the magic of broadcasting.

Adding user input

It might be fun to come back to our old idea of getting input from the user.

Then you can easily let the user decide which quadratic to plot.

Because user input is always understood to be a "string" (aka words) we need to convert to type "float" so that the system understands the input as numbers, not words.

In [16]:

a = float(input("Enter a"))
b = float(input("Enter b"))
c = float(input("Enter c"))

x = np.arange(-10,10,0.1)
y = a*x**2 + b*x + c

plt.plot(x,y)
plt.title(r"$y={}x^2+{}x+{}$".format(a,b,c))
plt.show()

Enter a

Enter b

Enter c

Making a real application

You might know from physics class that if you throw an object up or down in a vacuum then the height of the object $t$ seconds later is given by the function

s(t) = -\frac{9.8}{2}t^2+v_0t + s_0

where $v_0$ is the initial velocity (in m/s) and $s_0$ is the initial height (in m).

Let's make a short program that gets $v_0$ and $s_0$ from the user, and shows them how the object will move.

In the code I use the quadratic formula to find when the object lands.

We also use some new options, like showing a grid and annotating the $x$ and $y$ axes.

In [17]:

v0 = float(input("Enter initial velocity "))
s0 = float(input("Enter initial height"))

root1 = (-v0 + np.sqrt(v0**2+2*9.8*s0))/(-9.8)
root2 = (-v0 - np.sqrt(v0**2+2*9.8*s0))/(-9.8)
start = min(root1,root2)
stop = max(root1,root2)
step = (stop-start)/200

t = np.arange(start,stop,step)
y = -9.8/2*t**2+v0*t+s0

plt.plot(t,y)
plt.xlabel("time in seconds")
plt.ylabel("height in meters")
plt.grid(True, which='both')
plt.title(r"$v_0 = {}, s_0 = {}$".format(v0,s0))
plt.show()
print("The object will land at time t= {}.".format(max(root1,root2)))

Enter initial velocity 

Enter initial height

The object will land at time t= 4.050085429146913.

Another application: compound interest

You might know that if you invest $P$ dollars at an annual interest rate of $r$ percent, and the investment compounds $n$ times per year, then the value of the investment after $t$ years is given by the function

V(t) = P (1+\frac{r}{n})^{nt}

Let's make a simple app that helps an investor visualize their returns.

In [51]:

P = float(input("Please enter the initial investment "))
r = float(input("Please enter the annual interest rate "))
n = float(input("Please enter the number of compoundings per year "))
y = float(input("Please enter the number of years you want to invest "))
start = 0
stop = y
step = y/200

t = np.arange(start,stop,step)
v = P*(1+r/n)**(n*t)

plt.plot(t,v)
plt.title("The value of your investment")
plt.xlabel("time")
plt.ylabel("value")
plt.grid(True, which='both')

plt.show()

print("The value of the investment will be ${0:.2f}.".format(P*(1+r/n)**(n*y)))

Please enter the initial investment 

Please enter the annual interest rate 

Please enter the number of compoundings per year 

Please enter the number of years you want to invest 

The value of the investment will be $3310837.85.

Stock market investing

In some cases you may not know the average annual rate of return.

Let's do a version that shows the user a range of possibilities, depending on how the market does.

In [52]:

P = float(input("Please enter the initial investment "))
n = 1
y = float(input("Please enter the number of years you want to invest "))
start = 0
stop = y
step = y/200

t = np.arange(start,stop,step)
r=0.03
v_0 = P*(1+r/n)**(n*t)
r=0.06
v_1 = P*(1+r/n)**(n*t)
r=0.1
v_2 = P*(1+r/n)**(n*t)
r=0.15
v_3 = P*(1+r/n)**(n*t)
r=0.2
v_4 = P*(1+r/n)**(n*t)


plt.plot(t,v_0,label=r"$r=3$%")
plt.plot(t,v_1,label=r"$r=6$%")

plt.plot(t,v_2,label=r"$r=10$%")

plt.plot(t,v_3,label=r"$r=15$%")

plt.plot(t,v_4,label=r"$r=20$%")

plt.title("The value of your investment")
plt.xlabel("time")
plt.ylabel("value")
plt.grid(True, which='both')
plt.legend()

plt.show()

Please enter the initial investment 

Please enter the number of years you want to invest 

Visualizing complex arithmetic

Earlier in the semester you learned about complex numbers.

Each complex number has a real part and an imaginary part.

Python supports complex numbers. You use "j" to indicate the complex part of the number.

In [20]:

c1 = 12-0.3j

c2 = 0.1+1.0j

As you know you can add, subtract, multiply and divide complex numbers.

Let's try that out in python.

In [21]:

print("{} + {} = {}".format(c1,c2,c1+c2))
print("{} - {} = {}".format(c1,c2,c1-c2))
print("{} * {} = {}".format(c1,c2,c1*c2))
print("{} / {} = {}".format(c1,c2,c1/c2))

(12-0.3j) + (0.1+1j) = (12.1+0.7j)
(12-0.3j) - (0.1+1j) = (11.9-1.3j)
(12-0.3j) * (0.1+1j) = (1.5000000000000002+11.97j)
(12-0.3j) / (0.1+1j) = (0.8910891089108912-11.91089108910891j)

The complex plane

Any complex number can be thought of as a point in the complex plane.

Let's write a simple script to plot a complex number.

In [27]:

re = float(input("Enter a real part for a complex number: "))
im = float(input("Enter an imaginary part for a complex number: "))


plt.scatter([re],[im],label=r"${}+{}i$".format(re,im))
plt.scatter([0],[0],label="origin")

plt.xlabel("Real axis")
plt.ylabel("Imaginary axis")
plt.grid(True, which='both')
plt.ylim(min(im,-im)*1.2,max(im,-im)*1.2)
plt.xlim(min(re,-re)*1.2,max(re,-re)*1.2)
plt.arrow(0,0,re*.95,im*.95,head_width=0.08*abs(max(im,re)), head_length=0.08*abs(max(im,re)), fc='black', ec='black')
plt.title("Plotting a complex number")
plt.legend()

plt.show()

Enter a real part for a complex number: 

Enter an imaginary part for a complex number: 

Visualizing complex arithmetic

You might know that arithmetic on complex numbers has a geometric meaning.

Let's make a graphic that illustrates this.

First we'll do addition.

In [36]:

c1 = 2+3j
c2 = 5+1j

scale = max(np.abs(c1),np.abs(c2),np.abs(c1+c2))

plt.xlabel("Real axis")
plt.ylabel("Imaginary axis")
plt.grid(True, which='both')
plt.ylim(-scale,scale)
plt.xlim(-scale,scale)
plt.scatter([0],[0],label="origin")


re, im = c1.real,c1.imag

plt.arrow(0,0,re*.95,im*.95,head_width=0.08*abs(max(im,re)), head_length=0.08*abs(max(im,re)), fc='black', ec='black')
plt.scatter([re],[im],label=r"${}+{}i$".format(re,im))

re, im = c2.real,c2.imag

plt.arrow(0,0,re*.95,im*.95,head_width=0.08*abs(max(im,re)), head_length=0.08*abs(max(im,re)), fc='black', ec='black')
plt.scatter([re],[im],label=r"${}+{}i$".format(re,im))


c3 = c1 + c2


re, im = c3.real,c3.imag

plt.arrow(c1.real,c1.imag,(re-c1.real)*.95,(im-c1.imag)*.95,head_width=0.08*abs(max(im,re)), head_length=0.08*abs(max(im,re)), fc='black', ec='black')
plt.arrow(c2.real,c2.imag,(re-c2.real)*.95,(im-c2.imag)*.95,head_width=0.08*abs(max(im,re)), head_length=0.08*abs(max(im,re)), fc='black', ec='black')
plt.scatter([re],[im],label=r"${}+{}i$".format(re,im))


plt.title("Plotting the sum of complex numbers")
plt.legend()

plt.show()

Complex multiplication

Now we'll visualize complex multiplication. Note that the length of the product is the product of the lengths, and the argument of the product is the sum of the arguments.

In [38]:

c1 = 2+3j
c2 = 5+1j

scale = max(np.abs(c1),np.abs(c2),np.abs(c1*c2))

plt.xlabel("Real axis")
plt.ylabel("Imaginary axis")
plt.grid(True, which='both')
plt.ylim(-scale,scale)
plt.xlim(-scale,scale)
plt.scatter([0],[0],label="origin")


re, im = c1.real,c1.imag

plt.arrow(0,0,re*.95,im*.95,head_width=0.08*abs(max(im,re)), head_length=0.08*abs(max(im,re)), fc='black', ec='black')
plt.scatter([re],[im],label=r"${}+{}i$".format(re,im))

re, im = c2.real,c2.imag

plt.arrow(0,0,re*.95,im*.95,head_width=0.08*abs(max(im,re)), head_length=0.08*abs(max(im,re)), fc='black', ec='black')
plt.scatter([re],[im],label=r"${}+{}i$".format(re,im))


c3 = c1*c2


re, im = c3.real,c3.imag

plt.arrow(0,0,re*.95,im*.95,head_width=0.08*abs(max(im,re)), head_length=0.08*abs(max(im,re)), fc='black', ec='black')
plt.scatter([re],[im],label=r"${}+{}i$".format(re,im))


plt.title("Plotting the product of complex numbers")
plt.legend()

plt.show()

Complex exponents

Because of the way complex products work, exponentiating a complex number causes it to "spiral"...

In [50]:

c1 = 2+3j

scale = max(np.abs(c1),np.abs(c2),np.abs(c1**2.6))

plt.xlabel("Real axis")
plt.ylabel("Imaginary axis")
plt.grid(True, which='both')
plt.ylim(-scale,scale)
plt.xlim(-scale,scale)
plt.scatter([0],[0],label="origin")


re, im = c1.real,c1.imag

plt.arrow(0,0,re*.95,im*.95,head_width=0.08*abs(max(im,re)), head_length=0.08*abs(max(im,re)), fc='black', ec='black')
plt.scatter([re],[im],label=r"power = 1".format(re,im))

c2 = c1**1.2

re, im = c2.real,c2.imag

plt.arrow(0,0,re*.95,im*.95,head_width=0.08*abs(max(im,re)), head_length=0.08*abs(max(im,re)), fc='black', ec='black')
plt.scatter([re],[im],label="power = 1.2".format(re,im))


c3 = c1**1.4


re, im = c3.real,c3.imag

plt.arrow(0,0,re*.95,im*.95,head_width=0.08*abs(max(im,re)), head_length=0.08*abs(max(im,re)), fc='black', ec='black')
plt.scatter([re],[im],label="power = 1.4".format(re,im))


c4 = c1**1.8

re, im = c4.real,c4.imag

plt.arrow(0,0,re*.95,im*.95,head_width=0.08*abs(max(im,re)), head_length=0.08*abs(max(im,re)), fc='black', ec='black')
plt.scatter([re],[im],label="power = 1.8".format(re,im))

c5 = c1**2.2

re, im = c5.real,c5.imag

plt.arrow(0,0,re*.95,im*.95,head_width=0.08*abs(max(im,re)), head_length=0.08*abs(max(im,re)), fc='black', ec='black')
plt.scatter([re],[im],label="power = 2.2".format(re,im))

c6 = c1**2.5

re, im = c6.real,c6.imag

plt.arrow(0,0,re*.95,im*.95,head_width=0.08*abs(max(im,re)), head_length=0.08*abs(max(im,re)), fc='black', ec='black')
plt.scatter([re],[im],label="power = 2.5".format(re,im))



plt.title("Plotting the powers of complex numbers")
plt.legend()

plt.show()

In [0]: