Introduction to Sage

631 days ago by jason3

Introduction to Sage

Based on a similar introduction to Mathematica by Ben Woodruff

The Cartesian Plane

Review

Introduction

Getting started:

Shift-enter or clicking "Evaluate"

* Tab-completion

* Question mark (and double-questionmark)

Defining variables

Sage understands many different types of mathematical structures.  In order to tell Sage that x, y, z and t are variables, we must do the following.

Let's start with some basics.  To define an expression in mathematica, use the = sign.  Below I define f to be Sin[x].
var('x,y,z,t') 
       
(x, y, z, t)
(x, y, z, t)

Entering functions

Let's start with some basics.  To define an expression in mathematica, use the `=` sign.  Below I define f to be the function \sin(x).

f(x)=sin(x) 
       

Now whenever I type f(x), it is replaced with \sin(x)

3*f(x)+4 
       
3*sin(x) + 4
3*sin(x) + 4

To evaluate the function at a point, simply use a number or other variable instead of x.

f(pi/3) 
       
1/2*sqrt(3)
1/2*sqrt(3)
f(z^2) 
       
sin(z^2)
sin(z^2)
f(y+2) 
       
sin(y + 2)
sin(y + 2)
 
       

You can "undefine" the function (or any variable) by "deleting" it.

del f 
       

Now the computer no longer knows what f is.

       
Traceback (click to the left of this block for traceback)
...
NameError: name 'f' is not defined
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/sage/sagenb/sage_notebook/worksheets/jason3/208/code/15.py", line 6, in <module>
    f
  File "", line 1, in <module>
    
NameError: name 'f' is not defined

To reset everything back to how it was when Sage started, go to the "Action" menu above and click "Restart Sage".

 
       

Sage understands many mathematical functions.  Note that exp(x) is the same as e^x, and log(x) is the same as \ln(x).

print cos(x), sin(pi) print e^(2*x), exp(2*x) print arctan(1), sqrt(9) print log(x), ln(x) 
       
cos(x) 0
e^(2*x) e^(2*x)
1/4*pi 3
log(x) log(x)
cos(x) 0
e^(2*x) e^(2*x)
1/4*pi 3
log(x) log(x)

In order to get a numeric value, use the n() function.

print sqrt(2) print n(sqrt(2)) print e^3 print (e^3).n() 
       
sqrt(2)
1.41421356237310
e^3
20.0855369231877
sqrt(2)
1.41421356237310
e^3
20.0855369231877

Derivatives, graphing, integration, and differentials

Let's take a function, plot it, differentiate and integrate it, and find critical values.

First, let's define the function.

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

Now let's plot it between x=-5 and x=4.

plot(f(x), (x, -5, 4)) 
       

Find the derivative two different ways.

diff(f(x), x) 
       
3*x^2 + 8*x
3*x^2 + 8*x
f(x).diff(x) 
       
3*x^2 + 8*x
3*x^2 + 8*x

Now let's plot the function and its derivative on top of each other.  To do this, we can just put the two functions in a list.

plot( [f(x), diff(f(x),x)], (x, -5, 4)) 
       

Another way to combine plots is to just add two separate plots together.

plot(f(x), (x, -5, 4))+plot(diff(f(x), x), (x, -5, 4), color='green') 
       
 
       

Integrate the function with an indefinite integral, without bounds.

integrate(f(x), x) 
       
1/4*x^4 + 4/3*x^3 + 2*x
1/4*x^4 + 4/3*x^3 + 2*x
f(x).integrate(x) 
       
1/4*x^4 + 4/3*x^3 + 2*x
1/4*x^4 + 4/3*x^3 + 2*x

Integrate the function with bounds.

integrate(f(x), x, -10, 10) 
       
8120/3
8120/3
f(x).integrate(x, -10, 10) 
       
8120/3
8120/3

Let's find the critical points by setting the derivative to zero and solving for x.

solve(f(x)==0,x) 
       
[x == -1/2*(I*sqrt(3) + 1)*(1/9*sqrt(3)*sqrt(155) - 91/27)^(1/3) +
1/9*(8*I*sqrt(3) - 8)/(1/9*sqrt(3)*sqrt(155) - 91/27)^(1/3) - 4/3, x ==
-1/2*(-I*sqrt(3) + 1)*(1/9*sqrt(3)*sqrt(155) - 91/27)^(1/3) +
1/9*(-8*I*sqrt(3) - 8)/(1/9*sqrt(3)*sqrt(155) - 91/27)^(1/3) - 4/3, x ==
(1/9*sqrt(3)*sqrt(155) - 91/27)^(1/3) + 16/9/(1/9*sqrt(3)*sqrt(155) -
91/27)^(1/3) - 4/3]
[x == -1/2*(I*sqrt(3) + 1)*(1/9*sqrt(3)*sqrt(155) - 91/27)^(1/3) + 1/9*(8*I*sqrt(3) - 8)/(1/9*sqrt(3)*sqrt(155) - 91/27)^(1/3) - 4/3, x == -1/2*(-I*sqrt(3) + 1)*(1/9*sqrt(3)*sqrt(155) - 91/27)^(1/3) + 1/9*(-8*I*sqrt(3) - 8)/(1/9*sqrt(3)*sqrt(155) - 91/27)^(1/3) - 4/3, x == (1/9*sqrt(3)*sqrt(155) - 91/27)^(1/3) + 16/9/(1/9*sqrt(3)*sqrt(155) - 91/27)^(1/3) - 4/3]

In order to substitute the values back in, we ask solve to give us the solutions in a specific form.

critical_points=solve(diff(f(x),x)==0,x, solution_dict=True) critical_points 
       
[{x: -8/3}, {x: 0}]
[{x: -8/3}, {x: 0}]

Then we substitute in each value from the critical points.

[f(x).subs(cp) for cp in critical_points] 
       
[310/27, 2]
[310/27, 2]

Some more plotting

Sage lets you tweak plots even more.  For example, you can specify the minimum and maximum y-values.

plot( [sin(x), cos(x), tan(x), e^x], (x, -5, 5), ymin=-10, ymax=10) 
       

Sometimes we want to plot functions without having to solve for y.  We can do this with `implicit_plot`.

implicit_plot(x^2+y^2==1, (x, -2, 2), (y, -2, 2), aspect_ratio=1) 
       

To work with differentials, you just need to remember that dy = f' dx. The code below shows how to use differentials to estimate the change in area of a circle caused by a change in the radius.

 

dr=1/10 A(r)=pi*r^2 dA=A.diff(r)*dr dA(r) 
       
1/5*pi*r
1/5*pi*r
n(dA(r=4)) 
       
-2.51327412287183
-2.51327412287183
dr=-1/10 c=4 dA=diff(A,r)(r=c)*dr plot(A(r), (r, c-dr, c+2*dr))+line([(c,A(c)), (c+dr, A(c)+dA), (c+dr, A(c)), (c, A(c))])+text("dA",(3.88,49))+text("dr",(3.97,50.5))+text("slope=dA/dr",(4,47)) 
       

Here are some more examples with limits, derivatives, integrals, and graphing.

f(x)=sin(x)/x 
       
limit(f(x), x=0) 
       
1
1
d=diff(f,x) print d(x) print d(x=2) 
       
cos(x)/x - sin(x)/x^2
-1/4*sin(2) + 1/2*cos(2)
cos(x)/x - sin(x)/x^2
-1/4*sin(2) + 1/2*cos(2)
g(x)=arctan(x) integrate(g(x), x) 
       
x*arctan(x) - 1/2*log(x^2 + 1)
x*arctan(x) - 1/2*log(x^2 + 1)
integrate(g(x), x, 0, 3) 
       
-1/2*log(10) + 3*arctan(3)
-1/2*log(10) + 3*arctan(3)
plot([f, g], (x, -10, 10)) 
       

The interact command lets you add sliders to most things you create, which les you see the effects of changing a variable by a little bit.  Try adjusting the slider after evaluating the next command.

@interact def f(a=(1,10)): show(plot(sin(a*x), (x, 0, 2*pi))) 
       

Click to the left again to hide and once more to show the dynamic interactive window

Now let's take any function, create a tangent line at a point, and graph the function and tangent line.

f=sin(2*x) c=1 f_prime=diff(f,x) tangent_line=f_prime(x=c)*(x-c)+f(x=c) show(plot([f,tangent_line], (x, -2*pi, 2*pi), ymin=-2, ymax=2)) 
       

We can make the plot interactive too.

@interact def myplot(f=sin(2*x),c=(-2*pi,2*pi)): f_prime=diff(f,x) tangent_line=f_prime(x=c)*(x-c)+f(x=c) show(plot([f,tangent_line], (x, -2*pi, 2*pi), ymin=-2, ymax=2)) 
       

Click to the left again to hide and once more to show the dynamic interactive window

Sage can automate almost all the calculations you did in first-semester calculus.

f(x)=x^3-3*x+2 first_derivative=diff(f, x) critical_points=solve(first_derivative==0, x) print "critical points: ", critical_points second_derivative=diff(f, x, x) hyper_critical_points=solve(second_derivative==0, x) print "hyper critical points: ", hyper_critical_points plot([f, first_derivative, second_derivative], (x, -3, 3)) 
       
critical points:  [
x == -1,
x == 1
]
hyper critical points:  [
x == 0
]
critical points:  [
x == -1,
x == 1
]
hyper critical points:  [
x == 0
]

The nth derivative can be found easily too.

f(x)=x^8 diff(f(x), x, 8) 
       
40320
40320

Solving equations and systems

The solver

Let’s use Sage to solve equations. One key thing to note is that you must type a double equal sign to enter equations into
the solver, as a single equal sign is used to assign expressions.
Let’s use the solver to find maximums and minimums of a function. Here we have to solve for when the derivative is zero. I’ll
graph the function, and factor the derivative as well just to help us see how Sage can be used.

f(x)=x^3-3*x+5 plot(f, (x, -3, 3)) 
       
print factor(diff(f(x), x)) 
       
3*(x - 1)*(x + 1)
3*(x - 1)*(x + 1)
critical_values=solve(diff(f(x),x)==0, x, solution_dict=True) critical_values 
       
[{x: -1}, {x: 1}]
[{x: -1}, {x: 1}]

Solve just told us that the critical values are +1 and -1.  Next, we evaluate the function at these points to see what the function values are.

[f(x).subs(solution) for solution in critical_values] 
       
[7, 3]
[7, 3]

This shows us that f(-1)=7 and f(1)=3.

The second derivative helps us determine the concavity at each point. If we evaluate the second derivative at each of the critical
points, we can use this information to tell us which point is a max, and which is a min.

[diff(f(x),x,x).subs(solution) for solution in critical_values] 
       
[-6, 6]
[-6, 6]

The concavity of x=-1 is negative (so we are at a max) and the concavity is +6 at 1 (so we are at a min).

Systems of Equations

You can solve a sytem of equation using the solve feature. Equations must contain two equal signs.

First, let's solve for x and y, then for x and z, then y and z.  Then let's specifically say that z=t and solve for x, y, and z.

var('x,y,z,t') solve((x+y==1, 3*x+4*y+z==0),[x,y]) 
       
[[x == z + 4, y == -z - 3]]
[[x == z + 4, y == -z - 3]]
solve((x+y==1, 3*x+4*y+z==0),[x,z]) 
       
[[x == -y + 1, z == -y - 3]]
[[x == -y + 1, z == -y - 3]]
solve((x+y==1, 3*x+4*y+z==0),[y,z]) 
       
[[y == -x + 1, z == x - 4]]
[[y == -x + 1, z == x - 4]]
solve((x+y==1, 3*x+4*y+z==0,z==t),[x,y,z]) 
       
[[x == t + 4, y == -t - 3, z == t]]
[[x == t + 4, y == -t - 3, z == t]]

Let's plot these in 3d.  Notice that the intersection of these two planes is a line, which is exactly what is represented by the last solve command.

implicit_plot3d(x+y==1,(x,-5,5),(y, -5,5),(z,-5,5))+implicit_plot3d(3*x+4*y+z==0, (x,-5,5), (y,-5,5), (z,-5,5),color='green') 
       

Matrices

Sage does matrix arithmetic as well.  One way to create a matrix is to specify a list of rows

A=matrix([[1,2],[3,4]]) B=matrix([[2,0,1],[1,-2,4]]) 
       
A*B 
       
[ 4 -4  9]
[10 -8 19]
[ 4 -4  9]
[10 -8 19]

The product BA should give an error since the number of columns in the first matrix is not equal to the number of rows in the second.

B*A 
       
Traceback (click to the left of this block for traceback)
...
TypeError: unsupported operand parent(s) for '*': 'Full MatrixSpace of 2
by 3 dense matrices over Integer Ring' and 'Full MatrixSpace of 2 by 2
dense matrices over Integer Ring'
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/sage/sagenb/sage_notebook/worksheets/jason3/208/code/98.py", line 6, in <module>
    B*A
  File "", line 1, in <module>
    
  File "element.pyx", line 1925, in sage.structure.element.Matrix.__mul__ (sage/structure/element.c:13028)
  File "coerce.pyx", line 740, in sage.structure.coerce.CoercionModel_cache_maps.bin_op (sage/structure/coerce.c:6627)
TypeError: unsupported operand parent(s) for '*': 'Full MatrixSpace of 2 by 3 dense matrices over Integer Ring' and 'Full MatrixSpace of 2 by 2 dense matrices over Integer Ring'

We can also calculate the determinant.

 

A.det() 
       
-2
-2
det(A) 
       
-2
-2

Curves

Parametric curves

Use parametric_plot to plot parametric curves

var('t') parametric_plot( (cos(t), sin(t)), (t, 0, 3*pi/2),aspect_ratio=1) 
       

Often people use the variable r to store both x and y (we’ll see why as we introduce vectors later in the semester). You should
avoid using x and y as variable names, as it can cause problems in computations. Hence, r = [x, y] allows me to store both x
and y with the same variable. Then r[0] will give the first element of the list, and r[1] will give the second. Anytime you store
information in a list, you can access the n th element by using r[n-1].  Notice that the indexing is zero-based (i.e., the first element of the list is r[0], and the last element of a list with 10 elements is r[9]).

r=[cos(t), sin(t)] r[0] 
       
cos(t)
cos(t)
r[1] 
       
sin(t)
sin(t)
parametric_plot(r, (t, 0, 3*pi/2),axes_labels=["x","y"], aspect_ratio=1) 
       
var('t') tmin=0 tmax=3*pi/2 @interact def myplot(r=input_box([cos(t), sin(t)]), tt=slider(tmin, tmax,default=(tmin+tmax)/2)): p=parametric_plot(r, (t, tmin, tmax),axes_labels=["x","y"],alpha=0.5) background=parametric_plot(r, (t, tmin, tt),thickness=3,aspect_ratio=1) show(p+background) 
       
tt 

Click to the left again to hide and once more to show the dynamic interactive window

Now let's find tangent lines to parametric curves.  First, we'll calculate \frac{dy}{dx}

var('x,y,t') r(t)=(cos(t), sin(t)) c=pi/4 
       
dydt=diff(r[1],t) dxdt=diff(r[0],t) dydx=dydt/dxdt dydx(t) 
       
-cos(t)/sin(t)
-cos(t)/sin(t)
tangent_line= y-r[1]==dydx*(x-r[0]) 
       
tangent_line(t=c) 
       
y - 1/2*sqrt(2) == -x + 1/2*sqrt(2)
y - 1/2*sqrt(2) == -x + 1/2*sqrt(2)
implicit_plot(tangent_line(t=c), (x,-2,2), (y,-2,2)) + parametric_plot(r, (t,0,2*pi),aspect_ratio=1) 
       

Here is an interactive version of the above calculation.  Notice that in order to not divide by zero when t=0, the equation of our tangent line is \frac{dx}{dt}(y-r[1])==\frac{dy}{dt}(x-r[0]), which we can get from the definition of \frac{dy}{dx}.

var('x,y,t') @interact def tanline(r=input_box([cos(t), sin(t)]), c=(0,2*pi)): dydt=diff(r[1],t) dxdt=diff(r[0],t) dydx=dydt tangent_line= dxdt*(y-r[1])==dydt*(x-r[0]) show(implicit_plot(tangent_line(t=c), (x,-2,2), (y,-2,2)) + parametric_plot(r, (t,0,2*pi),aspect_ratio=1)+point([r[0](t=c), r[1](t=c)],pointsize=30)) 
       

Click to the left again to hide and once more to show the dynamic interactive window

Polar Equations

The basic plot command is polar_plot.

var('theta') 
       
theta
theta
polar_plot(cos(3*theta), (theta, 0, pi), aspect_ratio=1) 
       

Graphing in polar coordinates is just a special case of using parametric curves.  In fact, polar_plot is programmed using parametric curves like the following.

r(theta)=cos(theta) x(theta)=r(theta)*cos(theta) y(theta)=r(theta)*sin(theta) parametric_plot( (x(theta), y(theta)), (theta, 0, pi),aspect_ratio=1) 
       

We can also make an interactive plot that lets us see how the plot is traced out as theta increases.

var('theta') @interact def f(r=1-2*cos(theta),t=slider(0,2*pi.n())): p=polar_plot(r, (theta, 0, t),thickness=3)+polar_plot(r, (theta,0,2*pi),alpha=0.5) p+=line([(0,0), (4*cos(t),4*sin(t))]) p+=line([(0,0), (-4*cos(t),-4*sin(t))],linestyle='--') p+=point((r(theta=t)*cos(t), r(theta=t)*sin(t)), rgbcolor='red',pointsize=15) show(p,xmin=-3,xmax=3,ymin=-3,ymax=3,aspect_ratio=1) 
       

Click to the left again to hide and once more to show the dynamic interactive window

Area inside of polar curves

Now lets calculate the area inside of a polar curve r=1-\cos(\theta).

var('theta') r=cos(2*theta) a=0 b=pi dA=1/2*r^2 
       
area=integrate(dA, theta,a,b) area 
       
1/4*pi
1/4*pi

And we can graph the area too.

polar_plot(r, (theta,a,b),fill=True,aspect_ratio=1) 
       

If you need to find the points of intersection of two polar curves, then remember to both solve algebraically (when the two curves intersect each other) and graphically.  Not every solution will appear from the algebraic solution.  In the example below, you only get two of the points of intersection from the algebraic solution. The solution (0,0) can only be obtained from the graph.

r1=1-cos(theta) r2=cos(theta) a=0 b=2*pi print "Intersection: ", solve(r1==r2, theta) polar_plot([r1,r2], (theta,a,b),aspect_ratio=1) 
       
Intersection:  [
theta == 1/3*pi
]
Intersection:  [
theta == 1/3*pi
]

If you want to find arc length, then remember the formula for ds. Just add up little bits of arc length.

r=cos(theta) a=0 b=pi polar_plot(r, (theta, a, b),aspect_ratio=1) 
       
ds=sqrt(r^2+diff(r,theta)^2) arclength=integrate(ds,theta,a,b) arclength 
       
pi
pi