Navigation
- index
- modules |
- next |
- previous |
- Sage Tutorial, part 2.1 »
Now, let's enter a matrix. The function matrix() or Matrix() is used to do this. For now the
first argument to the function should be "QQ", this means that the matrix will
contain either integers or rational numbers (fractions). The second argument is
sage: M=matrix(QQ,[[2,4,0,8],[-1,3,3,-2],[0,1,1,0]]); M
Now, let's enter a matrix. The function matrix() or Matrix() is used to do this. For now the
first argument to the function should be "QQ", this means that the matrix will
contain either integers or rational numbers (fractions). The second argument is
sage: M=matrix(QQ,[[2,4,0,8],[-1,3,3,-2],[0,1,1,0]]); M
Perhaps the most confusing thing about using Sage for matrix work is that rows and columns are numbered starting at 0 rather than 1 as is usually done for matrices in mathematical work. This means that the rows in M are numbered 0, 1, and 2 while the columns are numbered 0, 1, 2, and 3. For example, we can access the first and second rows with:
sage: M.row(0); M.row(1)
Notice how the function row() is used; it is "attached" to the matrix varible with a dot. This means that the row function operates on the matrix M.
sage: A = matrix(QQ,3,3,[2,-1,1,1,-3,3,5,-2,-3])
sage: A
[ 2 -1 1]
[ 1 -3 3]
[ 5 -2 -3]
sage: matrix([[2,-1,1],[1,-3,3],[5,-2,-3]])
[ 2 -1 1]
[ 1 -3 3]
[ 5 -2 -3]
sage: Matrix([[2,-1,1],[1,-3,3],[5,-2,-3]])
[ 2 -1 1]
[ 1 -3 3]
[ 5 -2 -3]
Suppose M is the augmented matrix of a linear system. We can solve the linear system by performing elementary row operations on M. In Sage these row operations are implemented with the functions.
sage: Mswap_rows() # interchange two rows
sage: rescale_row() # scale a row by using a scale factor
sage: add_multiple_of_row() # add a multiple of one row to another row, replacing the row
Remember, row numbering starts at 0. Pay careful attention to the changes made to the matrix by each of the following commands.
We want to start with a 1 in the first position of the first column, so we begin by scaling the first row by 1/2.
sage: M.rescale_row(0,1/2); M
The first argument to rescale_row() is the index of the row to be scaled and the second argument is the scale factor. We could, of course, use 0.5 rather than 1/2 for the scale factor.
Now that there is a 1 in the first position of the first column we continue by eliminating the entry below it.
sage: M.add_multiple_of_row(1,0,1); M
The first argument is the index of the row to be replaced, the second argument is the row to form a multiple of, and the final argument is the scale factor. Thus, M.add_multiple_of_row(n,m,a) would replace row n with (row n)+a*(row m).
Since the last entry in the first column is already zero we can move on to the second column. Our first step is to get a 1 in the second position of the second column. Normally we would do this by scaling the row but in this case we can swap the second and third row
sage: M.swap_rows(1,2); M
The arguments to swap_rows() are fairly obvious, just remember that row 1 is the second row and row 2 is the third row.
Now we want to eliminate the 5 below the 1. This is done by multiplying the second row by −5 and adding it to the third row.
sage: M.add_multiple_of_row(2,1,-5); M
sage: M.rescale_row(2,-1/2); M
sage: M.add_multiple_of_row(1,2,-1); M
sage: M.add_multiple_of_row(0,1,-2); M
At this point our solution is complete. We see that the solution to our linear system is x1=2, x2=1, and x3=−1.
There is an easy way to check your work, or to carry out these steps in the future. First, let's reload the matrix M.
sage: M=matrix(QQ,[[2,4,0,8],[-1,3,3,-2],[0,1,1,0]]); M
sage: M.echelon_form()
sage: M
sage: M.echelonize(); M
There are several ways to solve system of linear algebraic equations:
\( {\bf A} {\bf x} = {\bf b} \)
It can be done by reduction of augmented matrix
sage: A = matrix(QQ,3,3,[2,-1,1,1,-3,3,5,-2,-3])
sage: b = vector((0,1,-1))
sage: C = A.augment(b)
sage: C.rref()
Another way is to use default commands:
sage: A = matrix(QQ,3,3,[2,-1,1,1,-3,3,5,-2,-3])
sage: b= vector([1,1,-4])
sage: A \ b
(2/5, 27/25, 32/25)
sage: x = A.solve_right(b)
sage: x
(2/5, 27/25, 32/25)
We show how to define a function of a square matrix using diagonalization procedure. This method is applicable only for such matrices, and is not suatable for defective matrices. Recall that a matrix A is called diagonalizable if there exists a nonsingular matrix S such that \( {\bf S}^{-1} {\bf A} {\bf S} = {\bf D} , \) a diagonal matrix. In other words, the matrix A is similar to a diagonal matrix. An \( n \times n \) square matrix is diagonalizable if and only if there exist n linearly independent eigenvectors, so geometrical multiplicities of each eigenvalue are the same as its algebraic multiplicities. Then the matrix S can be built from eigenvectors of A, column by column.
Let A be a square \( n \times n \) diagonalizable matrix, and let D be the corresponding diagonal matrix of its eigenvalues:
where \( \lambda_1 , \lambda_2 , \ldots , \lambda_n \) are eigenvalues (that may be equal) of the matrix A.
Let \( {\bf x}_1 , {\bf x}_2 , \ldots , {\bf x}_n \) be linearly independent eigenvectors, corresponding to the eigenvalues \( \lambda_1 , \lambda_2 , \ldots , \lambda_n .\) We build the nonsingular matrix S from these eigenvectors (every column is an eigenvector):
Numerical integration methods can prove useful if the integrand is known only at certain points or the antiderivate is very difficult or even mpossible to find. In education, calculating numerical methods by hand may be useful in some cases, but computer programs are usually better suited in finding patterns and comparing different methods. In the next example, three numerical integration methods are implemented in Sage: the midpoint rule, the trapezoidal rule and Simpson's rule. The differences between the exact value of integration and the approximation are tabulated by the number of subintervals n
sage: f(x) = x^2-3
sage: a = 0.0
sage: b = 2.0
sage: table = []
sage: exact = integrate(f(x), x, a, b)
sage: for n in [4, 10, 20, 50, 100]:
....: h = (b-a)/n
....: midpoint = sum([f(a+(i+1/2)*h)*h for i in range(n)])
....: trapezoid = h/2*(f(a) + 2*sum([f(a+i*h) for i in range(1,n)])
....: + f(b))
....: simpson = h/3*(f(a) + sum([4*f(a+i*h) for i in range(1,n,2)])
....: + sum([2*f(a+i*h) for i in range (2,n,2)]) + f(b))
....: table.append([n, h.n(digits=2), (midpoint-exact).n(digits=6),
....: (trapezoid-exact).n(digits=6), (simpson-exact).n(digits=6)])
....: html.table(table, header=["n", "h", "Midpoint rule",
....: "Trapezoidal rule", "Simpson's rule"])
There are also built-in methods for numerical integration in Sage.
For instance, it is possible to automatically produce piecewise-defined
line functions defined by the trapezoidal rule or the midpoint rule.
These functions can be used to visualize different geometric interpretations of the numerical integration methods. In the next example,
midpoint rule is used to calculate an approximation for the definite
integral of the function
f(x) =x^2 -5x + 10
over the interval [0;10] using six subintervals
sage: f(x) = x^2-5*x+10
sage: f = Piecewise([[(0,10), f]])
sage: g = f.riemann_sum(6, mode="midpoint")
sage: F = f.plot(color="blue")
sage: R = add([line([[a,0],[a,f(x=a)],[b,f(x=b)],[b,0]], color="red") for (a,b), f in g.list()]
sage: show(F+R)
The Trapezoid Rule approximates the area under a given curve by finding the area under a linear approximation to the curve. The linear segments are given by the secant lines through the endpoints of the subintervals: for f(x)=x^2+1 on [0,2] with 4 subintervals, this looks like:
sage: x = var('x')
sage: f1(x) = x^2 + 1
sage: f = Piecewise([[(0,2),f1]])
sage: trapezoid_sum = f.trapezoid(4)
sage: P = f.plot(rgbcolor=(0,0,1), plot_points=40)
sage: Q = trapezoid_sum.plot(rgbcolor=(1,0,0), plot_points=40)
sage: L = add([line([[a,0],[a,f(x=a)]], rgbcolor=(1,0,0)) for (a,b), f in trapezoid_sum.list()])
sage: M = line([[2,0],[2,f1(2)]], rgbcolor=(1,0,0))
sage: show(P + Q + L + M)
To integrate the function cos(2*x)^4 +sin(x) we do
sage: numerical_integral(lambda x: cos(2*x)^5 + sin(x), 0, 2*pi)
(2.3561944901923457, 5.322071926740605e-14)
The input can be any callable:
sage: numerical_integral(lambda x: cos(2*x)^5 + sin(x), 0, 2*pi)
(2.3561944901923457, 5.322071926740605e-14)
We check this with a symbolic integration:
sage: (cos(2*x)^4+sin(x)).integral(x,0,2*pi)
3/4*pi
It is possible to integrate on infinite intervals as well by using +Infinity or -Infinity in the interval argument. For example:
sage: f = exp(-2*x)
sage: numerical_integral(f, 0, +Infinity)
(0.4999999999993456, 2.5048509541974486e-07)
We can also numerically integrate symbolic expressions using either this function (which uses GSL) or the native integration (which uses Maxima):
sage: (x^3*sin(1/x)).nintegral(x,0,pi)
(9.874626194990983, 6.585622218041451e-08, 315, 0)
Regarding numerical approximation of \(\int_a^bf(x)\, dx\), where \(f\) is a piecewise defined function, can
sage: f1(x) = x^2
sage: f2(x) = 5-x^2
sage: f = Piecewise([[(0,1),f1],[(1,2),f2]])
sage: f.trapezoid(4)
Piecewise defined function with 4 parts, [[(0, 1/2), 1/2*x],
[(1/2, 1), 9/2*x - 2], [(1, 3/2), 1/2*x + 2],
[(3/2, 2), -7/2*x + 8]]
sage: f.riemann_sum_integral_approximation(6,mode="right")
19/6
sage: f.integral()
Piecewise defined function with 2 parts,
[[(0, 1), x |--> 1/3*x^3], [(1, 2), x |--> -1/3*x^3 + 5*x - 13/3]]
sage: f.integral(definite=True)
3
Differentiation:
sage: var('x k w')
(x, k, w)
sage: f = x^3 * e^(k*x) * sin(w*x); f
x^3*e^(k*x)*sin(w*x)
sage: f.diff(x)
w*x^3*cos(w*x)*e^(k*x) + k*x^3*e^(k*x)*sin(w*x) + 3*x^2*e^(k*x)*sin(w*x)
sage: latex(f.diff(x))
w x^{3} \cos\left(w x\right) e^{\left(k x\right)} + k x^{3} e^{\left(k x\right)} \sin\left(w x\right) + 3 \, x^{2} e^{\left(k x\right)} \sin\left(w x\right)
If you type view(f.diff(x)) another window will open up displaying the compiled output. In the notebook, you can enter
var('x k w')
f = x^3 * e^(k*x) * sin(w*x)
show(f)
show(f.diff(x))
into a cell and press shift-enter for a similar result. You can also differentiate and integrate using the commands
R = PolynomialRing(QQ,"x")
x = R.gen()
p = x^2 + 1
show(p.derivative())
show(p.integral())
in a notebook cell, or
sage: R = PolynomialRing(QQ,"x")
sage: x = R.gen()
sage: p = x^2 + 1
sage: p.derivative()
2*x
sage: p.integral()
1/3*x^3 + x
on the command line. At this point you can also type view(p.derivative()) or view(p.integral()) to open a new window with output typeset by LaTeX.
You can find critical points of a piecewise defined function:
sage: x = PolynomialRing(RationalField(), 'x').gen()
sage: f1 = x^0
sage: f2 = 1-x
sage: f3 = 2*x
sage: f4 = 10*x-x^2
sage: f = Piecewise([[(0,1),f1],[(1,2),f2],[(2,3),f3],[(3,10),f4]])
sage: f.critical_points()
[5.0]
Taylor series:
sage: var('f0 k x')
(f0, k, x)
sage: g = f0/sinh(k*x)^4
sage: g.taylor(x, 0, 3)
-62/945*f0*k^2*x^2 + 11/45*f0 - 2/3*f0/(k^2*x^2) + f0/(k^4*x^4)
sage: maxima(g).powerseries('_SAGE_VAR_x',0) # TODO: write this without maxima
16*_SAGE_VAR_f0*('sum((2^(2*i1-1)-1)*bern(2*i1)*_SAGE_VAR_k^(2*i1-1)*_SAGE_VAR_x^(2*i1-1)/factorial(2*i1),i1,0,inf))^4
Of course, you can view the LaTeX-ed version of this using view(g.powerseries('x',0)).
The Maclaurin and power series of \(\log({\frac{\sin(x)}{x}})\):
sage: f = log(sin(x)/x)
sage: f.taylor(x, 0, 10)
-1/467775*x^10 - 1/37800*x^8 - 1/2835*x^6 - 1/180*x^4 - 1/6*x^2
sage: [bernoulli(2*i) for i in range(1,7)]
[1/6, -1/30, 1/42, -1/30, 5/66, -691/2730]
sage: maxima(f).powerseries(x,0) # TODO: write this without maxima
'sum((-1)^i3*2^(2*i3-1)*bern(2*i3)*_SAGE_VAR_x^(2*i3)/(i3*factorial(2*i3)),i3,1,inf)
Numerical integration is discussed in Riemann and trapezoid sums for integrals below.
Sage can integrate some simple functions on its own:
sage: f = x^3
sage: f.integral(x)
1/4*x^4
sage: integral(x^3,x)
1/4*x^4
sage: f = x*sin(x^2)
sage: integral(f,x)
-1/2*cos(x^2)
Sage can also compute symbolic definite integrals involving limits.
sage: var('x, k, w')
(x, k, w)
sage: f = x^3 * e^(k*x) * sin(w*x)
sage: f.integrate(x)
((24*k^3*w - 24*k*w^3 - (k^6*w + 3*k^4*w^3 + 3*k^2*w^5 + w^7)*x^3 + 6*(k^5*w + 2*k^3*w^3 + k*w^5)*x^2 - 6*(3*k^4*w + 2*k^2*w^3 - w^5)*x)*cos(w*x)*e^(k*x) - (6*k^4 - 36*k^2*w^2 + 6*w^4 - (k^7 + 3*k^5*w^2 + 3*k^3*w^4 + k*w^6)*x^3 + 3*(k^6 + k^4*w^2 - k^2*w^4 - w^6)*x^2 - 6*(k^5 - 2*k^3*w^2 - 3*k*w^4)*x)*e^(k*x)*sin(w*x))/(k^8 + 4*k^6*w^2 + 6*k^4*w^4 + 4*k^2*w^6 + w^8)
sage: integrate(1/x^2, x, 1, infinity)
1
You can find the convolution of any piecewise defined function with another (off the domain of definition, they are assumed to be zero). Here is \(f\), \(f*f\), and \(f*f*f\), where \(f(x)=1\), \(0<x<1\):
sage: x = PolynomialRing(QQ, 'x').gen()
sage: f = Piecewise([[(0,1),1*x^0]])
sage: g = f.convolution(f)
sage: h = f.convolution(g)
sage: P = f.plot(); Q = g.plot(rgbcolor=(1,1,0)); R = h.plot(rgbcolor=(0,1,1))
To view this, type show(P+Q+R).
If you have a piecewise-defined polynomial function then there is a “native” command for computing Laplace transforms. This calls Maxima but it’s worth noting that Maxima cannot handle (using the direct interface illustrated in the last few examples) this type of computation.
sage: var('x s')
(x, s)
sage: f1(x) = 1
sage: f2(x) = 1-x
sage: f = Piecewise([[(0,1),f1],[(1,2),f2]])
sage: f.laplace(x, s)
-e^(-s)/s + (s + 1)*e^(-2*s)/s^2 + 1/s - e^(-s)/s^2
For other “reasonable” functions, Laplace transforms can be computed using the Maxima interface:
sage: var('k, s, t')
(k, s, t)
sage: f = 1/exp(k*t)
sage: f.laplace(t,s)
1/(k + s)
is one way to compute LT’s and
sage: var('s, t')
(s, t)
sage: f = t^5*exp(t)*sin(t)
sage: L = laplace(f, t, s); L
3840*(s - 1)^5/(s^2 - 2*s + 2)^6 - 3840*(s - 1)^3/(s^2 - 2*s + 2)^5 +
720*(s - 1)/(s^2 - 2*s + 2)^4
is another way.
Symbolically solving ODEs can be done using Sage interface with Maxima. See
sage:desolvers?
for available commands. Numerical solution of ODEs can be done using Sage interface with Octave (an experimental package), or routines in the GSL (Gnu Scientific Library).
An example, how to solve ODE’s symbolically in Sage using the Maxima interface (do not type the ...):
sage: y=function('y',x); desolve(diff(y,x,2) + 3*x == y, dvar = y, ics = [1,1,1])
3*x - 2*e^(x - 1)
sage: desolve(diff(y,x,2) + 3*x == y, dvar = y)
_K2*e^(-x) + _K1*e^x + 3*x
sage: desolve(diff(y,x) + 3*x == y, dvar = y)
(3*(x + 1)*e^(-x) + _C)*e^x
sage: desolve(diff(y,x) + 3*x == y, dvar = y, ics = [1,1]).expand()
3*x - 5*e^(x - 1) + 3
sage: f=function('f',x); desolve_laplace(diff(f,x,2) == 2*diff(f,x)-f, dvar = f, ics = [0,1,2])
x*e^x + e^x
sage: desolve_laplace(diff(f,x,2) == 2*diff(f,x)-f, dvar = f)
-x*e^x*f(0) + x*e^x*D[0](f)(0) + e^x*f(0)
If you have Octave and gnuplot installed,
sage: octave.de_system_plot(['x+y','x-y'], [1,-1], [0,2]) # optional - octave
yields the two plots \((t,x(t)), (t,y(t))\) on the same graph (the \(t\)-axis is the horizonal axis) of the system of ODEs
for \(0 <= t <= 2\). The same result can be obtained by using desolve_system_rk4:
sage: x, y, t = var('x y t')
sage: P=desolve_system_rk4([x+y, x-y], [x,y], ics=[0,1,-1], ivar=t, end_points=2)
sage: p1 = list_plot([[i,j] for i,j,k in P], plotjoined=True)
sage: p2 = list_plot([[i,k] for i,j,k in P], plotjoined=True, color='red')
sage: p1+p2
Graphics object consisting of 2 graphics primitives
Another way this system can be solved is to use the command desolve_system.
sage: t=var('t'); x=function('x',t); y=function('y',t)
sage: des = [diff(x,t) == x+y, diff(y,t) == x-y]
sage: desolve_system(des, [x,y], ics = [0, 1, -1])
[x(t) == cosh(sqrt(2)*t), y(t) == sqrt(2)*sinh(sqrt(2)*t) - cosh(sqrt(2)*t)]
The output of this command is not a pair of functions.
Finally, can solve linear DEs using power series:
sage: R.<t> = PowerSeriesRing(QQ, default_prec=10)
sage: a = 2 - 3*t + 4*t^2 + O(t^10)
sage: b = 3 - 4*t^2 + O(t^7)
sage: f = a.solve_linear_de(prec=5, b=b, f0=3/5)
sage: f
3/5 + 21/5*t + 33/10*t^2 - 38/15*t^3 + 11/24*t^4 + O(t^5)
sage: f.derivative() - a*f - b
O(t^4)
If \(f(x)\) is a piecewise-defined polynomial function on \(-L<x<L\) then the Fourier series
converges. In addition to computing the coefficients \(a_n,b_n\), it will also compute the partial sums (as a string), plot the partial sums (as a function of \(x\) over \((-L,L)\), for comparison with the plot of \(f(x)\) itself), compute the value of the FS at a point, and similar computations for the cosine series (if \(f(x)\) is even) and the sine series (if \(f(x)\) is odd). Also, it will plot the partial F.S. Cesaro mean sums (a “smoother” partial sum illustrating how the Gibbs phenomenon is mollified).
sage: f1 = lambda x: -1
sage: f2 = lambda x: 2
sage: f = Piecewise([[(0,pi/2),f1],[(pi/2,pi),f2]])
sage: f.fourier_series_cosine_coefficient(5,pi)
-3/5/pi
sage: f.fourier_series_sine_coefficient(2,pi)
-3/pi
sage: f.fourier_series_partial_sum(3,pi)
-3*cos(x)/pi - 3*sin(2*x)/pi + sin(x)/pi + 1/4
Type show(f.plot_fourier_series_partial_sum(15,pi,-5,5)) and show(f.plot_fourier_series_partial_sum_cesaro(15,pi,-5,5)) (and be patient) to view the partial sums.