Contact
CoCalc Logo Icon
StoreFeaturesDocsShareSupport News AboutSign UpSign In
| Download

cocalc mathday!

Views: 529

Welcome to the Sage Lab!

Today we'll show you some cool math things you can do using a computer, and teach you how to solve math problems using programming!

Part 1: Sage: The Final Frontier

We can use computers to do arithmetic. The cell below contains a basic sum. Click on the 2+22+2 and hit shift-enter to evaluate that sum (in general we will be executing snippets of code in this worksheet; to do so just click on the cell and hit shift-enter, or hit the green Run button at the top left of the page).

Example

Use shift-enter to evaluate the following cell.

2+2

FUN-ctions

We can define a variable xx and a function ff which operates on the input xx. In the code below we tell the computer that xx is a variable, and ff is the function f(x)=(x+2)(x34)f(x) = (x+2)(x^3-4). We then calculate f(3)f(3).

Example 1

Use shift-enter to evaluate the following cell.

x = var('x') f = (x + 2)*(x^3 - 4) f(x = 3)

Example 2

Alternatively we can define our function exactly the way we would write it out by hand.

f(x) = (x + 2)*(x^3 - 4) f(3)

Cool Pictures

With Sage we can also plot functions with 2D and 3D graphs.

Example

Plot the function we defined above by evaluating the following cell.

plot(f, xmin=-3, xmax=3)

Calculus 1

For those of you who know a bit of calculus, we can compute the derivative of f(x)f(x).

Example

Find the derivative of the function ff we defined above by evaluating the following cell.

f.derivative()

Calculus 2

For those with a bit more calculus, we can use Sage to both approximate and calculate (when possible) integrals.

Example

Evaluate the following cell to find an interactive Riemann sums application.

@interact(auto_update=True) def _( func = input_box(cos(x)*x + 3,label='Function',width=50), rectangle_type = selector(['Midpoint', 'Left', 'Right', 'Upper', 'Lower'], nrows=1, label="Endpoint rule"), N=slider([1..100],label='Number of rectangles',default=3), interval=range_slider([-50..50],default=(-10,10),label="Interval"), line_color=color_selector(widget='colorpicker', label="Line color"), bar_color=color_selector(widget='colorpicker', label="Rectangle color",default='#00b334'), dot_color=color_selector(widget='colorpicker',label='Dot color',default='#000000') ): import time start = time.time() a = interval[0] b = interval[1] dx = (b - a)/N area = 0 f = fast_callable(func,vars=['x']) def height(a,b,f): if rectangle_type == 'Left': c = a elif rectangle_type == 'Midpoint': c = (b+a)/2 elif rectangle_type == 'Right': c = b elif rectangle_type == 'Upper': #could be much better Y = [f(a + (b-a)*i/100.0) for i in range(101)] c = a + (b-a)*Y.index(max(Y))/100.0 elif rectangle_type == 'Lower': Y = [f(a + (b-a)*i/100.0) for i in range(101)] c = a + (b-a)*Y.index(min(Y))/100.0 return c,f(c) rects = Graphics() pts = Graphics() for i in range(N): xi = dx*i + a mx,my = height(xi,xi+dx,f) rects += polygon([(xi,0),(xi+dx,0),(xi+dx,my),(xi,my)],color=bar_color,alpha=0.6) rects += line([[xi,0],[xi,my]],thickness=1,color='#000000',alpha=1) #left side rects += line([[xi+dx,0],[xi+dx,my]],thickness=1,color='#000000',alpha=1) #right side rects += line([[xi,my],[xi+dx,my]],thickness=1,color='#000000',alpha=1) #top pts += point([mx,my],color=dot_color,size=50,alpha=1) area += dx * my p = plot(f,(x,a,b),color=line_color,thickness=2) end = time.time() print('Time taken: %s' % (end-start)) show(p + rects + pts) print("\n") show("$\\text{Estimated area} \\approx %s$" % area.n(100)) try: show("$\\text{Calculated integral} = %s $" % integrate(func,x,a,b)) show("$\qquad\qquad\qquad\qquad \\approx %s$" % integrate(func,x,a,b).n(100)) except Exception as e: print("Unable to integrate")

Even cooler pictures

Why restrict ourselves to 2-dimensional plots?

Example 1

Here's a 3d graph of the function z=ysin(xy)cos(x)z = y\sin(x-y)\cos(x):

x,y = var('x,y') plot3d(sin(x - y)*y*cos(x), (x,-3,3), (y,-3,3), mesh=True)

Example 2

Here's a plot of the famous Riemann zeta function ζ(s)=i=11ns\zeta(s) = \sum_{i=1}^\infty \frac{1}{n^s} that uses complex numbers (ask your presenters about this one)!

complex_plot(zeta, (-20,10), (-30,30))

Example 3

Ever heard of the Mandelbrot Set? It is an example of some really cool mathematical objects called fractals.

Try move some of the sliders and see what happens...

@interact def mandel_plot(expo = slider(-10,10,0.1,2), formula = ['mandel','ff'], iterations=slider(1,100,1,30), zoom_x = range_slider(-2,2,0.01,(-2,1)), zoom_y = range_slider(-2,2,0.01,(-1.5,1.5))): var('z c') f(z,c) = z^expo + c ff_m = fast_callable(f, vars=[z,c], domain=CDF) # messing around with fast_callable for i in range(int(iterations)/3): f(z,c) = f(z,c)^expo+c ff = fast_callable(f, vars=[z,c], domain=CDF) def mandel(z): c = z for i in range(iterations): z = ff_m(z,c) if abs(z) > 2: return z return z print 'z <- z^%s + c' % expo # calling ff three times, otherwise it fast_callable exceeds a recursion limit if formula is 'ff': func = lambda z: ff(ff(ff(z,z),z),z) elif formula is 'mandel': func = mandel complex_plot(func, zoom_x,zoom_y, plot_points=200, dpi=150).show(frame=True, aspect_ratio=1)

Codes and Secrets, i.e., the cool stuff

Let's do some cryptography. Enter a message in the text box below, then use the slider to encrypt it! This is an example of a Caesar cipher, where we take a message and shift all the letters to the left or right by a fixed amount (for example a shift by 1 would send aa to bb, bb to cc etc.).

Example

Can you create a secret message and pass it to the person next to you to decode?

Try to intercept someone elses message and decoding it.

@interact def _( message=input_box(default='sage is cool',type=str,label='Message'), key=slider([0..25],0,label='Shift by:',width=35), method=selector(["Encrypt","Decrypt"],label="Method",buttons=True) ): S = ShiftCryptosystem(AlphabeticStrings()) k = key m = S.encoding(message) print "Message = %s" % m print "Key = %s" % k if method is "Encrypt": print "Encrypted = %s" % S.enciphering(k,m) if method is "Decrypt": print "Decrypted = %s" % S.deciphering(k,m)
############################KEEP####################################################
##############################GOING#################################################
###############################YOU'RE###############################################
##################################DOING#############################################
####################################GREAT###########################################

Part 2: Mathematical Programming

There are thousands of ways we can use computers to help us do math. We're going to take a closer look at how to program a computer to effortlessly do tedious calculations that would take a human ages to do. This will allow us to explore a beautiful unsolved problem in mathematics relating to prime numbers.

A Modest Goal

What is the sum of the first 100100 whole numbers? That is, can you calculate 1+2+3+4+...+99+1001 + 2 + 3 + 4 + ... + 99 + 100?

Example

In the cell below we've typed out a simple command in sage: 1+2+31 + 2 + 3. To evaluate this expression, click on the text and hit shift-enter.

1 + 2 + 3

Challenge

Can you modify the above code to compute the sum of the first 1010 whole numbers, i.e. 1+2+3+4+5+6+7+8+9+101 + 2 + 3 + 4 + 5 + 6 + 7 + 8 + 9 + 10?

Feel free to use copy and paste if you want.

# put code here

A LOOPy solution

We could also modify the above code to sum the first 100100 whole numbers, but this would be pretty boring and tedious. Luckily, we have a much better solution, the for loop. It allows us to specify an instruction that repeats in a systematic way. Intelligent loops are one of the most fundamental programming ideas.

The code below defines a number xx, which we initially set to be zero. Note that here "==" means "set the variable on the left to be the value on the right". Then starting at n=1n=1 we add nn to our number xx. Then we increment nn to be equal to 22, and add it to xx again, then we increment nn to 33 and add it to xx again etc. We keep doing this all the way up to n=10n=10. In so doing we've added the numbers 1,2,3,4,5,6,7,8,91, 2, 3, 4, 5, 6, 7, 8, 9 and 1010 to 00. Finally we print out the number xx.

Hit shift-enter after clicking on the code below, and you should find that the number printed out should equal the number you calculated above.

#Adds the numbers 1, 2, 3, ..., 10 together x = 0 for n in range(1, 11): x = x+n print(x)

Challenge

Copy and paste the code in the following cell into the blank one below. Can you modify it to calculate the sum of the first 100 whole numbers? How about the first 10000 whole numbers?

#insert code here x = 0 for n in range(1, 101): x = x+n print(x)

All in all

See how easy that was? Much faster than doing it in your head or typing out the full instruction manually!

For loops are cool.

Level 2

Above we found out how to add all the numbers 1,2,3,,1001, 2, 3, \dots, 100. What if we only want to add the even numbers, i.e., 2,4,6,8,,96,98,1002, 4, 6, 8, \dots, 96, 98, 100?

Don't Reinvent the Wheel

We can just modify our code to do exactly this. Take a look at the code below.

#Adds the numbers 2, 4, 6, ..., 98, 100 together x = 0 for n in range(1,101): if 2.divides(n): x = x+n print(x)

Explanation

What we've done above is to first test if the number nn is even or odd by checking if it's divisible by 2. If it is, then we add it to our number xx. To do this we've used something called a if statement. This tests a specified condition and only executes the following line or lines of code if that statement is true.

For us, we only add nn to xx if 2 divides nn.

FUN-ctions again, but in computers

Now we're going to use a for loop to do something even more awesome. Remember the definition of a prime number:

A positive whole number pp is prime if the only positive whole numbers that divide pp are 1 and pp itself.

For example 2,3,5,72, 3, 5, 7 and 1111 are all prime, but 66 and 1212 are not (as both 66 and 1212 are divisible by 22 and 33).

Sometimes we want to repeat the same snippet of code over and over again, in this case we made a function. Just above, we used the function divides. We are going to define a function that uses a for loop to test if a given integer is prime or not. Below is some code that creates a function is_it_prime, that takes as input a positive whole number pp and returns the value True if pp is a prime number, or False if it isn't.

def is_it_prime(p): primality = True if 2.divides(p): primality = False return primality

Try it out

Click on the code above and hit shift-enter. This evaluates the code; the computer then learns the function is_it_prime, so we can then use it. In the empty cell below type is_it_prime(2) and hit shift-enter to see if the code works. Does it work for n=3n=3? How about n=4,5,6,7,8n=4, 5, 6, 7, 8 and 99?

#insert code here

Bugs in the code

Oh no, our function doesn't work properly!

You should have found that it returns True for n=9n=9. But 99 isn't prime - it's divisible by 33! It also returns False for n=2n=2 - but 22 is a prime number.

Take a closer look at the function. For input number pp we are only testing if 22 divides pp. If it does then we return False; otherwise return True. However there are plenty of odd numbers that aren't prime - for example 9, 15, 21 and 45. And there is one even number that is prime - namely the number 22.

In order to correctly ascertain if a number pp is prime we need to check if any numbers less than pp divide into pp.

Challenge

Can you see why we don't need to check any numbers larger than pp?

# code here

Debugging

Let's fix this so that it does check primality. First, we only need to check if numbers smaller than the input divide the input.

Second, we need to check for more numbrs than just 22!

def is_it_prime(p): primality = True for n in srange(2,p): if n.divides(p): primality = False return primality

Did it work?

Now we are checking if nn divides pp for all whole numbers nn starting with 22 and going up to p1p-1. We should still do more testing. Go ahead, try this for your favoite numbers!

# code here

Into the wild

One of the most intriguing parts of math is that there are still unsolved problems! We often write programs to help us gather data to see what should be true and to help us solve these problems.

One open problem is:

Are there infinitely many twin primes?

Twin primes are pairs of prime numbers that are only two apart. For example, 33 and 55 are twin primes and 1111 and 1313 are twin primes.

That is, twin primes are defined are pairs of numbers pp and p+2p+2 where both pp and p+2p+2 are prime. Can you think of any more?

Let's think about it

It's pretty easy to think for a minute and find all twin primes less than 3030, or even 5050. But this isn't enough data to draw any conclusions. Let's write a function that prints out all pairs of twin primes up to a number MM.

To test for twin primes, we need to check that both pp and p+2p+2 are prime. We can use our is_it_prime function to do this!

def twin_primes(M): for n in srange(2,M+1): if is_it_prime(n) and is_it_prime(n+2): print(n,n+2)

Challenge

Can you use the above function to find out how many twin prime pairs there are up to 100? How about 1000? Do you think there are infinitely many twin primes?

#Try out you're code here

Another Challenge

Can you write a function that finds all pairs of primes that differ by 4? These are called cousin primes. Do you think there are infinitely many pairs of cousin primes?

How about pairs of primes that differ by 33? How many such pairs can you find? Can you think as to why the data for primes that differ by 33 looks so different than that for primes that differ by 22 or 44?