Stat/Math
Software Support
Software Consulting
Software Availability
Software Price
Contact

User Support
Documentation
Knowledge Base
Education
Consulting
Podcasts

Systems & Services
Cyberinfrastructure
Supercomputers
Grid Computing
Storage
Visualization
Digital Libraries & Data

Results & Impact
Publications
Grants & Grant Info
Events & Outreach
Economic Impact
Survey Results

Vision & Planning
News & Features

# Linear Algebra with Maple

The output here is in character mode, rather than the default "high-resolution" mode. It's still formatted in two dimensions, so you need a monospaced font (for preformatted material) and a "sufficiently" wide window.

> with(linalg);
```BlockDiagonal, GramSchmidt, JordanBlock, Wronskian, add,
basis, bezout, blockmatrix, charmat, charpoly, col, coldim
colspace, colspan, companion, concat, cond, copyinto,
crossprod, curl, definite, delcols, delrows, det, diag,
diverge, dotprod, eigenvals, eigenvects, entermatrix, equal,
exponential, extend, ffgausselim, fibonacci, frobenius,
hessian, hilbert, htranspose, ihermite, indexfunc,
innerprod, intbasis, inverse, ismith, iszero, jacobian,
jordan, kernel, laplacian, leastsqrs, linsolve, matrix,
minor, minpoly, mulcol, mulrow, multiply, norm, normalize,
nullspace, orthog, permanent, pivot, potential, randmatrix,
randvector, rank, ratform, row, rowdim, rowspace, rowspan,
rref, scalarmul, singularvals, smith, stack, submatrix,
subvector, sumbasis, swapcol, swaprow, sylvester, toeplitz,
trace, transpose, vandermonde, vecpotent, vectdim, vector
```

This is a list of the functions added by loading the "linalg" package. To find out more about vector, for example, type ?vector

> v := vector(3, [a, b, c]); m := matrix(2,2, [[a, b], [c, d]]);
```
v := [ a, b, c ]

[ a  b ]
m := [      ]
[ c  d ]
```
> evalm( m &* [2,3] );
```
[2 a + 3 b, 2 c + 3 d]
```
> linsolve(m,[2,3]);
```
- 2 d + 3 b    - 2 c + 3 a
[ -----------, - ----------- ]
b c - a d      b c - a d
```
> det(m);
```                            a d - b c
```
> eigenvals(m);
```
2            2         1/2
1/2 a + 1/2 d + 1/2 (a  - 2 a d + d  + 4 b c)   ,

2            2         1/2
1/2 a + 1/2 d - 1/2 (a  - 2 a d + d  + 4 b c)
```
```                                 - 1/2 a + 1/2 d - 1/2 %1
[1/2 a + 1/2 d + 1/2 %1, 1, {[ - ------------------------, 1 ]}]
c

,

- 1/2 a + 1/2 d + 1/2 %1
[1/2 a + 1/2 d - 1/2 %1, 1, {[ - ------------------------, 1 ]}]
c

2            2         1/2
%1 :=            (a  - 2 a d + d  + 4 b c)

```

A common error is trying to compute the determinant, or eigenvalues, of a large symbolic n by n matrix. The problem is called combinatorial explosion -- there might be n2 expressions in the answer which are each n! terms long!

The Hilbert matrix is a favorite example of numeric (as opposed to symbolic) linear algebra.

> H5 := hilbert(5);
```
[  1   1/2  1/3  1/4  1/5 ]
[                         ]
[ 1/2  1/3  1/4  1/5  1/6 ]
[                         ]
H5 := [ 1/3  1/4  1/5  1/6  1/7 ]
[                         ]
[ 1/4  1/5  1/6  1/7  1/8 ]
[                         ]
[ 1/5  1/6  1/7  1/8  1/9 ]
```

The Hilbert matrix is badly conditioned, which means that it corresponds to a system of linear equations which are hard to solve accurately.

> cond(H5);
```
943656
```

This means that errors [for example in the measured values of the right-hand side of the equations] may be multiplied by a factor of nearly a million!

A curious feature of the Hilbert matrix is that the entries of its inverse are all integers. Maple computes them exactly:

> inverse(H5);
```
[   25    -300     1050    -1400     630  ]
[                                         ]
[  -300   4800    -18900   26880   -12600 ]
[                                         ]
[  1050  -18900   79380   -117600   56700 ]
[                                         ]
[ -1400   26880  -117600   179200  -88200 ]
[                                         ]
[  630   -12600   56700    -88200   44100 ]
```

The method for finding eigenvectors of m that Maple used above would not be appropriate for this matrix; if it were symbolic, this matrix would be too big; but more importantly, since it's numeric, there is another method [the QR method] which is less sensative to the numeric instability implied by a big condition number, and is also computed more quickly. Maple uses that method, if the inputs are floating-point [decimal] numbers.

> eigsys := eigenvects( map( evalf, H5) );
```
eigsys :=

[.01140749151, 1,
{ [ .2142136245, -.7241021311, -.1204532771, .3095739659, .5651934124 ] }],

[1.567050691, 1,
{ [ .7678547356, .4457910609, .3215782948, .2534389442, .2098226368 ] }],

[.2085342186, 1,
{ [-.6018714786, .2759134174, .4248766220, .4439030401, .4290133540 ] }],
-5
[.3288052655*10  , 1,
{ [ .006173874786, -.1166928633, .5061638445, -.7671911293, .3762453886 ] }],

[.0003058977746, 1,
{ [ .04716180523, -.4326673028, .6673503032, .2330247284, -.5576000510 ] }]
```

"eigsys" is a list of compound objects: the eigenvalues, each with its multiplicity and corresponding eigenvector.

To select the second component of the second eigenvector, do this:

> eigsys[2][3][1][2];
```
.4457910609
```

To sort this out, examine in turn

• eigsys[2](the "system:" second eigenvector, its eigenvalue and multiplicity),
• eigsys[2][3](the set consisting of the second eigenvector system's vector),
• eigsys[2][3][1](the vector itself).

CAUTION: The order of the eigensystem cannot be predicted, from one session to the next.

The next command finds the singular values to 4 floating-point digits, among other things.

> svdlist := evalf(Svd(H5, 'U', 'V'), 4);
```
svdlist := [ 1.567, .2086, .01136, .0002933, .00001142
```
> d := diag( seq(svdlist[k], k=1..5));
```
[ 1.567    0       0        0         0     ]
[                                           ]
[   0    .2086     0        0         0     ]
[                                           ]
d := [   0      0    .01136      0         0     ]
[                                           ]
[   0      0       0    .0002933      0     ]
[                                           ]
[   0      0       0        0     .00001142 ]
```

This next command checks the accuracy of the result just computed, carrying 20 digits:

> evalf(evalm( (transpose(U) &* H5 &* V) - d), 20);
```
[                                                                        ]
[       -18            -19           -19            -19           -19    ]
[ -.1*10         -.5*10         .4*10         -.2*10        -.3*10       ]
[                                                                        ]
[       -20            -19           -20            -20           -20    ]
[ -.9*10         -.1*10         .2*10          .4*10         .5*10       ]
[                                                                        ]
[         -19           -20                         -21            -20   ]
[ -.108*10        .61*10        0              .9*10        -.49*10      ]
[                                                                        ]
[           -19           -20           -20           -20            -20 ]
[ -.27040*10     -.3301*10      .6147*10      -.973*10      -.5038*10    ]
[                                                                        ]
[            -20            -20           -20            -20         -21 ]
[ -.395257*10    -.619981*10    .750459*10    -.182306*10    .7701*10    ]
[                                                                        ]
```