Basic Algebra and Calculus¶
Sage can perform various computations related to basic algebra and calculus: for example, finding solutions to equations, differentiation, integration, and Laplace transforms. See the Sage Constructions documentation for more examples.
In all these examples, it is important to note that the variables in the
functions are defined to be var(...)
. As an example:
sage: u = var('u')
sage: diff(sin(u), u)
cos(u)
>>> from sage.all import *
>>> u = var('u')
>>> diff(sin(u), u)
cos(u)
u = var('u') diff(sin(u), u)
If you get a NameError
, check to see if you misspelled something,
or forgot to define a variable with var(...)
.
Solving Equations¶
Solving Equations Exactly¶
The solve
function solves equations. To use it, first specify
some variables; then the arguments to solve
are an equation (or a
system of equations), together with the variables for which to
solve:
sage: x = var('x')
sage: solve(x^2 + 3*x + 2, x)
[x == -2, x == -1]
>>> from sage.all import *
>>> x = var('x')
>>> solve(x**Integer(2) + Integer(3)*x + Integer(2), x)
[x == -2, x == -1]
x = var('x') solve(x^2 + 3*x + 2, x)
You can solve equations for one variable in terms of others:
sage: x, b, c = var('x b c')
sage: 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)]
>>> from sage.all import *
>>> x, b, c = var('x b c')
>>> solve([x**Integer(2) + b*x + c == Integer(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, b, c = var('x b c') solve([x^2 + b*x + c == 0],x)
You can also solve for several variables:
sage: x, y = var('x, y')
sage: solve([x+y==6, x-y==4], x, y)
[[x == 5, y == 1]]
>>> from sage.all import *
>>> x, y = var('x, y')
>>> solve([x+y==Integer(6), x-y==Integer(4)], x, y)
[[x == 5, y == 1]]
x, y = var('x, y') solve([x+y==6, x-y==4], x, y)
The following example of using Sage to solve a system of non-linear equations was provided by Jason Grout: first, we solve the system symbolically:
sage: var('x y p q')
(x, y, p, q)
sage: eq1 = p+q==9
sage: eq2 = q*y+p*x==-6
sage: eq3 = q*y^2+p*x^2==24
sage: solve([eq1,eq2,eq3,p==1],p,q,x,y)
[[p == 1, q == 8, x == -4/3*sqrt(10) - 2/3, y == 1/6*sqrt(10) - 2/3], [p == 1, q == 8, x == 4/3*sqrt(10) - 2/3, y == -1/6*sqrt(10) - 2/3]]
>>> from sage.all import *
>>> var('x y p q')
(x, y, p, q)
>>> eq1 = p+q==Integer(9)
>>> eq2 = q*y+p*x==-Integer(6)
>>> eq3 = q*y**Integer(2)+p*x**Integer(2)==Integer(24)
>>> solve([eq1,eq2,eq3,p==Integer(1)],p,q,x,y)
[[p == 1, q == 8, x == -4/3*sqrt(10) - 2/3, y == 1/6*sqrt(10) - 2/3], [p == 1, q == 8, x == 4/3*sqrt(10) - 2/3, y == -1/6*sqrt(10) - 2/3]]
var('x y p q') eq1 = p+q==9 eq2 = q*y+p*x==-6 eq3 = q*y^2+p*x^2==24 solve([eq1,eq2,eq3,p==1],p,q,x,y)
For numerical approximations of the solutions, you can instead use:
sage: solns = solve([eq1,eq2,eq3,p==1],p,q,x,y, solution_dict=True)
sage: [[s[p].n(30), s[q].n(30), s[x].n(30), s[y].n(30)] for s in solns]
[[1.0000000, 8.0000000, -4.8830369, -0.13962039],
[1.0000000, 8.0000000, 3.5497035, -1.1937129]]
>>> from sage.all import *
>>> solns = solve([eq1,eq2,eq3,p==Integer(1)],p,q,x,y, solution_dict=True)
>>> [[s[p].n(Integer(30)), s[q].n(Integer(30)), s[x].n(Integer(30)), s[y].n(Integer(30))] for s in solns]
[[1.0000000, 8.0000000, -4.8830369, -0.13962039],
[1.0000000, 8.0000000, 3.5497035, -1.1937129]]
solns = solve([eq1,eq2,eq3,p==1],p,q,x,y, solution_dict=True) [[s[p].n(30), s[q].n(30), s[x].n(30), s[y].n(30)] for s in solns]
(The function n
prints a numerical approximation, and the
argument is the number of bits of precision.)
Solving Equations Numerically¶
Often times, solve
will not be able to find an exact solution to
the equation or equations specified. When it fails, you can use
find_root
to find a numerical solution. For example, solve does
not return anything interesting for the following equation:
sage: theta = var('theta')
sage: solve(cos(theta)==sin(theta), theta)
[sin(theta) == cos(theta)]
>>> from sage.all import *
>>> theta = var('theta')
>>> solve(cos(theta)==sin(theta), theta)
[sin(theta) == cos(theta)]
theta = var('theta') solve(cos(theta)==sin(theta), theta)
On the other hand, we can use find_root
to find a solution to the
above equation in the range
sage: phi = var('phi')
sage: find_root(cos(phi)==sin(phi),0,pi/2)
0.785398163397448...
>>> from sage.all import *
>>> phi = var('phi')
>>> find_root(cos(phi)==sin(phi),Integer(0),pi/Integer(2))
0.785398163397448...
phi = var('phi') find_root(cos(phi)==sin(phi),0,pi/2)
Differentiation, Integration, etc.¶
Sage knows how to differentiate and integrate many functions. For
example, to differentiate
sage: u = var('u')
sage: diff(sin(u), u)
cos(u)
>>> from sage.all import *
>>> u = var('u')
>>> diff(sin(u), u)
cos(u)
u = var('u') diff(sin(u), u)
To compute the fourth derivative of
sage: diff(sin(x^2), x, 4)
16*x^4*sin(x^2) - 48*x^2*cos(x^2) - 12*sin(x^2)
>>> from sage.all import *
>>> diff(sin(x**Integer(2)), x, Integer(4))
16*x^4*sin(x^2) - 48*x^2*cos(x^2) - 12*sin(x^2)
diff(sin(x^2), x, 4)
To compute the partial derivatives of
sage: x, y = var('x,y')
sage: f = x^2 + 17*y^2
sage: f.diff(x)
2*x
sage: f.diff(y)
34*y
>>> from sage.all import *
>>> x, y = var('x,y')
>>> f = x**Integer(2) + Integer(17)*y**Integer(2)
>>> f.diff(x)
2*x
>>> f.diff(y)
34*y
x, y = var('x,y') f = x^2 + 17*y^2 f.diff(x) f.diff(y)
We move on to integrals, both indefinite and definite. To compute
sage: integral(x*sin(x^2), x)
-1/2*cos(x^2)
sage: integral(x/(x^2+1), x, 0, 1)
1/2*log(2)
>>> from sage.all import *
>>> integral(x*sin(x**Integer(2)), x)
-1/2*cos(x^2)
>>> integral(x/(x**Integer(2)+Integer(1)), x, Integer(0), Integer(1))
1/2*log(2)
integral(x*sin(x^2), x) integral(x/(x^2+1), x, 0, 1)
To compute the partial fraction decomposition of
sage: f = 1/((1+x)*(x-1))
sage: f.partial_fraction(x)
-1/2/(x + 1) + 1/2/(x - 1)
>>> from sage.all import *
>>> f = Integer(1)/((Integer(1)+x)*(x-Integer(1)))
>>> f.partial_fraction(x)
-1/2/(x + 1) + 1/2/(x - 1)
f = 1/((1+x)*(x-1)) f.partial_fraction(x)
Solving Differential Equations¶
You can use Sage to investigate ordinary differential equations. To
solve the equation
sage: t = var('t') # define a variable t
sage: x = function('x')(t) # define x to be a function of that variable
sage: DE = diff(x, t) + x - 1
sage: desolve(DE, [x,t])
(_C + e^t)*e^(-t)
>>> from sage.all import *
>>> t = var('t') # define a variable t
>>> x = function('x')(t) # define x to be a function of that variable
>>> DE = diff(x, t) + x - Integer(1)
>>> desolve(DE, [x,t])
(_C + e^t)*e^(-t)
t = var('t') # define a variable t x = function('x')(t) # define x to be a function of that variable DE = diff(x, t) + x - 1 desolve(DE, [x,t])
This uses Sage’s interface to Maxima [Max], and so its output may be
a bit different from other Sage output. In this case, this says
that the general solution to the differential equation is
You can compute Laplace transforms also; the Laplace transform of
sage: s = var("s")
sage: t = var("t")
sage: f = t^2*exp(t) - sin(t)
sage: f.laplace(t,s)
-1/(s^2 + 1) + 2/(s - 1)^3
>>> from sage.all import *
>>> s = var("s")
>>> t = var("t")
>>> f = t**Integer(2)*exp(t) - sin(t)
>>> f.laplace(t,s)
-1/(s^2 + 1) + 2/(s - 1)^3
s = var("s") t = var("t") f = t^2*exp(t) - sin(t) f.laplace(t,s)
Here is a more involved example. The displacement from equilibrium (respectively) for a coupled spring attached to a wall on the left
|------\/\/\/\/\---|mass1|----\/\/\/\/\/----|mass2|
spring1 spring2
is modeled by the system of 2nd order differential equations
where
Example: Use Sage to solve the above problem with
Solution: Take the Laplace transform of the first equation (with
the notation
sage: t,s = SR.var('t,s')
sage: x = function('x')
sage: y = function('y')
sage: f = 2*x(t).diff(t,2) + 6*x(t) - 2*y(t)
sage: f.laplace(t,s)
2*s^2*laplace(x(t), t, s) - 2*s*x(0) + 6*laplace(x(t), t, s) - 2*laplace(y(t), t, s) - 2*D[0](x)(0)
>>> from sage.all import *
>>> t,s = SR.var('t,s')
>>> x = function('x')
>>> y = function('y')
>>> f = Integer(2)*x(t).diff(t,Integer(2)) + Integer(6)*x(t) - Integer(2)*y(t)
>>> f.laplace(t,s)
2*s^2*laplace(x(t), t, s) - 2*s*x(0) + 6*laplace(x(t), t, s) - 2*laplace(y(t), t, s) - 2*D[0](x)(0)
t,s = SR.var('t,s') x = function('x') y = function('y') f = 2*x(t).diff(t,2) + 6*x(t) - 2*y(t) f.laplace(t,s)
This is hard to read, but it says that
(where the Laplace transform of a lower case function like
sage: de2 = maxima("diff(y(t),t, 2) + 2*y(t) - 2*x(t)")
sage: lde2 = de2.laplace("t","s"); lde2.sage()
s^2*laplace(y(t), t, s) - s*y(0) - 2*laplace(x(t), t, s) + 2*laplace(y(t), t, s) - D[0](y)(0)
>>> from sage.all import *
>>> de2 = maxima("diff(y(t),t, 2) + 2*y(t) - 2*x(t)")
>>> lde2 = de2.laplace("t","s"); lde2.sage()
s^2*laplace(y(t), t, s) - s*y(0) - 2*laplace(x(t), t, s) + 2*laplace(y(t), t, s) - D[0](y)(0)
de2 = maxima("diff(y(t),t, 2) + 2*y(t) - 2*x(t)") lde2 = de2.laplace("t","s"); lde2.sage()
This says
Plug in the initial conditions for
sage: var('s X Y')
(s, X, Y)
sage: eqns = [(2*s^2+6)*X-2*Y == 6*s, -2*X +(s^2+2)*Y == 3*s]
sage: solve(eqns, X,Y)
[[X == 3*(s^3 + 3*s)/(s^4 + 5*s^2 + 4),
Y == 3*(s^3 + 5*s)/(s^4 + 5*s^2 + 4)]]
>>> from sage.all import *
>>> var('s X Y')
(s, X, Y)
>>> eqns = [(Integer(2)*s**Integer(2)+Integer(6))*X-Integer(2)*Y == Integer(6)*s, -Integer(2)*X +(s**Integer(2)+Integer(2))*Y == Integer(3)*s]
>>> solve(eqns, X,Y)
[[X == 3*(s^3 + 3*s)/(s^4 + 5*s^2 + 4),
Y == 3*(s^3 + 5*s)/(s^4 + 5*s^2 + 4)]]
var('s X Y') eqns = [(2*s^2+6)*X-2*Y == 6*s, -2*X +(s^2+2)*Y == 3*s] solve(eqns, X,Y)
Now take inverse Laplace transforms to get the answer:
sage: var('s t')
(s, t)
sage: inverse_laplace((3*s^3 + 9*s)/(s^4 + 5*s^2 + 4),s,t)
cos(2*t) + 2*cos(t)
sage: inverse_laplace((3*s^3 + 15*s)/(s^4 + 5*s^2 + 4),s,t)
-cos(2*t) + 4*cos(t)
>>> from sage.all import *
>>> var('s t')
(s, t)
>>> inverse_laplace((Integer(3)*s**Integer(3) + Integer(9)*s)/(s**Integer(4) + Integer(5)*s**Integer(2) + Integer(4)),s,t)
cos(2*t) + 2*cos(t)
>>> inverse_laplace((Integer(3)*s**Integer(3) + Integer(15)*s)/(s**Integer(4) + Integer(5)*s**Integer(2) + Integer(4)),s,t)
-cos(2*t) + 4*cos(t)
var('s t') inverse_laplace((3*s^3 + 9*s)/(s^4 + 5*s^2 + 4),s,t) inverse_laplace((3*s^3 + 15*s)/(s^4 + 5*s^2 + 4),s,t)
Therefore, the solution is
This can be plotted parametrically using
sage: t = var('t')
sage: P = parametric_plot((cos(2*t) + 2*cos(t), 4*cos(t) - cos(2*t) ),
....: (t, 0, 2*pi), rgbcolor=hue(0.9))
sage: show(P)
>>> from sage.all import *
>>> t = var('t')
>>> P = parametric_plot((cos(Integer(2)*t) + Integer(2)*cos(t), Integer(4)*cos(t) - cos(Integer(2)*t) ),
... (t, Integer(0), Integer(2)*pi), rgbcolor=hue(RealNumber('0.9')))
>>> show(P)
t = var('t') P = parametric_plot((cos(2*t) + 2*cos(t), 4*cos(t) - cos(2*t) ), (t, 0, 2*pi), rgbcolor=hue(0.9)) show(P)
The individual components can be plotted using
sage: t = var('t')
sage: p1 = plot(cos(2*t) + 2*cos(t), (t,0, 2*pi), rgbcolor=hue(0.3))
sage: p2 = plot(4*cos(t) - cos(2*t), (t,0, 2*pi), rgbcolor=hue(0.6))
sage: show(p1 + p2)
>>> from sage.all import *
>>> t = var('t')
>>> p1 = plot(cos(Integer(2)*t) + Integer(2)*cos(t), (t,Integer(0), Integer(2)*pi), rgbcolor=hue(RealNumber('0.3')))
>>> p2 = plot(Integer(4)*cos(t) - cos(Integer(2)*t), (t,Integer(0), Integer(2)*pi), rgbcolor=hue(RealNumber('0.6')))
>>> show(p1 + p2)
t = var('t') p1 = plot(cos(2*t) + 2*cos(t), (t,0, 2*pi), rgbcolor=hue(0.3)) p2 = plot(4*cos(t) - cos(2*t), (t,0, 2*pi), rgbcolor=hue(0.6)) show(p1 + p2)
For more on plotting, see Plotting. See section 5.5 of [NagleEtAl2004] for further information on differential equations.
Euler’s Method for Systems of Differential Equations¶
In the next example, we will illustrate Euler’s method for first and second order ODEs. We first recall the basic idea for first order equations. Given an initial value problem of the form
we want to find the approximate value of the solution at
Recall from the definition of the derivative that
where
If we call
If we break the interval from
… |
||
… |
||
… |
||
??? |
… |
The goal is to fill out all the blanks of the table, one row at a
time, until we reach the ??? entry, which is the
Euler’s method approximation for
The idea for systems of ODEs is similar.
Example: Numerically approximate
We must reduce the 2nd order ODE down to a system of two first
order DEs (using
sage: t,x,y = PolynomialRing(RealField(10),3,"txy").gens()
sage: f = y; g = -x - y * t
sage: eulers_method_2x2(f,g, 0, 1, 0, 1/4, 1)
t x h*f(t,x,y) y h*g(t,x,y)
0 1 0.00 0 -0.25
1/4 1.0 -0.062 -0.25 -0.23
1/2 0.94 -0.12 -0.48 -0.17
3/4 0.82 -0.16 -0.66 -0.081
1 0.65 -0.18 -0.74 0.022
>>> from sage.all import *
>>> t,x,y = PolynomialRing(RealField(Integer(10)),Integer(3),"txy").gens()
>>> f = y; g = -x - y * t
>>> eulers_method_2x2(f,g, Integer(0), Integer(1), Integer(0), Integer(1)/Integer(4), Integer(1))
t x h*f(t,x,y) y h*g(t,x,y)
0 1 0.00 0 -0.25
1/4 1.0 -0.062 -0.25 -0.23
1/2 0.94 -0.12 -0.48 -0.17
3/4 0.82 -0.16 -0.66 -0.081
1 0.65 -0.18 -0.74 0.022
t,x,y = PolynomialRing(RealField(10),3,"txy").gens() f = y; g = -x - y * t eulers_method_2x2(f,g, 0, 1, 0, 1/4, 1)
Therefore,
We can also plot the points eulers_method_2x2_plot
will
do this; in order to use it, we need to define functions
sage: f = lambda z: z[2] # f(t,x,y) = y
sage: g = lambda z: -sin(z[1]) # g(t,x,y) = -sin(x)
sage: P = eulers_method_2x2_plot(f,g, 0.0, 0.75, 0.0, 0.1, 1.0)
>>> from sage.all import *
>>> f = lambda z: z[Integer(2)] # f(t,x,y) = y
>>> g = lambda z: -sin(z[Integer(1)]) # g(t,x,y) = -sin(x)
>>> P = eulers_method_2x2_plot(f,g, RealNumber('0.0'), RealNumber('0.75'), RealNumber('0.0'), RealNumber('0.1'), RealNumber('1.0'))
f = lambda z: z[2] # f(t,x,y) = y g = lambda z: -sin(z[1]) # g(t,x,y) = -sin(x) P = eulers_method_2x2_plot(f,g, 0.0, 0.75, 0.0, 0.1, 1.0)
At this point, P
is storing two plots: P[0]
, the plot of P[1]
, the plot of
sage: show(P[0] + P[1])
>>> from sage.all import *
>>> show(P[Integer(0)] + P[Integer(1)])
show(P[0] + P[1])
(For more on plotting, see Plotting.)
Special functions¶
Several orthogonal polynomials and special functions are implemented, using both PARI [GAP] and Maxima [Max]. These are documented in the appropriate sections (“Orthogonal polynomials” and “Special functions”, respectively) of the Sage reference manual.
sage: x = polygen(QQ, 'x')
sage: chebyshev_U(2,x)
4*x^2 - 1
sage: bessel_I(1,1).n(250)
0.56515910399248502720769602760986330732889962162109200948029448947925564096
sage: bessel_I(1,1).n()
0.565159103992485
sage: bessel_I(2,1.1).n()
0.167089499251049
>>> from sage.all import *
>>> x = polygen(QQ, 'x')
>>> chebyshev_U(Integer(2),x)
4*x^2 - 1
>>> bessel_I(Integer(1),Integer(1)).n(Integer(250))
0.56515910399248502720769602760986330732889962162109200948029448947925564096
>>> bessel_I(Integer(1),Integer(1)).n()
0.565159103992485
>>> bessel_I(Integer(2),RealNumber('1.1')).n()
0.167089499251049
x = polygen(QQ, 'x') chebyshev_U(2,x) bessel_I(1,1).n(250) bessel_I(1,1).n() bessel_I(2,1.1).n()
At this point, Sage has only wrapped these functions for numerical use. For symbolic use, please use the Maxima interface directly, as in the following example:
sage: maxima.eval("f:bessel_y(v, w)")
'bessel_y(v,w)'
sage: maxima.eval("diff(f,w)")
'(bessel_y(v-1,w)-bessel_y(v+1,w))/2'
>>> from sage.all import *
>>> maxima.eval("f:bessel_y(v, w)")
'bessel_y(v,w)'
>>> maxima.eval("diff(f,w)")
'(bessel_y(v-1,w)-bessel_y(v+1,w))/2'
maxima.eval("f:bessel_y(v, w)") maxima.eval("diff(f,w)")
Vector calculus¶
See the Vector Calculus Tutorial.