2018-01-24.ipynb

Kernel: SageMath 8.1

Math 157: Intro to Mathematical Software

UC San Diego, winter 2018

January 24, 2018: SageMath (part 2 of 3)

Administrivia:

• Homework 2 is due Friday, January 26 at 8pm. If you do not have it, contact course staff immediately.

• We are continuing to monitor the course waitlist. If you are on it and still want to join the course, please continue to keep up with lectures and homework, and watch your email for further instructions.

Finding Roots: symbolic

You can use the solve command to solve for zeroes of a function.

x^2+3 = 5  # This is an error; a Python assignment cannot have a complex expression on its left-hand side.

  File "<ipython-input-67-d6dc05fa28e5>", line 1
    x**Integer(2)+Integer(3) = Integer(5)  # This is an error; a Python assignment cannot have a complex expression on its left-hand side.
SyntaxError: can't assign to operator

x^2 + 3 == 5 ## This doesn't make sense in Python, but in Sage it is a valid object!

x^2 + 3 == 5

eqn = x^2 + 3 == 5
show(eqn)

f(x) = sin(x)+exp(x)

show(f)

eqn.add_to_both_sides(-3)

x^2 == 2

v = [(0,0), (1,1), (1,2), (2,0), (3,1)]
line(v, thickness=3, color='blue', zorder=-1) + points(v, pointsize=200, color='grey', zorder=1)

pi == pi

pi == pi

N(pi)

3.14159265358979

bool(pi == pi)

True

bool(eqn) # Should return False, because this is not an identity.

False

solve(x^4 + 2*x + 3 == 0, x) # There is an analogue of the quadratic formula for polynomials of degree 3 and 4...

[x == -1/2*sqrt(((2*I*sqrt(15) + 2)^(2/3) + 4)/(2*I*sqrt(15) + 2)^(1/3)) - 1/2*sqrt(-(2*I*sqrt(15) + 2)^(1/3) - 4/(2*I*sqrt(15) + 2)^(1/3) + 4/sqrt(((2*I*sqrt(15) + 2)^(2/3) + 4)/(2*I*sqrt(15) + 2)^(1/3))), x == -1/2*sqrt(((2*I*sqrt(15) + 2)^(2/3) + 4)/(2*I*sqrt(15) + 2)^(1/3)) + 1/2*sqrt(-(2*I*sqrt(15) + 2)^(1/3) - 4/(2*I*sqrt(15) + 2)^(1/3) + 4/sqrt(((2*I*sqrt(15) + 2)^(2/3) + 4)/(2*I*sqrt(15) + 2)^(1/3))), x == 1/2*sqrt(((2*I*sqrt(15) + 2)^(2/3) + 4)/(2*I*sqrt(15) + 2)^(1/3)) - 1/2*sqrt(-(2*I*sqrt(15) + 2)^(1/3) - 4/(2*I*sqrt(15) + 2)^(1/3) - 4/sqrt(((2*I*sqrt(15) + 2)^(2/3) + 4)/(2*I*sqrt(15) + 2)^(1/3))), x == 1/2*sqrt(((2*I*sqrt(15) + 2)^(2/3) + 4)/(2*I*sqrt(15) + 2)^(1/3)) + 1/2*sqrt(-(2*I*sqrt(15) + 2)^(1/3) - 4/(2*I*sqrt(15) + 2)^(1/3) - 4/sqrt(((2*I*sqrt(15) + 2)^(2/3) + 4)/(2*I*sqrt(15) + 2)^(1/3)))]

show(solve(x^4 + 2*x + 3 == 0, x)[0]) 

solve(x^5+x+1 == 0, x) # ... but unless you get lucky...

[x == -1/2*(1/18*sqrt(23)*sqrt(3) - 25/54)^(1/3)*(I*sqrt(3) + 1) - 1/18*(-I*sqrt(3) + 1)/(1/18*sqrt(23)*sqrt(3) - 25/54)^(1/3) + 1/3, x == -1/2*(1/18*sqrt(23)*sqrt(3) - 25/54)^(1/3)*(-I*sqrt(3) + 1) - 1/18*(I*sqrt(3) + 1)/(1/18*sqrt(23)*sqrt(3) - 25/54)^(1/3) + 1/3, x == (1/18*sqrt(23)*sqrt(3) - 25/54)^(1/3) + 1/9/(1/18*sqrt(23)*sqrt(3) - 25/54)^(1/3) + 1/3, x == -1/2*I*sqrt(3) - 1/2, x == 1/2*I*sqrt(3) - 1/2]

solve(x^5+3*x+1 == 0, x) # ... not in degree 5. I'll cover this in Math 100C next quarter!

[0 == x^5 + 3*x + 1]

solve(sqrt(x) == 2, x)

[x == 4]

solve(sin(x) == 0, x) # This doesn't give *all* of the solutions, but you can infer the rest.

[x == 0]

solve(e^(3*x) == 5, x)

[x == log(1/2*I*5^(1/3)*sqrt(3) - 1/2*5^(1/3)), x == log(-1/2*I*5^(1/3)*sqrt(3) - 1/2*5^(1/3)), x == 1/3*log(5)]

v = solve(e^(3*x) == 5, x, solution_dict=True); v

[{x: log(1/2*I*5^(1/3)*sqrt(3) - 1/2*5^(1/3))},
 {x: log(-1/2*I*5^(1/3)*sqrt(3) - 1/2*5^(1/3))},
 {x: 1/3*log(5)}]

var('x, y')
solve([x^2 == y^2, x^3 == y^3 + 5], [x,y])

[[x == 1.35720887245841, y == -1.35720887245841], [x == (-0.6786044041487247 + 1.175377306225595*I), y == (0.6786044041487285 - 1.175377306225602*I)], [x == (-0.6786044041487247 - 1.175377306225595*I), y == (0.6786044041487285 + 1.175377306225602*I)]]

solve([x^2 == y^2, x^3 == y^3 + 5], [x,y], solution_dict=True)

[{y: -1.35720887245841, x: 1.35720887245841},
 {y: 0.6786044041487285 - 1.175377306225602*I,
  x: -0.6786044041487247 + 1.175377306225595*I},
 {y: 0.6786044041487285 + 1.175377306225602*I,
  x: -0.6786044041487247 - 1.175377306225595*I}]

g = implicit_plot(x^3 == y^3 + 5, (x, -3, 3), (y,-3,3))
g += implicit_plot(x^2 == y^2, (x, -3, 3), (y,-3,3),color='red')
g

show(v[0][x])

e^(v[0][x]*3)

e^(3*log(1/2*I*5^(1/3)*sqrt(3) - 1/2*5^(1/3)))

(e^(v[0][x]*3)).simplify_full()

5

# or use roots:
v = (e^(3*x) - 5).roots()
v

[(log(1/2*I*5^(1/3)*sqrt(3) - 1/2*5^(1/3)), 1),
 (log(-1/2*I*5^(1/3)*sqrt(3) - 1/2*5^(1/3)), 1),
 (1/3*log(5), 1)]

show(v[0][0])

Here we numerically find a single zero of a function in a given interval.

f(x) = e^(3*x) - 5
plot(f, -2, 2, ymax=1)

f.find_root(0, 1)

0.5364793041447001

One can also turn a symbolic expression into a numerical approximation, to any desired accuracy (with some caution).

f = pi/10^30 + pi - e/10^30 - acos(0)
f

500000000000000000000000000001/1000000000000000000000000000000*pi - 1/1000000000000000000000000000000*e

numerical_approx(f)

1.57079632679490

alpha = log(1/2*I*5^(1/3)*sqrt(3) - 1/2*5^(1/3))
show(alpha)

numerical_approx(alpha)

0.536479304144700 + 2.09439510239320*I

numerical_approx(alpha, prec=200)

0.53647930414470012486691977774206254650853378475617257397088 + 2.0943951023931954923084289221863352561314462662500705473166*I

numerical_approx(alpha, digits=20)

0.53647930414470012487 + 2.0943951023931954923*I

# This is the number of digits used in computing the result, NOT the number of correct digits in output!
numerical_approx(alpha, digits=3)

0.536 + 2.09*I

N(9/2, prec=100)

4.5000000000000000000000000000

In order to obtain an answer which is guaranteed to have all reported digits be accurate (assuming no bugs in CoCalc, cosmic rays, or other irreproducible phenomena), one can use interval arithmetic. This is like keeping track of "error bars". It is somewhat less efficient, and so not commonly used in "real-life" applications; but it absolutely essentially if you want to use real-number arithmetic as the basis of a mathematical proof, as in the solution of the Kepler conjecture.

N is numerical_approx # Handy shorthand!

True

alpha.N()

alpha.n()

alpha.numerical_approx()

# Interval arithmetic --
# Every displayed digit except last in the output is definitely right:
ComplexIntervalField(20)(alpha)

0.53648? + 2.09440?*I

RIF = RealIntervalField()

u = RIF(sqrt(2))+RIF(sqrt(3))

u.lower(), u.upper()

(3.14626436994197, 3.14626436994198)

More about 2d plots

You can do much more than just plot(function...). E.g.,

• a point

• a bunch of points

• text

• "line" through a bunch of points

• polygon

• ellipse

• implicit plot

• contour plot

• vector field

point((1,2))

point2d((1,2), pointsize=150, color='red')

point([(random(), random()) for i in range(100)])

point([(x,x^2) for x in range(50)])

g =  point((1,2), pointsize=300, color='red')
g += text(r"Get $\int_0^{\pi} \sin(x) dx$ points!", (1,2), alpha=0.5, fontsize=30, fontweight='bold', color='black', rotation=20, zorder=1)
g.show(frame=True, axes=False)

line([(0,0), (1,1), (1,2), (2,0), (3,1)], thickness=3, color='blue')

v = [(0,0), (1,1), (1,2), (2,0), (3,1)]
line(v, thickness=3, color='blue', zorder=-1) + points(v, pointsize=200, color='grey', zorder=1)

var('x, y')
pl = implicit_plot(y^2 + cos(y*e^x) == x^3 - 2*x + 3, (x,-2,2), (y,-3,3))

pl.save('test.pdf') # To preserve your handiwork!

var('x,y')
plot3d(x^2+y^2, (x,-2,2),(y,-2,2))