# MATH 3600

## Introduction to Computational Mathematics

### Daniel D. Warner, Modified by Neil J. Calkin

Watch the video on Spirals, Fibonacci, and Being a Plant.

Is Ms. Vihart doing Mathematics?

What is the 100th Fibonacci number . . . the 1,000th . . . the 1,000,000 . . . the $10^{100}$?

Can we use the computer to calculate these numbers?

Can a computer count to $10^{100}$?

Are Ms. Vihart's rectangles related to the Golden Rectangle that was so prized by the ancient Greek architects?

If so, how?

What about drawing those spirals, can we get the computer to draw them?

Understanding how to effectively use a computer as a tool in doing mathematics is the subject of Computational Mathematics.

The first step is to learn how to program, and for this class we will use today's most powerful "platform" for programming in mathematics.

print 'Hello World!' print 'Cruel World!!'
 Hello World! Cruel World!! Hello World! Cruel World!!

# Introduction to Sage

## WWW.SAGEMATH.ORG

Sage is a platform that brings together many powerful open source software packages for doing computational mathematics, and we talk to these packages using the Python programming language.

Python is a high-level scripting language that is particularly well suited for supporting the exploratory nature of Computational Mathematics.

## We start by making sure that everyone can access Sage.

1. Access the Sage Server at  https://sage.math.clemson.edu:34567
2. Set up an account and click on the "Published Button."  All worksheets for this course will start with MATH 3600.
3. You can easily download a copy of Sage from WWW.SAGEMATH.ORG, and install it on your personal computer, but be sure that you have high speed connection and at least 3 GB of space available on your hard drive.
4. You can also create an account at sagemath.org and run Sage on their cloud servers.

## And now a short demo...

# This is a simple introduction to Sage. # # The Sage notebook uses blocks of text, which are essentially Python statements that get executed. # Lines that start with the pound sign are comments and do not get echoed when the block is evaluated. # # A block is evaluated by clicking on the evaluate link or by typing shift-enter. # # Pressing tab will give auto-complete options, or function descriptions after a parens ( # A '?' after a command gives detailed information, e.g. 'plot?' tells you about plot # Sage can do simple calculations 3+2*(6^2)
 75 75
# Sage can use variables googol = 10^100

Note the slightly different spelling here: explore why it is spelled that way.

2*googol googol 3*googol
 300000000000000000000000000000000000000000000000000000000000000000000000\ 00000000000000000000000000000 30000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
print(googol)
 100000000000000000000000000000000000000000000000000000000000000000000000\ 00000000000000000000000000000 10000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
show(googol)
 \newcommand{\Bold}[1]{\mathbf{#1}}10000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
print(2*googol+1)
 200000000000000000000000000000000000000000000000000000000000000000000000\ 00000000000000000000000000001 20000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000001
# and more elaborate calculations factorial(20)
 2432902008176640000 2432902008176640000

# Sage can factor integers factor(42627277772672262487989356835685788743568583563568538356835683568356857645627277)
 563 * 102299 * 305687981 * 1271181968890993 * 61409798072746003 * 31015876512649230749968872235979 563 * 102299 * 305687981 * 1271181968890993 * 61409798072746003 * 31015876512649230749968872235979
log(42627277772672262487989356835685788743568583563568538356835683568356857645627277.)/log(10.)
 79.6296875989841 79.6296875989841
# The same operations can be applied to algebraic expressions factor(x^6-1)
 (x^2 + x + 1)*(x^2 - x + 1)*(x + 1)*(x - 1) (x^2 + x + 1)*(x^2 - x + 1)*(x + 1)*(x - 1)
type(10)
  
a=10/3 print(a) type(a)
 10/3  10/3 
type((1-x)^5)
  
p = x^5-1 print(type(p)) print(p)
  x^5 - 1  x^5 - 1
print(p^3) print(p^4)
 (x^5 - 1)^3 (x^5 - 1)^4 (x^5 - 1)^3 (x^5 - 1)^4
plot(p)
# Sage can plot these algebraic expressions P1 = plot( p , -0.5, 0.5, rgbcolor=(1,0,0) ) P2 = plot( p^3, -0.5, 0.5, rgbcolor=(0,1,0) ) show(P1+P2)
factor(p^3)
 (x^4 + x^3 + x^2 + x + 1)^3*(x - 1)^3 (x^4 + x^3 + x^2 + x + 1)^3*(x - 1)^3
# Sage handles rational functions with ease rf = x / (x^2 - 1) print(rf) show(rf)
 x/(x^2 - 1) \newcommand{\Bold}[1]{\mathbf{#1}}\frac{x}{x^{2} - 1} x/(x^2 - 1)
P3 = plot( rf , -1.5, 1.5,ymin=-10,ymax=+10, rgbcolor=(1,0,1) ) show(P3)
P3 = plot( rf , -1.5, 1.5,ymin=-10,ymax=+10, color="cyan" ) show(P3)
plot(
 Traceback (click to the left of this block for traceback) ... SyntaxError: invalid syntax Traceback (most recent call last): File "", line 1, in File "_sage_input_51.py", line 10, in exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("cGxvdCg="),globals())+"\\n"); execfile(os.path.abspath("___code___.py")) File "", line 1, in File "/tmp/tmpubsbt1/___code___.py", line 3 ^ SyntaxError: invalid syntax
# What are x, p, and rf? print type(googol) print type(x) print type(p) print type(rf) # and what can we do with them?
    
# Partial fraction expansions print(rf.partial_fraction())
# Sage can display results in the more familiar mathematical format show(rf.partial_fraction())
rf2=y/(y^2-1)
 Traceback (click to the left of this block for traceback) ... NameError: name 'y' is not defined Traceback (most recent call last): File "", line 1, in File "_sage_input_3.py", line 10, in exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("cmYyPXkvKHleMi0xKQ=="),globals())+"\\n"); execfile(os.path.abspath("___code___.py")) File "", line 1, in File "/tmp/tmpHpNQCa/___code___.py", line 3, in exec compile(u'rf2=y/(y**_sage_const_2 -_sage_const_1 ) File "", line 1, in NameError: name 'y' is not defined
var('y') rf2=y/(y^2-1)
show(rf2)
 \newcommand{\Bold}[1]{\mathbf{#1}}\frac{y}{y^{2} - 1}

## Solving Algebraic Equations

# A quadratic example eq1 = x^2 - 2 == 3*x show(eq1) show(solve([eq1],x))
# Cubic polynomials are just as easy eq2 = x^3 + 2*x^2 - 5*x - 5 == 0 show(eq2) show(solve([eq2],x))
# Sage can even find the ugly answers eq3 = x^3 - 2 == x print(eq3) show(eq3) sol = solve([eq3],x) show(sol)
for s in sol: show(s)
# A linear systems of equations var('x,y,z') eq4 = x-1 == y+1+z eq5 = x+1 == 2*(y-1)+z show(eq4) show(eq5) eq4;eq5
sol = solve([eq4,eq5],x,y,z) print(sol)
 [ [x == r1 + 7, y == 5, z == r1] ] [ [x == r1 + 7, y == 5, z == r1] ]
# A nonlinear system of equations var('x y') eq6 = (x-1)^2 + (y-1)^2 == 4; show(eq6) eq7 = 4*y - 3*x == -1; show(eq7) sol3 = solve([eq6,eq7],x,y,solution_dict=True) show(sol3) for s in sol3: show(s[x])

## Now back to Fibonacci...

1 + 1 = 2

1 + 2 = 3

2 + 3 = 5

3 + 5 = 8

5 + 8 = 13

and so on . . .

# One approach # Use 3 variables and shift just like the example a = 1 b = 1 c = b + a print a,'+',b,'=',c
 1 + 1 = 2 1 + 1 = 2
# Do it again and again a = b b = c c = b + a print a,'+',b,'=',c a = b b = c c = b + a print a,'+',b,'=',c a = b b = c c = b + a print a,'+',b,'=',c a = b b = c c = b + a print a,'+',b,'=',c a = b b = c c = b + a print a,'+',b,'=',c a = b b = c c = b + a print a,'+',b,'=',c a = b b = c c = b + a print a,'+',b,'=',c a = b b = c c = b + a print a,'+',b,'=',c a = b b = c c = b + a print a,'+',b,'=',c a = b b = c c = b + a print a,'+',b,'=',c a = b b = c c = b + a print a,'+',b,'=',c a = b b = c c = b + a print a,'+',b,'=',c a = b b = c c = b + a print a,'+',b,'=',c a = b b = c c = b + a print a,'+',b,'=',c a = b b = c c = b + a print a,'+',b,'=',c a = b b = c c = b + a print a,'+',b,'=',c
 1 + 2 = 3 2 + 3 = 5 3 + 5 = 8 5 + 8 = 13 8 + 13 = 21 13 + 21 = 34 21 + 34 = 55 34 + 55 = 89 55 + 89 = 144 89 + 144 = 233 144 + 233 = 377 233 + 377 = 610 377 + 610 = 987 610 + 987 = 1597 987 + 1597 = 2584 1597 + 2584 = 4181 1 + 2 = 3 2 + 3 = 5 3 + 5 = 8 5 + 8 = 13 8 + 13 = 21 13 + 21 = 34 21 + 34 = 55 34 + 55 = 89 55 + 89 = 144 89 + 144 = 233 144 + 233 = 377 233 + 377 = 610 377 + 610 = 987 610 + 987 = 1597 987 + 1597 = 2584 1597 + 2584 = 4181
# Use a loop -- but make it stop a = 1 b = 1 c = b + a while c < 20000000: print a,'+',b,'=',c a = b b = c c = b + a
 1 + 1 = 2 1 + 2 = 3 2 + 3 = 5 3 + 5 = 8 5 + 8 = 13 8 + 13 = 21 13 + 21 = 34 21 + 34 = 55 34 + 55 = 89 55 + 89 = 144 89 + 144 = 233 144 + 233 = 377 233 + 377 = 610 377 + 610 = 987 610 + 987 = 1597 987 + 1597 = 2584 1597 + 2584 = 4181 2584 + 4181 = 6765 4181 + 6765 = 10946 6765 + 10946 = 17711 10946 + 17711 = 28657 17711 + 28657 = 46368 28657 + 46368 = 75025 46368 + 75025 = 121393 75025 + 121393 = 196418 121393 + 196418 = 317811 196418 + 317811 = 514229 317811 + 514229 = 832040 514229 + 832040 = 1346269 832040 + 1346269 = 2178309 1346269 + 2178309 = 3524578 2178309 + 3524578 = 5702887 3524578 + 5702887 = 9227465 5702887 + 9227465 = 14930352 1 + 1 = 2 1 + 2 = 3 2 + 3 = 5 3 + 5 = 8 5 + 8 = 13 8 + 13 = 21 13 + 21 = 34 21 + 34 = 55 34 + 55 = 89 55 + 89 = 144 89 + 144 = 233 144 + 233 = 377 233 + 377 = 610 377 + 610 = 987 610 + 987 = 1597 987 + 1597 = 2584 1597 + 2584 = 4181 2584 + 4181 = 6765 4181 + 6765 = 10946 6765 + 10946 = 17711 10946 + 17711 = 28657 17711 + 28657 = 46368 28657 + 46368 = 75025 46368 + 75025 = 121393 75025 + 121393 = 196418 121393 + 196418 = 317811 196418 + 317811 = 514229 317811 + 514229 = 832040 514229 + 832040 = 1346269 832040 + 1346269 = 2178309 1346269 + 2178309 = 3524578 2178309 + 3524578 = 5702887 3524578 + 5702887 = 9227465 5702887 + 9227465 = 14930352

# Suppose we want to generate the first 25 Fibonacci numbers # We can put in a counter a = 1 b = 1 c = b + a count = 3 while count <= 10^2: #print count, a,'+',b,'=',c a = b b = c c = b + a count = count + 1 print(c)
 573147844013817084101 573147844013817084101
# Same code but with some shorthand notation a = 1 b = 1 count = 3 while count <= 25: print a,'+',b,'=',(a+b) a, b = b, (a+b) count += 1
 1 + 1 = 2 1 + 2 = 3 2 + 3 = 5 3 + 5 = 8 5 + 8 = 13 8 + 13 = 21 13 + 21 = 34 21 + 34 = 55 34 + 55 = 89 55 + 89 = 144 89 + 144 = 233 144 + 233 = 377 233 + 377 = 610 377 + 610 = 987 610 + 987 = 1597 987 + 1597 = 2584 1597 + 2584 = 4181 2584 + 4181 = 6765 4181 + 6765 = 10946 6765 + 10946 = 17711 10946 + 17711 = 28657 17711 + 28657 = 46368 28657 + 46368 = 75025 1 + 1 = 2 1 + 2 = 3 2 + 3 = 5 3 + 5 = 8 5 + 8 = 13 8 + 13 = 21 13 + 21 = 34 21 + 34 = 55 34 + 55 = 89 55 + 89 = 144 89 + 144 = 233 144 + 233 = 377 233 + 377 = 610 377 + 610 = 987 610 + 987 = 1597 987 + 1597 = 2584 1597 + 2584 = 4181 2584 + 4181 = 6765 4181 + 6765 = 10946 6765 + 10946 = 17711 10946 + 17711 = 28657 17711 + 28657 = 46368 28657 + 46368 = 75025
range(10)
 [0, 1, 2, 3, 4, 5, 6, 7, 8, 9] [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
range?
range(2,25)
 [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24] [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24]
a=range(10) print(type(a[3])) b=srange(10) print(type(b[3]))
    
print(a) print(b) a==b
 [0, 1, 2, 3, 4, 5, 6, 7, 8, 9] [0, 1, 2, 3, 4, 5, 6, 7, 8, 9] True [0, 1, 2, 3, 4, 5, 6, 7, 8, 9] [0, 1, 2, 3, 4, 5, 6, 7, 8, 9] True
a = 1 b = 1 for count in range(2,25): c = b + a print a,'+',b,'=',c a = b b = c
 1 + 1 = 2 1 + 2 = 3 2 + 3 = 5 3 + 5 = 8 5 + 8 = 13 8 + 13 = 21 13 + 21 = 34 21 + 34 = 55 34 + 55 = 89 55 + 89 = 144 89 + 144 = 233 144 + 233 = 377 233 + 377 = 610 377 + 610 = 987 610 + 987 = 1597 987 + 1597 = 2584 1597 + 2584 = 4181 2584 + 4181 = 6765 4181 + 6765 = 10946 6765 + 10946 = 17711 10946 + 17711 = 28657 17711 + 28657 = 46368 28657 + 46368 = 75025 1 + 1 = 2 1 + 2 = 3 2 + 3 = 5 3 + 5 = 8 5 + 8 = 13 8 + 13 = 21 13 + 21 = 34 21 + 34 = 55 34 + 55 = 89 55 + 89 = 144 89 + 144 = 233 144 + 233 = 377 233 + 377 = 610 377 + 610 = 987 610 + 987 = 1597 987 + 1597 = 2584 1597 + 2584 = 4181 2584 + 4181 = 6765 4181 + 6765 = 10946 6765 + 10946 = 17711 10946 + 17711 = 28657 17711 + 28657 = 46368 28657 + 46368 = 75025

Print the 1,000th Fibonacci number, but not all the others.

a = 1 b = 1 for count in range(2,1000): c = b + a a,b = b,c print a,'+',b,'=',c
 268638100244853593861467272021429239676166093189869523401231759976179817\ 002478816893383696544833565641918278561614433563129766736422103503246348\ 50410377680367334151172899169723197082763985615764450078474174626 + 434665576869374564356885276750406258025646605173717804024817290895365554\ 179490518904038798400792551692959225930803226347752096896232398733224711\ 61642996440906533187938298969649928516003704476137795166849228875 = 434665576869374564356885276750406258025646605173717804024817290895365554\ 179490518904038798400792551692959225930803226347752096896232398733224711\ 61642996440906533187938298969649928516003704476137795166849228875 26863810024485359386146727202142923967616609318986952340123175997617981700247881689338369654483356564191827856161443356312976673642210350324634850410377680367334151172899169723197082763985615764450078474174626 + 43466557686937456435688527675040625802564660517371780402481729089536555417949051890403879840079255169295922593080322634775209689623239873322471161642996440906533187938298969649928516003704476137795166849228875 = 43466557686937456435688527675040625802564660517371780402481729089536555417949051890403879840079255169295922593080322634775209689623239873322471161642996440906533187938298969649928516003704476137795166849228875

How long do you think it would take to compute the 1,000,000th Fibonacci number?

Give it a shot!

How long do you think it would take to compute the $10^{100}$th Fibonacci number?

10^90/(3.15*10^7)
 3.17460317460317e82 3.17460317460317e82

A fast computer --- 10 GHz has $10^{10}$ cycles per second.  Computing the  $10^{100}$th Fibonacci number

takes $10^{100}$ additions.  Suppose each addition mod 20 takes one cycle.  That means we take $10^{90}$ seconds.

That's about $3\times 10^{80}$ centuries.

## Arrays and Lists

An array is an ordered homogeneous collection.  A list is an ordered inhomogeneous collection.  In low level languages like C when an array is set up it has a fixed size.  In higher level languages like Matlab, arrays are dynamic.  That is, they can grow as needed.  In Python lists are both inhomogeneous and dynamic.  They generalize dynamic arrays.

We will start by considering some array operations that are common in most languages, and then explore the additional features provided in Python.

# The range command gives us an list of integers a = range(15)
type(a)
  
 [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14] [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14]
len(a)
 15 15
a[8] = 2*a[7] a[9] = 3*a[7] a
 [0, 1, 2, 3, 4, 5, 6, 7, 14, 21, 10, 11, 12, 13, 14] [0, 1, 2, 3, 4, 5, 6, 7, 14, 21, 10, 11, 12, 13, 14]
# Generate an list of the first 25 Fibonacci numbers a = [0 for k in range(25)] print a a[0] = 1 a[1] = 1 for k in range(2,25): a[k] = a[k-1] + a[k-2] print a
 [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0] [1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765, 10946, 17711, 28657, 46368, 75025] [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0] [1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765, 10946, 17711, 28657, 46368, 75025]
# What can we learn about a? type(a)
  
a = [1,1] for k in range(2,25): a.append(a[k-1]+a[k-2]) a
 [1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765, 10946, 17711, 28657, 46368, 75025] [1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765, 10946, 17711, 28657, 46368, 75025]
# Entries in a list can be any object that Sage can work with. a = [0, -15, sqrt(2), 3/5, 'John Doe', [1,2,3]] a
 [0, -15, sqrt(2), 3/5, 'John Doe', [1, 2, 3]] [0, -15, sqrt(2), 3/5, 'John Doe', [1, 2, 3]]
for k in range(len(a)): print a[k], type(a[k])
 0 -15 sqrt(2) 3/5 John Doe [1, 2, 3]  0 -15 sqrt(2) 3/5 John Doe [1, 2, 3] 
for k in range(len(a)): show(a[k], type(a[k]))
a+b
 [1, 2, 3, 5, 6, 'Philadelphia'] [1, 2, 3, 5, 6, 'Philadelphia']
b+a
 [5, 6, 'Philadelphia', 1, 2, 3] [5, 6, 'Philadelphia', 1, 2, 3]
a = [1,1] for k in range(2,25): a=a+[(a[k-1]+a[k-2])] print(a)
 [1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765, 10946, 17711, 28657, 46368, 75025] [1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765, 10946, 17711, 28657, 46368, 75025]
a = [1,1] for k in range(2,25): a += [(a[k-1]+a[k-2])] print(a)
def fib(n): if n==0: return(1) elif n==1: return(1) else: return(fib(n-1)+fib(n-2))

fib(25)
 121393 121393
fib(50)
fibonacci(2500)
 131709051675194962952276308712531641206660696499250714188774693672753087\ 040503842576450313012318640774657086218587192595276683635211911952815631\ 558263246079038383460565488061265771846563256883924597824847305817942207\ 073555312471638545088664055239227385677067223979716426435692766130834967\ 194167364320573334359270171671578825517067957550027918605331636558325918\ 692735935102338729837168622286082741537144355375995365951412088276380814\ 259336640247225134836000891558521529150498437169752387119955393571405695\ 9634778700594751875 1317090516751949629522763087125316412066606964992507141887746936727530870405038425764503130123186407746570862185871925952766836352119119528156315582632460790383834605654880612657718465632568839245978248473058179422070735553124716385450886640552392273856770672239797164264356927661308349671941673643205733343592701716715788255170679575500279186053316365583259186927359351023387298371686222860827415371443553759953659514120882763808142593366402472251348360008915585215291504984371697523871199553935714056959634778700594751875