Calculus in SAGE: A Quickstart Guide
Jamie Mulholland
Department of Mathematics
Simon Fraser University
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
- Differentiation
- Limits
- Solving Equations
- Newtons Method
- Integration
- Sequences
- Series
- Taylor Series
- Applied Project: Where to Sit at the Movies
Differential Calculus:
Integral Calculus:
Applications:
Let us start with the basics.
1) Syntax and Arthmetic
2+5
7 7 |
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:
2/3+5
17/3 17/3 |
(4/7)*(5/3)
20/21 20/21 |
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:
sqrt(2)+1
sqrt(2) + 1 sqrt(2) + 1 |
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".
n(sqrt(2)+1)
2.41421356237309 2.41421356237309 |
numerical_approx(sqrt(2)+1)
2.41421356237309 2.41421356237309 |
If you would like more digits then include the optional argument "digits" in the function call.
n(sqrt(2)+1,digits=20)
2.4142135623730950488 2.4142135623730950488 |
Sage also knows what \pi is:
pi
pi pi |
Notice it represents it exactly. If you want an approximation then you must use the n function.
n(pi)
3.14159265358979 3.14159265358979 |
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
a=3
|
|
a
3 3 |
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.
b=5
|
|
a+b
8 8 |
Look at that, Sage pulled a and b out of memory and then added them together.
We can overwrite variable names too.
We can overwrite variable names too.
a=1/2
a
1/2 1/2 |
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.
f(x)=x^2
f
x |--> x^2 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:
f(5)
25 25 |
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).
g(x)=a*x+1;
g
x |--> 1/2*x + 1 x |--> 1/2*x + 1 |
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
x |--> a*x + 1 x |--> a*x + 1 |
Notice we can evaluate f at different points as follows.
f(sqrt(2))
2 2 |
f(pi)
pi^2 pi^2 |
Again, Sage returns the exact values, no approximations here unless you ask for them.
Another example:
g(x)=sin(x)
|
|
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.
plot(f,(-1,1))
|
|
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
x |--> 4*x^3 x |--> 4*x^3 |
We can now evaluate the derivative at a point, say x=3.
f(3); df(3);
81 108 81 108 |
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).
diff(f,u)
x |--> 0 x |--> 0 |
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
x |--> 2*cos(x) 0 x |--> 2*cos(x) 0 |
higher derivatives:
To compute the n^{th} derivative of a function f(x) use the command
diff(f,x,n)
diff(f,x,2);
x |--> 12*x^2 x |--> 12*x^2 |
diff(f,x,3);
x |--> 24*x x |--> 24*x |
Sage handles complicated derivatives easily. Recall that we would do this example by using logarithmic differentiation.
diff(x^x,x)
(log(x) + 1)*x^x (log(x) + 1)*x^x |
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:
limit(x^2,x=2)
4 4 |
Indeterminate Form Examples:
limit((x^2-1)/(x-1),x=1)
2 2 |
limit(sin(x)/x,x=0)
1 1 |
limit((x^3+3)/(2*x^3+2*x-1),x=infinity)
1/2 1/2 |
limit(x^2/exp(x),x=infinity);
0 0 |
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')
+Infinity +Infinity |
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')
-Infinity -Infinity |
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)
[x == -2, x == -1] [x == -2, x == -1] |
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)
[[x == 5, y == 1]] [[x == 5, y == 1]] |
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)
0.78539816339744839 0.78539816339744839 |
See the Applied Project: Where to Sit at the Movies below for another example of where we use the find_root command.
[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
0.870073695796209 0.870073695796209 |
The third approximation is then:
x3=N(NewtonIt(x2)); x3
0.871221414262919 0.871221414262919 |
and the fourth is
x4=N(NewtonIt(x3)); x4
0.871221514495088 0.871221514495088 |
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
1/5*x^5 1/5*x^5 |
We can add limits of integration to compute the definite integral over an interval.
integral(f(x),x,1,3)
242/5 242/5 |
Here is another example, one that requires the method of integration by parts if we were to do by hand.
integral(x*cos(x),x)
x*sin(x) + cos(x) x*sin(x) + cos(x) |
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.)
n(h)
0.42467909778730545 0.42467909778730545 |
Alternatively, we could do this all in one step as follows.
n(integral(sin(x)/x, x, 1, pi/2))
0.42467909778730545 0.42467909778730545 |
[back to top]
Definite Integration: Interactive Applet
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.
%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]
Approximate (Numerical) Integration: Interactive Applet
Right-endpoint, Left-enpoint, Midpoint, and Trapezoid.
%hide
# By Graham Enos (based on the work of Nick Alexander (based on the work of
# Marshall Hampton))
var('x')
@interact
def riemann_rects(f = input_box(default = sin(x), type = SR,),
a = input_box(default = "0", type = str),
b = input_box(default = "pi/2", type = str),
n = slider(1, 20, 1, default=4, label="n"),
endpoint_rule = selector(['Left', 'Right', 'Midpoint', 'Trapezoid'],
nrows=1, label="Endpoint rule")):
# General Setup
a, b = sage_eval(a), sage_eval(b)
t = sage.calculus.calculus.var('t')
func = fast_callable(f(x=t), vars=[t])
dx = (b-a)/(n)
rects = Graphics()
# Exact answer (from Sage)
exact_answer = integral(func(t),t,a,b)
# Print only a few digits of precision
def cap(x):
return RealField(20)(x)
# Lists of x and y coordinates, depending on endpoint_rule
if endpoint_rule == 'Left':
xs = [q*dx + a for q in range(n)]
elif endpoint_rule == 'Right':
xs = [(q+1)*dx + a for q in range(n)]
elif endpoint_rule == 'Midpoint':
xs = [q*dx + a + dx/2 for q in range(n)]
else:
xs = [q*dx + a for q in range(n+1)]
ys = [func(x) for x in xs]
################################################
# The rest depends on rectangles vs trapezoids #
################################################
### Rectangles
if endpoint_rule in ['Left', 'Right', 'Midpoint']:
# Constructing rectangles
for q in range(n):
xm = q*dx + dx/2 + a
x = xs[q]
y = ys[q]
rects += line([[xm-dx/2,0],[xm-dx/2,y],[xm+dx/2,y],[xm+dx/2,0]],
rgbcolor = (1,0,0))
# Plotting: find max_y and min_y to determine height of plot
min_y = min(0, find_minimum_on_interval(func,a,b)[0])
max_y = max(0, find_maximum_on_interval(func,a,b)[0])
show(plot(func,a,b,thickness=2) + rects, xmin = a, xmax = b, ymin =
min_y, ymax = max_y)
# Pretty printing setup for sum and numerical values
sum_html = "%s • [%s]" % (cap(dx), ' + '.join([" f(%s)"
% cap(i) for i in xs]))
num_html = "%s • [%s]" % (cap(dx), ' + '.join([str(
cap(i)) for i in ys]))
# Estimated answer via Riemann Sums
estimated_answer = dx * sum([ys[q] for q in range(n)])
# Pretty printing, using \LaTeX
html(r'''
Integral = %s
Sum from 1 to %s of f(x_k) Δx
= %s
= %s
= %s
Absolute Error = %s
''' % (exact_answer, n, sum_html, num_html, estimated_answer.n(),
abs(exact_answer - estimated_answer).n()))
### Trapezoids
else:
# Constructing trapezoids
for q in range(n):
xL, xR = xs[q], xs[q+1]
yL, yR = ys[q], ys[q+1]
rects += line([[xL,0],[xL,yL],[xR,yR],[xR,0]], rgbcolor = (1,0,0))
# Plotting: find max_y and min_y to determine height of plot
min_y = min(0, find_minimum_on_interval(func,a,b)[0])
max_y = max(0, find_maximum_on_interval(func,a,b)[0])
show(plot(func,a,b,thickness=2) + rects, xmin = a, xmax = b, ymin =
min_y, ymax = max_y)
# Pretty printing setup for sum and numerical values; depends on endpoint_rule
sum_html = "%s/2 • [%s]" % (cap(dx), ' + '.join
(["[f(%s) + f(%s)] " % (cap(xs[i]), cap(xs[i+1])) for i in
range(n)]))
num_html = "%s/2 • [%s]" % (cap(dx), ' + '.join
([str(cap(i)) for i in ys]))
# Estimated answer via Riemann Sums
estimated_answer = dx/2 * sum(ys[i] + ys[i+1] for i in range(n))
# Pretty printing, using LaTeX
html(r'''
Integral = %s
Sum from 1 to %s of [f(x_i) + f(x_{i+1})]/2 • Δx
= %s
= %s
= %s
Absolute Error = %s
''' % (exact_answer, n, sum_html, num_html, estimated_answer.n(),
abs(exact_answer - estimated_answer).n()))
Click to the left again to hide and once more to show the dynamic interactive window |
Simpson's Rule:
def Simpson(f,a,b,n):
deltax=var("deltax"); simpsonarea=var("simpsonarea");
deltax=(b-a)/n;
simpsonarea=deltax/3*(f(x=a)+f(x=b)+2*sum([f(x=a+(2*j)*deltax) for j in range(1,n/2)])+4*sum([f(x=a+(2*j-1)*deltax) for j in range(1,n/2+1)]))
return N(simpsonarea);
|
|
Simpson(1/x,1,2,4), N(ln(2))
(0.693253968253968, 0.693147180559945) (0.693253968253968, 0.693147180559945) |
range(1,4)
[1, 2, 3] [1, 2, 3] |
n=var("n");
N(sum(1/n,n,1,100000))
12.0901461298634 12.0901461298634 |
var("n,ii");
for ii in [1000..1005]:
print N(sum(1/n,n,1,ii));
7.48547086055035 7.48646986154935 7.48746786554136 7.48846487451444 7.48946089045070 7.49045591532632 7.48547086055035 7.48646986154935 7.48746786554136 7.48846487451444 7.48946089045070 7.49045591532632 |
|
|
|
|
[back to top]
Applied Project: Where to Sit at the Movies
This example is taken from the Applied Project: Where to Sit at the Movies on page 446 of James Stewart's Calculus textbook (6th Ed.).
alpha,a,b=var('alpha,a,b');
alpha = pi/9;
a=(9+x*cos(alpha))^2+(31-x*sin(alpha))^2;
b=(9+x*cos(alpha))^2+(x*sin(alpha)-6)^2;
theta(x)=arccos((a+b-625)/(2*sqrt(a*b)));
theta
x |--> arccos(1/2*(2*(x*cos(1/9*pi) + 9)^2 + (x*sin(1/9*pi) - 31)^2 + (x*sin(1/9*pi) - 6)^2 - 625)/sqrt(((x*cos(1/9*pi) + 9)^2 + (x*sin(1/9*pi) - 6)^2)*((x*cos(1/9*pi) + 9)^2 + (x*sin(1/9*pi) - 31)^2))) x |--> arccos(1/2*(2*(x*cos(1/9*pi) + 9)^2 + (x*sin(1/9*pi) - 31)^2 + (x*sin(1/9*pi) - 6)^2 - 625)/sqrt(((x*cos(1/9*pi) + 9)^2 + (x*sin(1/9*pi) - 6)^2)*((x*cos(1/9*pi) + 9)^2 + (x*sin(1/9*pi) - 31)^2))) |
plot(theta,(0,60));
|
|
It looks like the maximum occurs somewhere slightly less than 10, say around 8 or 8.5.
We can now find this maximum by determining where the derivative is 0: \theta^{\prime}(x) = 0.
solve(theta.diff()==0,x)
[x == -9/cos(1/9*pi), x == -(9*sin(1/9*pi)^2 + 9*cos(1/9*pi)^2 + sqrt(27*sin(1/9*pi)^4 + 111*sin(1/9*pi)^3*cos(1/9*pi) + 89*sin(1/9*pi)^2*cos(1/9*pi)^2 + 111*sin(1/9*pi)*cos(1/9*pi)^3 + 62*cos(1/9*pi)^4)*sqrt(3))/(sin(1/9*pi)^2*cos(1/9*pi) + cos(1/9*pi)^3), x == -(9*sin(1/9*pi)^2 + 9*cos(1/9*pi)^2 - sqrt(27*sin(1/9*pi)^4 + 111*sin(1/9*pi)^3*cos(1/9*pi) + 89*sin(1/9*pi)^2*cos(1/9*pi)^2 + 111*sin(1/9*pi)*cos(1/9*pi)^3 + 62*cos(1/9*pi)^4)*sqrt(3))/(sin(1/9*pi)^2*cos(1/9*pi) + cos(1/9*pi)^3)] [x == -9/cos(1/9*pi), x == -(9*sin(1/9*pi)^2 + 9*cos(1/9*pi)^2 + sqrt(27*sin(1/9*pi)^4 + 111*sin(1/9*pi)^3*cos(1/9*pi) + 89*sin(1/9*pi)^2*cos(1/9*pi)^2 + 111*sin(1/9*pi)*cos(1/9*pi)^3 + 62*cos(1/9*pi)^4)*sqrt(3))/(sin(1/9*pi)^2*cos(1/9*pi) + cos(1/9*pi)^3), x == -(9*sin(1/9*pi)^2 + 9*cos(1/9*pi)^2 - sqrt(27*sin(1/9*pi)^4 + 111*sin(1/9*pi)^3*cos(1/9*pi) + 89*sin(1/9*pi)^2*cos(1/9*pi)^2 + 111*sin(1/9*pi)*cos(1/9*pi)^3 + 62*cos(1/9*pi)^4)*sqrt(3))/(sin(1/9*pi)^2*cos(1/9*pi) + cos(1/9*pi)^3)] |
Well, it looks like the solve command returned the exact solutions. This isn't too helpful here, since these look a little messy, so instead we compute the solutions numerically. This can be done using the find_root command. We can also specify the interval to look in: x in (0,60).
find_root(theta.diff()==0,0,60)
8.253062089733648 8.253062089733648 |
And there we have it!
|
|
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!
|
|
|
|
|
|
|