Calculus in SAGE: A Quickstart Guide
Jamie Mulholland
Department of Mathematics
Simon Fraser University
Sage is a free open-source mathematics software system licensed under the GPL. It combines the power of many existing open-source packages into a common Python-based interface. Their mission is to create a viable free open source alternative to Magma, Maple, Mathematica and Matlab. Sage is extremely powerful and is extensively used at the research level in mathematics.
To download a free copy of Sage point your browser to
http://www.sagemath.org/. Choose to download from the SFU mirror site. Installation instructions are included in the download. I will point out that Sage can run from the command line, or through a notebook (what you are reading right now). The notebook runs in your web browser.
I'll assume now that you have opened the Sage program and are using the notebook interface (this includes having just opened this page in your web browser). You can view the interactive Sage tutorial by clicking Help, then click Tutorial. However, for those who wish for an even quicker tutorial the document you are reading now will give you a quick look at how Sage can tackle exercises in Calculus.
Contents:
Basics:
- Syntax and Arthmetic
- Getting Help
- Storing Items in Memory
- Functions and Graphs
Differential Calculus:
- Differentiation
- Limits
- Solving Equations
- Newtons Method
Integral Calculus:
- Integration
- Sequences
- Series
- Taylor Series
Let us start with the basics.
1) Syntax and Arthmetic
The text in black (i.e. the text next to the red vertical bar in the box, a.k.a execution block, is the "input". This is a command that you are giving
to Sage. Pressing "shift-enter" then gets Sage to execute the command and the output is returned in blue.
Here are some more examples:
Notice Sage returns the exact value (i.e a fraction) not just a decimal approximation, which is what
most calculators would return. So Sage works with exact expressions, not just approximations.
Note: " * " represents multiplication
" / " represents division
" ^ " represents exponentiation
Sage works with exact expression as seen here:
If you would like a decimal approximation of \sqrt{2}+1 then use the "n" function. Or you could use its longer name "numerical_approx".
numerical_approx(sqrt(2)+1)
If you would like more digits then include the optional argument "digits" in the function call.
Sage also knows what \pi is:
Notice it represents it exactly. If you want an approximation then you must use the n function.
For you to try:
1) Evaluate
3^2+\frac{\pi}{2} to 10 digits. To 20 digits.
[
back to top]
2) Getting Help:
If you ever need help on a topic just type
[topic name]?
For example
sin?
would give you help on the sine function.
[
back to top]
3) Storing Items in Memory
In the first statement we stored the value 3 in the variable name a. In the next execution block we
called the value of a back out of memory.
Let's do another.
Look at that, Sage pulled a and b out of memory and then added them together.
We can overwrite variable names too.
For you to try:
2) Store the value
\frac{3}{4} under the name
a and the value
-\frac{9}{2} under the name
b. Then get Sage to compute
a+b,
a-b,
ab, and
\frac{a}{b}.
[
back to top]
4) Functions and Graphs:
Let us declare
f to be the function given by
f(x)=x^2.
Perfect, we can just use good old function notation.
(Note the second f is there just so it will display the value of f. I've just done this so you can see it read in and stored the value of f properly. We don't need to do this from now on.)
Now what if we want to determine f(5). Simple, this is done just as we would expect:
We can also input expression into the function.
f(x+h);
expand(f(x+h)); #(this is a comment) here we ask Sage to expand the expression for f(x+h)
|
|
(h + x)^2
h^2 + 2*h*x + x^2
(h + x)^2
h^2 + 2*h*x + x^2
|
Let us create a new function with the parameter a, say g(x)=a*x+1. This is just a family of linear functions with the same y-intercept but different slopes (the slope is a which is a parameter).
This may not be what we expected. What happened? Oh right, I already stored a value for the variable "a" above. It seems to have just used it here. Ok, so how can we clear the value from a variable and use the variable in a function? We do this as follows:
a=var('a') # here we declared a to be a variable and hence cleared its value.
g(x)=a*x+1
g
Notice we can evaluate f at different points as follows.
Again, Sage returns the exact values, no approximations here unless you ask for them.
Another example:
g(0); g(pi/3); g(pi/4); g(pi/6);
|
|
0
1/2*sqrt(3)
1/2*sqrt(2)
1/2
0
1/2*sqrt(3)
1/2*sqrt(2)
1/2
|
Wonderful, it even knows the exact values of trigonometric functions at special angles.
Hmmm... what about the sum of two functions? Or the product? Or the composition?
(f+g)(1); (f*g)(2); f(g(x)); g(f(x));
|
|
sin(1) + 1
4*sin(2)
sin(x)^2
sin(x^2)
sin(1) + 1
4*sin(2)
sin(x)^2
sin(x^2)
|
Yes! Sage can do these too. Notice the syntax for composition is exactly what we would expect:
(f\circ g)(x) = f(g(x)).
Some functions Sage knows:
| function |
Sage notation |
| \sqrt{x} |
sqrt(x) |
| e^x |
exp(x) |
| \ln{x} |
ln(x) or log(x) |
| \log_n(x) (base n logarithm) |
log(x,n) |
| sin(x) |
sin(x) |
| cos(x) |
cos(x) |
| tangent |
tan(x) |
It also knows arcsin, arccsos, acrtan, sinh, cosh, and many many more.
Remember, if you ever need help on a topic just type:
[topic name]?
For you to try:
3) Define
h(x)=\cos{x}+e^{x^2}+x+3 and evaluate
h(0) and
h(-\frac{\pi}{3}).
How about plotting the functions
f and
g that we defined above?
This is done using the "plot" command.
Notice below that I have specified that the function only be plotted for x values from -1 to 1. You can
input whatever domain you like.
Remember, type "plot?" to get more information on plotting. There are a lot of options you can use to
customize your plots.
Another couple of examples are:
plot(g,(-pi,2*pi),rgbcolor='red')
plot(x^3+x+1,(-2,2),rgbcolor='#ff00ff',)
We can plot two functions on the same coordinate system.
plot([x^2,0.5*x+1],(-2,2))
However, if you want to color the graphs differently, or use different domains for each function, then you can do this by creating the plots individually, then creating an new plot which is the join of the two plots (using the + operator), and then use the 'show' method to display the plots.
p1=plot(x^2,(-2,2));
p2=plot(0.5*x+1,(-2,4),rgbcolor='green');
p=p1+p2;
p.show()
For you to try:
4) Plot the function
h(x)=\cos{x}+e^{x^2}+x+3 on the interval
[-1,1]. Then plot the function
f(x)=x^2 on the same set of axis but color it purple.
[
back to top]
5) Differentiation:
After all this is Calculus, so we should show how to use Sage to do differentiation.
f(x)=x^4;
df=diff(f,x); # here we take the derivative of f with respect to x and call it df
df
We can now evaluate the derivative at a point, say x=3.
Notice that the command to differentiate explicitly mentions the variable x. This is the variable we used to define the function f. If we differentiate with respect to another variable, say u, then we would expect to get 0 as a result (since f is constant with respect to u).
Another example:
h(x)=2*sin(x);
dh=diff(h,x); # take derivative and store it under the name dh
dh; dh(pi/2); # call to print dh and the value of dh at pi/2
higher derivatives:
To compute the n^{th} derivative of a function f(x) use the command
diff(f,x,n)
Sage handles complicated derivatives easily. Recall that we would do this example by using logarithmic differentiation.
For you to try:
5) Pick some derivative questions from the textbook and have Sage compute them.
[
back to top]
6) Limits:
Sage can certainly do limits, and quite well I might add.
The command to compute
\lim_{x \to a} f(x) is given by:
limit(f,x=a)
Here are some examples:
Indeterminate Form Examples:
limit((x^3+3)/(2*x^3+2*x-1),x=infinity)
limit(x^2/exp(x),x=infinity);
To compute the right-hand limit. lim_{x \to a^+} f(x) use the argument dir='plus' as follows:
limit((x+1)/(x-7),x = 7, dir='plus')
To compute the left-hand limit. lim_{x \to a^-} f(x) use the argument dir='minus' as follows:
limit(cos(x)/(sin(x)-1),x=pi/2,dir='minus')
For you to try:
6) Pick some limit questions from the textbook and have Sage compute them. In particular choose some questions from the section on L'Hospital's Rule.
[
back to top]
7) Solving Equations:
Solving Equations Exactly
x=var('x')
solve(x^2+3*x+2==0,x)
Solve for one variable in terms of others.
x, b, c, = var('x, b, c')
solve(x^2+b*x+c==0,x)
|
|
[x == -1/2*b - 1/2*sqrt(b^2 - 4*c), x == -1/2*b + 1/2*sqrt(b^2 - 4*c)]
[x == -1/2*b - 1/2*sqrt(b^2 - 4*c), x == -1/2*b + 1/2*sqrt(b^2 - 4*c)]
|
Solve for several variables.
x, y = var('x, y')
solve([x+y==6,x-y==4],x,y)
Solving Equations Numerically
When solve can't find the exact solution, find_root can approximate the solution. But you need to specify a interval in which to look. In the next example we look in 0< \theta < \frac{\pi}{2}.
theta=var('theta')
find_root(cos(theta)==sin(theta),0,pi/2)
[
back to top]
8) Newton's Method:
Let us see how we can use out newfound knowledge to get Sage to solve a Newton's method problem.
Problem: Use Newton's Method with an initial guess
x_1=1 to approximate the root of
5x+\cos{x}=5 which lies in the interval
[0,\frac{\pi}{2}] to four decimal places.
The function we with to find the root of is
f(x)=\cos{x}+5x-5. So lets put this into Sage and also input Newton's Iteration formula:
x, xn = var('x, xn');
f(x)=cos(x)+5*x-5;
df=diff(f,x);
NewtonIt(x)=x-(f/df)(x); #Newton's Iterative Formula
Now that Newton's Iterative formula is ready to go, we can start to compute the sequence of approximations to the root.
x1=1;
x2=N(NewtonIt(x1)); x2
The third approximation is then:
and the fourth is
and so on.
We could also use a loop to save us from having to do a lot of typing.
Sage is implemented using Python so we use Python commands to construct loops. So two quick points about loops and list in Python you should know are: In the for loop below i runs over integer values in the range [0,5), i.e. over the number 0, 1, 2, 3, 4. The elements of a list in Python are indexed by 0, 1, 2, etc. (i.e. the first element has index 0, not 1).
xnlist=[1]
for i in range(0,5):
xnlist.append(N(NewtonIt(xnlist[i])))
xnlist
|
|
[1, 0.870073695796209, 0.871221414262919, 0.871221514495088,
0.871221514495089, 0.871221514495089]
[1, 0.870073695796209, 0.871221414262919, 0.871221514495088, 0.871221514495089, 0.871221514495089]
|
Since the fifth and sixth approximations agree to 14 decimal places then we have found the solution to the equation to 14 decimal places:
0.87122151449508.
[
back to top]
9) Integration
The command for taking the integral of a function is simply "integral".
f(x)=x^4;
integral(f(x),x); # here we take the indefinite integral (antidervivative) of f with respect to x
We can add limits of integration to compute the definite integral over an interval.
Here is another example, one that requires the method of integration by parts if we were to do by hand.
As we have seen, some functions do not have antiderivatives expressible in terms of elementary functions. For example the antiderivative of e^{x^2} is not able to be expressed using elementary functions: polynomials, exponentials, logarithms, trigonometric functions.
To evaluate definite integrals of such functions we use numerical techniques: Riemann sums, Trapezoid Rule, Simpson's Rule, etc.
Here is an example where we compute the integral of a function numerically. Notice the function "n" returns the numerical value of the integral.
h = integral(sin(x)/x, x, 1, pi/2); h
|
|
integrate(sin(x)/x, x, 1, 1/2*pi)
integrate(sin(x)/x, x, 1, 1/2*pi)
|
(Notice in returns exactly what we entered. This indicates the antiderivative is not expressible in terms of elementary functions.)
Alternatively, we could do this all in one step as follows.
n(integral(sin(x)/x, x, 1, pi/2))
Here is an interactive applet in which you can input one, or two, functions and it will compute the definite integral and display a nice plot.
Interactive: Definite Integral:
%hide
# Definite Integral
# Lauri Ruotsalainen, 2010
@interact
def definite_integral(
f = input_box(default = 3*x),
g = input_box(default = x^2),
interval = input_box(default=(0,3), label="Interval"),
x_range = input_box(default=(0,3), label = "Plotting range (x)"),
selection = selector(["f", "g", "f and g", "f - g"], label="Select")
):
f(x) = f; g(x) = g
f_plot = Graphics(); g_plot = Graphics(); h_plot = Graphics();
text = ""
# Plot function f.
if selection != "g":
f_plot = plot(f(x), x, x_range, color="blue", thickness=1.5)
# Color and calculate the area between f and the horizontal axis.
if selection == "f" or selection == "f and g":
f_plot += plot(f(x), x, interval, color="blue", fill=True, fillcolor="blue", fillalpha=0.15)
text += "$\int_%s^%s(\color{Blue}{f(x)})\,dx=\int_%s^%s(%s)\,dx=%s$" % (
interval[0], interval[1],
interval[0], interval[1],
latex(f(x)), #produces latex code for expression f
f(x).nintegrate(x, interval[0], interval[1])[0]
)
# Plot function g. Also color and calculate the area between g and the horizontal axis.
if selection == "g" or selection == "f and g":
g_plot = plot(g(x), x, x_range, color="green", thickness=1.5)
g_plot += plot(g(x), x, interval, color="green", fill=True, fillcolor="yellow", fillalpha=0.5)
text += "\n$\int_%s^%s(\color{Green}{g(x)})\,dx=\int_%s^%s(%s)\,dx=%s$" % (
interval[0], interval[1],
interval[0], interval[1],
latex(g(x)), #produces latex code for expression g
g(x).nintegrate(x, interval[0], interval[1])[0]
)
# Plot function f-g. Also color and calculate the area between f-g and the horizontal axis.
if selection == "f - g":
g_plot = plot(g(x), x, x_range, color="green", thickness=1.5)
g_plot += plot(g(x), x, interval, color="green", fill=f(x), fillcolor="red", fillalpha=0.15)
h_plot = plot(f(x)-g(x), x, interval, color="red", thickness=1.5, fill=True, fillcolor="red", fillalpha=0.15)
text = "$\int_%s^%s(\color{Red}{f(x)-g(x)})\,dx=\int_%s^%s(%s)\,dx=%s$" % (
interval[0], interval[1],
interval[0], interval[1],
latex(f(x)-g(x)), #produces latex code for expression f-g
(f(x)-g(x)).nintegrate(x, interval[0], interval[1])[0]
)
show(f_plot + g_plot + h_plot, gridlines=True)
html(text)
|
|
Click to the left again to hide and once more to show the dynamic interactive window
|
[
back to top]
11) Series
...coming soon...
This has been a very brief introduction to Sage to get you up and running (ok, maybe a slow jog). You can find out a lot more about Sage and Python by working through the interactive Sage tutorial by clicking Help, then click Tutorial. Have fun!