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?

1

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.

2

In [1]:

x = 17 print(x)

3

17

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.

4

In [2]:

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

5

Type a number:

5

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.

6

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

7

Type a number

Type an animal

Type a location

I want 6 zebra's delivered to zanzibar immediately.

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.

8

In [4]:

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

9

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.

10

In [5]:

import matplotlib.pyplot as plt

11

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

12

In [6]:

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

13

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.

14

In [7]:

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

15

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.

16

In [8]:

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

17

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

18

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

19

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.

20

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

21

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.

22

In [11]:

import numpy as np

23

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.

24

In [12]:

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

25

[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!

26

In [13]:

ran = dom**2 print(ran)

27

[0.000e+00 1.000e-02 4.000e-02 9.000e-02 1.600e-01 2.500e-01 3.600e-01
4.900e-01 6.400e-01 8.100e-01 1.000e+00 1.210e+00 1.440e+00 1.690e+00
1.960e+00 2.250e+00 2.560e+00 2.890e+00 3.240e+00 3.610e+00 4.000e+00
4.410e+00 4.840e+00 5.290e+00 5.760e+00 6.250e+00 6.760e+00 7.290e+00
7.840e+00 8.410e+00 9.000e+00 9.610e+00 1.024e+01 1.089e+01 1.156e+01
1.225e+01 1.296e+01 1.369e+01 1.444e+01 1.521e+01 1.600e+01 1.681e+01
1.764e+01 1.849e+01 1.936e+01 2.025e+01 2.116e+01 2.209e+01 2.304e+01
2.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.

28

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

29

That's starting to look pretty professional!

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.

30

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

31

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.

32

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.

33

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

34

Enter a

Enter b

Enter c

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.

35

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

36

Enter initial velocity

Enter initial height

The object will land at time t= 4.050085429146913.

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.

37

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

38

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.

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.

39

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

40

Please enter the initial investment

Please enter the number of years you want to invest

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.

41

In [20]:

c1 = 12-0.3j c2 = 0.1+1.0j

42

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

Let's try that out in python.

43

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

44

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

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.

45

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

46

Enter a real part for a complex number:

Enter an imaginary part for a complex number:

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.

47

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

48

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.

49

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

50

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

51

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

52

In [ ]:

53