rlopez

3010 Reputation

14 Badges

19 years, 282 days

Dr. Robert J. Lopez, Emeritus Professor of Mathematics at the Rose-Hulman Institute of Technology in Terre Haute, Indiana, USA, is an award winning educator in mathematics and is the author of several books including Advanced Engineering Mathematics (Addison-Wesley 2001). For over two decades, Dr. Lopez has also been a visionary figure in the introduction of Maplesoft technology into undergraduate education. Dr. Lopez earned his Ph.D. in mathematics from Purdue University, his MS from the University of Missouri - Rolla, and his BA from Marist College. He has held academic appointments at Rose-Hulman (1985-2003), Memorial University of Newfoundland (1973-1985), and the University of Nebraska - Lincoln (1970-1973). His publication and research history includes manuscripts and papers in a variety of pure and applied mathematics topics. He has received numerous awards for outstanding scholarship and teaching.

MaplePrimes Activity


These are Posts that have been published by rlopez

This post describes how Maple was used to investigate the Givens rotation matrix, and to answer a simple question about its behavior. The "Givens" part is the medium, but the message is that it really is better to teach, learn, and do mathematics with a tool like Maple.

The question: If Givens rotations are used to take the vector Y = <5, -2, 1> to Y2 = , about what axis and through what angle will a single rotation accomplish the same thing?

The Givens matrix G21 takes Y to the vector Y1 =, and the Givens matrix G31 takes Y1 to Y2. Graphing the vectors Y, Y1, and Y2 reveals that Y1 lies in the xz-plane and that Y2 is parallel to the x-axis. (These geometrical observations should have been obvious, but the typical usage of the Givens technique to "zero-out" elements in a vector or matrix obscured this, at least for me.)

The matrix G = G31 G21 rotates Y directly to Y2; is the axis of rotation the vector W = Y x Y2, and is the angle of rotation the angle  between Y and Y2? To test these hypotheses, I used the RotationMatrix command in the Student LinearAlgebra package to build the corresponding rotation matrix R. But R did not agree with G. I had either the axis or the angle (actually both) incorrect.

The individual Givens rotation matrices are orthogonal, so G, their product is also orthogonal. It will have 1 as its single real eigenvalue, and the corresponding eigenvector V is actually the direction of the axis of the rotation. The vector W is a multiple of <0, 1, 2> but V = <a, b, 1>, where . Clearly, W  V.

The rotation matrix that rotates about the axis V through the angle  isn't the matrix G either. The correct angle of rotation about V turns out to be

the angle between the projections of Y and Y2 onto the plane orthogonal to V. That came as a great surprise, one that required a significant adjustment of my intuition about spatial rotations. So again, the message is that teaching, learning, and doing mathematics is so much more effective and efficient when done with a tool like Maple.

A discussion of the Givens rotation, and a summary of the actual computations described above are available in the attached worksheet, What Gives with Givens.mw.

I have two linear algebra texts [1, 2]  with examples of the process of constructing the transition matrix Q that brings a matrix A to its Jordan form J. In each, the authors make what seems to be arbitrary selections of basis vectors via processes that do not seem algorithmic. So recently, while looking at some other calculations in linear algebra, I decided to revisit these calculations in as orderly a way as possible.

 

First, I needed a matrix A with a prescribed Jordan form. Actually, I started with a Jordan form, and then constructed A via a similarity transform on J. To avoid introducing fractions, I sought transition matrices P with determinant 1.

 

Let's begin with J, obtained with Maple's JordanBlockMatrix command.

 

• 

Tools_Load Package: Linear Algebra

Loading LinearAlgebra

J := JordanBlockMatrix([[2, 3], [2, 2], [2, 1]])

Matrix([[2, 1, 0, 0, 0, 0], [0, 2, 1, 0, 0, 0], [0, 0, 2, 0, 0, 0], [0, 0, 0, 2, 1, 0], [0, 0, 0, 0, 2, 0], [0, 0, 0, 0, 0, 2]])

 

``

The eigenvalue lambda = 2 has algebraic multiplicity 6. There are sub-blocks of size 3×3, 2×2, and 1×1. Consequently, there will be three eigenvectors, supporting chains of generalized eigenvectors having total lengths 3, 2, and 1. Before delving further into structural theory, we next find a transition matrix P with which to fabricate A = P*J*(1/P).

 

The following code generates random 6×6 matrices of determinant 1, and with integer entries in the interval [-2, 2]. For each, the matrix A = P*J*(1/P) is computed. From these candidates, one A is then chosen.

 

L := NULL:

 

 

After several such trials, the matrix A was chosen as

 

A := Matrix(6, 6, {(1, 1) = -8, (1, 2) = -8, (1, 3) = 4, (1, 4) = -8, (1, 5) = -1, (1, 6) = 5, (2, 1) = -1, (2, 2) = 3, (2, 3) = 1, (2, 4) = -2, (2, 5) = 2, (2, 6) = -1, (3, 1) = -13, (3, 2) = -9, (3, 3) = 8, (3, 4) = -11, (3, 5) = 1, (3, 6) = 5, (4, 1) = 3, (4, 2) = 3, (4, 3) = -1, (4, 4) = 4, (4, 5) = 1, (4, 6) = -2, (5, 1) = 7, (5, 2) = 5, (5, 3) = -3, (5, 4) = 6, (5, 5) = 2, (5, 6) = -3, (6, 1) = -6, (6, 2) = -2, (6, 3) = 3, (6, 4) = -7, (6, 5) = 2, (6, 6) = 3})

 

 

for which the characteristic and minimal polynomials are

 

factor(CharacteristicPolynomial(A, lambda))

(lambda-2)^6

factor(MinimalPolynomial(A, lambda))

(lambda-2)^3

 

 

So, if we had started with just A, we'd now know that the algebraic multiplicity of its one eigenvalue lambda = 2 is 6, and there is at least one 3×3 sub-block in the Jordan form. We would not know if the other sub-blocks were all 1×1, or a 1×1 and a 2×2, or another 3×3. Here is where some additional theory must be invoked.

``

The null spaces M[k] of the matrices (A-2*I)^k are nested: `&sub;`(`&sub;`(M[1], M[2]), M[3]) .. (), as depicted in Figure 1, where the vectors a[k], k = 1, () .. (), 6, are basis vectors.

 

Figure 1   The nesting of the null spaces M[k] 

 

 

The vectors a[1], a[2], a[3] are eigenvectors, and form a basis for the eigenspace M[1]. The vectors a[k], k = 1, () .. (), 5, form a basis for the subspace M[2], and the vectors a[k], k = 1, () .. (), 6, for a basis for the space M[3], but the vectors a[4], a[5], a[6] are not yet the generalized eigenvectors. The vector a[6] must be replaced with a vector b[6] that lies in M[3] but is not in M[2]. Once such a vector is found, then a[4] can be replaced with the generalized eigenvector `&equiv;`(b[4], (A-2*I)^2)*b[6], and a[1] can be replaced with `&equiv;`(b[1], A-2*I)*b[4]. The vectors b[1], b[4], b[6] are then said to form a chain, with b[1] being the eigenvector, and b[4] and b[6] being the generalized eigenvectors.

 

If we could carry out these steps, we'd be in the state depicted in Figure 2.

 

Figure 2   The null spaces M[k] with the longest chain determined

 

 

Next, basis vector a[5] is to be replaced with b[5], a vector in M[2] but not in M[1], and linearly independent of b[4]. If such a b[5] is found, then a[2] is replaced with the generalized eigenvector `&equiv;`(b[2], A-2*I)*b[5]. The vectors b[2] and b[5] would form a second chain, with b[2] as the eigenvector, and b[5] as the generalized eigenvector.

``

Define the matrix C = A-2*I by the Maple calculation

 

C := A-2

Matrix([[-10, -8, 4, -8, -1, 5], [-1, 1, 1, -2, 2, -1], [-13, -9, 6, -11, 1, 5], [3, 3, -1, 2, 1, -2], [7, 5, -3, 6, 0, -3], [-6, -2, 3, -7, 2, 1]])

 

``

and note

 

N := convert(NullSpace(C), list)

[Vector(6, {(1) = 1/2, (2) = 1/2, (3) = 1, (4) = 0, (5) = 0, (6) = 1}), Vector(6, {(1) = -1/2, (2) = -1/2, (3) = -2, (4) = 0, (5) = 1, (6) = 0}), Vector(6, {(1) = -2, (2) = 1, (3) = -1, (4) = 1, (5) = 0, (6) = 0})]

NN := convert(LinearAlgebra:-NullSpace(C^2), list)

[Vector(6, {(1) = 2/5, (2) = 0, (3) = 0, (4) = 0, (5) = 0, (6) = 1}), Vector(6, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 1, (6) = 0}), Vector(6, {(1) = -1, (2) = 0, (3) = 0, (4) = 1, (5) = 0, (6) = 0}), Vector(6, {(1) = 2/5, (2) = 0, (3) = 1, (4) = 0, (5) = 0, (6) = 0}), Vector(6, {(1) = -3/5, (2) = 1, (3) = 0, (4) = 0, (5) = 0, (6) = 0})]

 

``

The dimension of M[1] is 3, and of M[2], 5. However, the basis vectors Maple has chosen for M[2] do not include the exact basis vectors chosen for M[1].

 

We now come to the crucial step, finding b[6], a vector in M[3] that is not in M[2] (and consequently, not in M[1] either). The examples in [1, 2] are simple enough that the authors can "guess" at the vector to be taken as b[6]. What we will do is take an arbitrary vector in M[3] and project it onto the 5-dimensional subspace M[2], and take the orthogonal complement as b[6].

``

A general vector in M[3] is

 

Z := `<,>`(u || (1 .. 6))

Vector[column]([[u1], [u2], [u3], [u4], [u5], [u6]])

 

``

A matrix that projects onto M[2] is

 

P := ProjectionMatrix(NN)

Matrix([[42/67, -15/67, 10/67, -25/67, 0, 10/67], [-15/67, 58/67, 6/67, -15/67, 0, 6/67], [10/67, 6/67, 63/67, 10/67, 0, -4/67], [-25/67, -15/67, 10/67, 42/67, 0, 10/67], [0, 0, 0, 0, 1, 0], [10/67, 6/67, -4/67, 10/67, 0, 63/67]])

 

``

The orthogonal complement of the projection of Z onto M[2] is then -P*Z+Z. This vector can be simplified by choosing the parameters in Z appropriately. The result is taken as b[6].

 

b[6] := 67*(eval(Z-Typesetting:-delayDotProduct(P, Z), Equate(Z, UnitVector(1, 6))))*(1/5)

Vector[column]([[5], [3], [-2], [5], [0], [-2]])

NULL

 

``

The other two members of this chain are then

 

b[4] := Typesetting:-delayDotProduct(C, b[6])

Vector[column]([[-132], [-12], [-169], [40], [92], [-79]])

b[1] := Typesetting:-delayDotProduct(C, b[4])

Vector[column]([[-67], [134], [67], [67], [0], [134]])

 

``

A general vector in M[2] is a linear combination of the five vectors that span the null space of C^2, namely, the vectors in the list NN. We obtain this vector as

 

ZZ := add(u || k*NN[k], k = 1 .. 5)

Vector[column]([[(2/5)*u1-u3+(2/5)*u4-(3/5)*u5], [u5], [u4], [u3], [u2], [u1]])

 

``

A vector in M[2] that is not in M[1] is the orthogonal complement of the projection of ZZ onto the space spanned by the eigenvectors spanning M[1] and the vector b[4]. This projection matrix is

 

PP := LinearAlgebra:-ProjectionMatrix(convert(`union`(LinearAlgebra:-NullSpace(C), {b[4]}), list))

Matrix([[69/112, -33/112, 19/112, -17/56, 0, 19/112], [-33/112, 45/112, 25/112, 13/56, 0, 25/112], [19/112, 25/112, 101/112, 1/56, 0, -11/112], [-17/56, 13/56, 1/56, 5/28, 0, 1/56], [0, 0, 0, 0, 1, 0], [19/112, 25/112, -11/112, 1/56, 0, 101/112]])

 

``

The orthogonal complement of ZZ, taken as b[5], is then

 

b[5] := 560*(eval(ZZ-Typesetting:-delayDotProduct(PP, ZZ), Equate(`<,>`(u || (1 .. 5)), LinearAlgebra:-UnitVector(4, 5))))

Vector[column]([[-9], [-59], [17], [58], [0], [17]])

 

``

Replace the vector a[2] with b[2], obtained as

 

b[2] := Typesetting:-delayDotProduct(C, b[5])

Vector[column]([[251], [-166], [197], [-139], [-112], [-166]])

 

 

The columns of the transition matrix Q can be taken as the vectors b[1], b[4], b[6], b[2], b[5], and the eigenvector a[3]. Hence, Q is the matrix

 

Q := `<|>`(b[1], b[4], b[6], b[2], b[5], N[3])

Matrix([[-67, -132, 5, 251, -9, -2], [134, -12, 3, -166, -59, 1], [67, -169, -2, 197, 17, -1], [67, 40, 5, -139, 58, 1], [0, 92, 0, -112, 0, 0], [134, -79, -2, -166, 17, 0]])

 

``

Proof that this matrix Q indeed sends A to its Jordan form consists in the calculation

 

1/Q.A.Q = Matrix([[2, 1, 0, 0, 0, 0], [0, 2, 1, 0, 0, 0], [0, 0, 2, 0, 0, 0], [0, 0, 0, 2, 1, 0], [0, 0, 0, 0, 2, 0], [0, 0, 0, 0, 0, 2]])``

 

NULL

The bases for M[k], k = 1, 2, 3, are not unique. The columns of the matrix Q provide one set of basis vectors, but the columns of the transition matrix generated by Maple, shown below, provide another.

 

JordanForm(A, output = 'Q')

Matrix([[-5, -43/5, -9/5, 7/5, -14/5, -3/5], [10, -4/5, -6/25, 1/5, -6/25, -3/25], [5, -52/5, -78/25, 13/5, -78/25, -39/25], [5, 13/5, 38/25, -2/5, 38/25, 4/25], [0, 6, 42/25, -1, 42/25, 21/25], [10, -29/5, -11/25, 1/5, -11/25, 7/25]])

 

``

I've therefore added to my to-do list the investigation into Maple's algorithm for determining an appropriate set of basis vectors that will support the Jordan form of a matrix.

 

References

 

NULL

[1] Linear Algebra and Matrix Theory, Evar Nering, John Wiley and Sons, Inc., 1963

[2] Matrix Methods: An Introduction, Richard Bronson, Academic Press, 1969

 

NULL

``

Download JordanForm_blog.mw

A wealth of knowledge is on display in MaplePrimes as our contributors share their expertise and step up to answer others’ queries. This post picks out one such response and further elucidates the answers to the posted question. I hope these explanations appeal to those of our readers who might not be familiar with the techniques embedded in the original responses.

Before I begin, a quick note that the content below was primarily created by one of our summer interns, Pia, with guidance and advice from me.

The Question: Rearranging the expression of equations

SY G wanted to be able to re-write an equation in terms of different variables.  SY G presented this example: 

I have the following two equations:

x1 = a-y1-d*y2;
x2 = a-y2-d*y1;

I wish to express the first equation in terms of y1 and x2, so that

x1 = c - b*y1+d*x2;

where c=a-a*d and b=1-d^2. How can I get Maple to rearrange the original equation x1 in term of y1, x2, c and b?

This question was answered by nm who provided code with a systematic approach:

restart;
eq1:=x1=a-y1-d*y2:
eq2:=x2=a-y2-d*y1:
z:=expand(subs(y2=solve(eq2,y2),eq1)):
z:=algsubs((a-a*d)=c,z):
algsubs((1-d^2)=b,z);

On the other hand, Carl Love answered this enquiry using a more direct and simple code:

simplify(x1=a-y1-d*y2, {a-y2-d*y1= x2, 1-d^2= b, a-a*d= c});

Let’s talk more about the expand, algsubs, subs, and simplify commands

First let’s take a look at the method nm used to solve the problem using the commands expand, subs, solve and algsubs.

The expand command, expand(expr, expr1, expr2, ..., exprn), distributes products over sums. This is done for all polynomials. For quotients of polynomials, only sums in the numerator are expanded; products and powers are left alone.

The solve command, solve(equations, variables), solves one or more equations or inequalities for their unknowns.

The subs command, subs(x=a,expr), substitutes a for x in the expression expr.

The function algsubs, algsubs(a = b, f),performs an algebraic substitution, replacing occurrences of a with b in the expression f.  It is a generalization of the subs command, which only handles syntactic substitution.

Let’s tackle the Maple code written by nm step by step:

1) restart;
The restart command is used to clear Maple’s internal memory

2)  eq1:=x1=a-y1-d*y2:
      eq2:=x2=a-y2-d*y1:
The names eq1 and eq2 were assigned to the equations SY G provided.

3) z:=expand(subs(y2=solve(eq2,y2),eq1)):
A new variable, z, was created, which will end up being x1 written in the terms SY G wanted.

  • solve(eq2,y2)
    • the solve command was used to solve the expression eq2 for the variable y2.

  • subs(y2=solve(eq2,y2),eq1)
    • The subs command was used to replace in expression eq1, y2 as determined by the solve step. 

  • expand(subs(y2=solve(eq2,y2),eq1))
    • The expand command was used to distribute products over sums. Note: this step served to ensure that the final output looked exactly how SY G wanted.

4) z:=algsubs((a-a*d)=c,z):
First, nm equated a-a*d to c, so later the algsubs command could be applied to substitute the new variable c into the expression z.

5) algsubs((1-d^2)=b,z);
Again, nm equated 1-d^2 to b, so later the algsubs command could be applied to substitute the new variable b into the expression z.

An alternate approach

Now let us check out Carl Love’s approach. Carl Love uses the simplify command in conjunction with side relations.

The simplify command has many calling sequences and one of them is the simplify(expr,eqns), that is known as simplify/siderels. A simplification of expr with respect to the side relations eqns is performed. The result is an expression which is mathematically equivalent toexpr but which is in normal form with respect to the specified side relations. Basically you are telling Maple to simplify the expression (expr) using the parameters (eqns) you gave to it.

 

I hope that you find this useful. If there is a particular question on MaplePrimes that you would like further explained, please let me know. 

A wealth of knowledge is on display in MaplePrimes as our contributors share their expertise and step up to answer others’ queries. This post picks out one such response and further elucidates the answers to the posted question. I hope these explanations appeal to those of our readers who might not be familiar with the techniques embedded in the original responses.

Before I begin, a quick note that the content below was primarily created by one of our summer interns, Pia, with guidance and advice from me.

The Question: Source Code of Math Apps

Eberch, a new Maple user, was interested in learning how to build his own Math Apps by looking at the source code of some of the already existing Math Apps that Maple offers.

Acer helpfully suggested that he look into the Startup Code of a Math App, in order to see definitions of procedures, modules, etc. He also recommended Eberch take a look at the “action code” that most of the Math Apps have which consist of function calls to procedures or modules defined in the Startup Code. The Startup Code can be accessed from the Edit menu. The function calls can be seen by right-clicking on the relevant component and selecting Edit Click Action.

Acer’s answer is correct and helpful. But for those just learning Maple, I wanted to provide some additional explanation.

Let’s talk more about building your own Math Apps

Building your own Math Apps can seem like something that involves complicated code and rare commands, but Daniel Skoog perfectly portrays an easy and straightforward method to do this in his latest webinar. He provides a clear definition of a Math App, a step-by-step approach to creating a Math App using the explore and quiz commands, and ways to share your applications with the Maple community. It is highly recommended that you watch the entire webinar if you would like to learn more about the core concepts of working with Maple, but you can find the Math App information starting at the 33:00 mark.

I hope that you find this useful. If there is a particular question on MaplePrimes that you would like further explained, please let me know. 

A wealth of knowledge is on display in MaplePrimes as our contributors share their expertise and step up to answer others’ queries. This post picks out one such response and further elucidates the answers to the posted question. I hope these explanations appeal to those of our readers who might not be familiar with the techniques embedded in the original responses.

Before I begin, a quick note that the content below was primarily created by one of our summer interns, Pia, with guidance and advice from me.

MaplePrimes member Thomas Dean wanted 1/2*x^(1/2) + 1/13*x^(1/3) + 1/26*x^(45/37)  to become  0.5*x^0.500000 + 0.07692307692*x^0.333333 + 0.03846153846*x^1.216216216  using the evalf command.

Here you can see the piece of code that Thomas Dean wrote in Maple:

eq:=1/2*x^(1/2) + 1/13*x^(1/3) + 1/26*x^(45/37);
evalf(eq);

Carl Love replied simply and effectively with a piece of code, using the evalindets command instead:

evalindets(eq, fraction, evalf);

As always, Love provided an accurate response, and it is absolutely correct. But for those just learning Maple, I wanted to provide some additional explanation.

The evalindets command, evalindets( expr, atype, transformer, rest ), is a particular combination of calls to eval and indets that allows you to efficiently transform all subexpressions of a given type by some algorithm. It encapsulates a common "pattern" used in expression manipulation and transformation.

Each subexpression of type atype is transformed by the supplied transformer procedure. Then, each subexpression is replaced in the original expression, using eval, with the corresponding transformed expression.

 

Note: the parameter restis an optional expression sequence of extra arguments to be passed to transformer. In this example it was not used.

I hope that you find this useful. If there is a particular question on MaplePrimes that you would like further explained, please let me know. 

1 2 3 4 5 6 7 Last Page 3 of 11