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

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

Python

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

Sage Live

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

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

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

Python

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

Sage Live

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, \(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

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)

Python

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

Sage Live

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 \(x\) vs. \(t\), and P[1], the plot of \(y\) vs. \(t\). We can plot both of these as follows:

Sage

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

Python

>>> from sage.all import *
>>> show(P[Integer(0)] + P[Integer(1)])

Sage Live

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

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

Python

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

Sage Live

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

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'

Python

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

Sage Live

maxima.eval("f:bessel_y(v, w)")
maxima.eval("diff(f,w)")

Vector calculus¶

See the Vector Calculus Tutorial.

Basic Algebra and Calculus¶

Solving Equations¶

Solving Equations Exactly¶

Solving Equations Numerically¶

Differentiation, Integration, etc.¶

Solving Differential Equations¶

Euler’s Method for Systems of Differential Equations¶

Special functions¶

Vector calculus¶