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)

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]

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)]

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]]

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]]

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]]

(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)]

On the other hand, we can use find_root to find a solution to the above equation in the range \(0 < \phi < \pi/2\):

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...

Differentiation, Integration, etc.#

Sage knows how to differentiate and integrate many functions. For example, to differentiate \(\sin(u)\) with respect to \(u\), do the following:

sage: u = var('u')
sage: diff(sin(u), u)
cos(u)

>>> from sage.all import *
>>> u = var('u')
>>> diff(sin(u), u)
cos(u)

To compute the fourth derivative of \(\sin(x^2)\):

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)

To compute the partial derivatives of \(x^2+17y^2\) with respect to \(x\) and \(y\), respectively:

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

We move on to integrals, both indefinite and definite. To compute \(\int x\sin(x^2)\, dx\) and \(\int_0^1 \frac{x}{x^2+1}\, dx\)

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)

To compute the partial fraction decomposition of \(\frac{1}{x^2-1}\):

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)

Solving Differential Equations#

You can use Sage to investigate ordinary differential equations. To solve the equation \(x'+x-1=0\):

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)

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 \(x(t) = e^{-t}(e^{t}+c)\).

You can compute Laplace transforms also; the Laplace transform of \(t^2e^t -\sin(t)\) is computed as follows:

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

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

\[ \begin{align}\begin{aligned}m_1 x_1'' + (k_1+k_2) x_1 - k_2 x_2 = 0\\m_2 x_2''+ k_2 (x_2-x_1) = 0,\end{aligned}\end{align} \]

where \(m_{i}\) is the mass of object i, \(x_{i}\) is the displacement from equilibrium of mass i, and \(k_{i}\) is the spring constant for spring i.

Example: Use Sage to solve the above problem with \(m_{1}=2\), \(m_{2}=1\), \(k_{1}=4\), \(k_{2}=2\), \(x_{1}(0)=3\), \(x_{1}'(0)=0\), \(x_{2}(0)=3\), \(x_{2}'(0)=0\).

Solution: Take the Laplace transform of the first equation (with the notation \(x=x_{1}\), \(y=x_{2}\)):

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)

This is hard to read, but it says that

\[-2x'(0) + 2s^2 \cdot X(s) - 2sx(0) - 2Y(s) + 6X(s) = 0\]

(where the Laplace transform of a lower case function like \(x(t)\) is the upper case function \(X(s)\)). Take the Laplace transform of the second equation:

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)

This says

\[-Y'(0) + s^2Y(s) + 2Y(s) - 2X(s) - sy(0) = 0.\]

Plug in the initial conditions for \(x(0)\), \(x'(0)\), \(y(0)\), and \(y'(0)\), and solve the resulting two equations:

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)]]

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)

Therefore, the solution is

\[x_1(t) = \cos(2t) + 2\cos(t), \quad x_2(t) = 4\cos(t) - \cos(2t).\]

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)

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)

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

\[y'=f(x,y), \quad y(a)=c,\]

we want to find the approximate value of the solution at \(x=b\) with \(b>a\).

Recall from the definition of the derivative that

\[y'(x) \approx \frac{y(x+h)-y(x)}{h},\]

where \(h>0\) is given and small. This and the DE together give \(f(x,y(x))\approx \frac{y(x+h)-y(x)}{h}\). Now solve for \(y(x+h)\):

\[y(x+h) \approx y(x) + h\cdot f(x,y(x)).\]

If we call \(h \cdot f(x,y(x))\) the “correction term” (for lack of anything better), call \(y(x)\) the “old value of \(y\)”, and call \(y(x+h)\) the “new value of \(y\)”, then this approximation can be re-expressed as

\[y_{new} \approx y_{old} + h\cdot f(x,y_{old}).\]

If we break the interval from \(a\) to \(b\) into \(n\) steps, so that \(h=\frac{b-a}{n}\), then we can record the information for this method in a table.

\(x\)

\(y\)

\(h\cdot f(x,y)\)

\(a\)

\(c\)

\(h\cdot f(a,c)\)

\(a+h\)

\(c+h\cdot f(a,c)\)

…

\(a+2h\)

…

…

\(b=a+nh\)

???

…

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 \(y(b)\).

The idea for systems of ODEs is similar.

Example: Numerically approximate \(z(t)\) at \(t=1\) using 4 steps of Euler’s method, where \(z''+tz'+z=0\), \(z(0)=1\), \(z'(0)=0\).

We must reduce the 2nd order ODE down to a system of two first order DEs (using \(x=z\), \(y=z'\)) and apply Euler’s method:

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

Therefore, \(z(1)\approx 0.65\).

We can also plot the points \((x,y)\) to get an approximate picture of the curve. The function eulers_method_2x2_plot will do this; in order to use it, we need to define functions \(f\) and \(g\) which takes one argument with three coordinates: (\(t\), \(x\), \(y\)).

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'))

At this point, P is storing two plots: P[0], the plot of \(x\) vs. \(t\), and P[1], the plot of \(y\) vs. \(t\). We can plot both of these as follows:

sage: show(P[0] + P[1])

>>> from sage.all import *
>>> show(P[Integer(0)] + P[Integer(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

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'

Vector calculus#

See the Vector Calculus Tutorial.