# MAPLE ASIIGNMENT 5
# The first command loads the matrix algebra operations into Maple.
# Change the : to ; if you wish to see all the commands in this package.
> with(linalg):
Warning, new definition for norm
Warning, new definition for trace
# Next define a few matrices. You can change the numbers to any numbers
# you like.
> M1:=matrix(3,3,[1,0,-1,1,5,7,3,3,9]);

                               [1    0    -1]
                               [            ]
                         M1 := [1    5     7]
                               [            ]
                               [3    3     9]

> M2:=matrix(3,3,[1,0,-1,3,6,21,3,3,9]);

                               [1    0    -1]
                               [            ]
                         M2 := [3    6    21]
                               [            ]
                               [3    3     9]

# The next commands create some other matrices from M1 and M2.
> M3:=submatrix(M2,2..3,1..2);

                                  [3    6]
                            M3 := [      ]
                                  [3    3]

> M4:=submatrix(M1,2..3,1..3);

                               [1    5    7]
                         M4 := [           ]
                               [3    3    9]

> M5:=stack(M2,M4);

                               [1    0    -1]
                               [            ]
                               [3    6    21]
                               [            ]
                         M5 := [3    3     9]
                               [            ]
                               [1    5     7]
                               [            ]
                               [3    3     9]

> M6:=concat(M3,M4);

                          [3    6    1    5    7]
                    M6 := [                     ]
                          [3    3    3    3    9]

#  The next command illustrates multiplication by a scalar.  Evalm is
# needed for the answer to appear 
# on the screen. Refer back to M6 to see what has happened
> evalm(3*M6);

                      [9    18    3    15    21]
                      [                        ]
                      [9     9    9     9    27]

# 
# The aritmetic operations for matrices are + , - , &* , ^ . For easier
# reference the next command 
# brings M1 and M2 back on to the screen side by side in a 3X6 matrix.
# Execute the next 
# command and after you are sure you understand what has happened
# execute the command agaiin 
# with + changed first to - and then to &* .
> concat(M1,M2);

                    [1    0    -1    1    0    -1]
                    [                            ]
                    [1    5     7    3    6    21]
                    [                            ]
                    [3    3     9    3    3     9]

> evalm(M1+M2);

                          [-2    -3    -10]
                          [               ]
                          [37    51    167]
                          [               ]
                          [39    45    141]

# If you wish to check your understanding of matrix multiplication the
# following command calculates the 
# [2,2] entry of the product matrix. Change the numbers in this command
# to calculate some other 
# entries in the product matrix. If one of the numbers is negative put
# it in brackets.
> (1*0)+(5*6)+(7*3);

                                  51

# The next comand does a power of a matrix. You can check the answers
# out in the same way the 
# product was checked.
> evalm(M1^2);

                          [-2    -3    -10]
                          [               ]
                          [27    46     97]
                          [               ]
                          [33    42     99]

# The following illustrate some of the properties of matrix algebra. A
# matrix is transposed when its 
# rows and columns are switched..  The properties are the associative ,
# commutative and distributive 
# properties.
> M7:=transpose(M4)&*M4;

                               [1    3]
                               [      ]
                         M7 := [5    3] &* M4
                               [      ]
                               [7    9]

> evalm(M4);

                            [1    5    7]
                            [           ]
                            [3    3    9]

> evalm(M7);

                          [10    14     34]
                          [               ]
                          [14    34     62]
                          [               ]
                          [34    62    130]

> evalm((M1+M2)+M7);

                          [12    14     32]
                          [               ]
                          [18    45     90]
                          [               ]
                          [40    68    148]

> evalm(M1+(M2+M7));

                          [12    14     32]
                          [               ]
                          [18    45     90]
                          [               ]
                          [40    68    148]

> evalm((M1&*M2)&*M7);

                       [-402     -750    -1554]
                       [                      ]
                       [6762    12606    26130]
                       [                      ]
                       [5814    10818    22446]

> evalm(M1&*(M2&*M7));

                       [-402     -750    -1554]
                       [                      ]
                       [6762    12606    26130]
                       [                      ]
                       [5814    10818    22446]

> evalm((M1+M2)-(M2+M1));

                                  0

# The commutative law does not however work with multiplication.
> evalm((M1&*M2)-(M2&*M1));

                         [  0      0      0]
                         [                 ]
                         [-35    -42    -61]
                         [                 ]
                         [  6      3     42]

# The final property illustrated here is one of the distributive laws.
> evalm((M1+M2)&*M7);

                        [ -48     -96    -192]
                        [                    ]
                        [1146    2166    4458]
                        [                    ]
                        [ 756    1404    2916]

> evalm((M1&*M7)+(M2&*M7));

                        [ -48     -96    -192]
                        [                    ]
                        [1146    2166    4458]
                        [                    ]
                        [ 756    1404    2916]

# Diagonal matrices are created with the diag command.  This can be used
# to create an identity matrix.
> diag(-3,2,Pi); 

                           [-3    0    0 ]
                           [             ]
                           [ 0    2    0 ]
                           [             ]
                           [ 0    0    Pi]

> iden:=diag(1,1,1);

                                [1    0    0]
                                [           ]
                        iden := [0    1    0]
                                [           ]
                                [0    0    1]

# Inverse matrices are either A^(-1)  or inverse(A).  The following
# calculate the inverse, check that it 
# works and illustrate the properties of inverses.  You should also try
# to calculate the inverse of M2 .
> M8:=M1^(-1);

                                     1
                              M8 := ----
                                     M1

> evalm(");

                        [        -1         ]
                        [2/3     --     5/36]
                        [        12         ]
                        [                   ]
                        [1/3     1/3    -2/9]
                        [                   ]
                        [        -1         ]
                        [-1/3    --     5/36]
                        [        12         ]

> evalm(M1&*M8-iden);

                            [0    0    0]
                            [           ]
                            [0    0    0]
                            [           ]
                            [0    0    0]

> evalm(M8&*M1-iden);

                            [0    0    0]
                            [           ]
                            [0    0    0]
                            [           ]
                            [0    0    0]

> evalm(M8^(-1));

                            [1    0    -1]
                            [            ]
                            [1    5     7]
                            [            ]
                            [3    3     9]

> M9:=evalm(M1+M7);

                             [11    14     33]
                             [               ]
                       M9 := [15    39     69]
                             [               ]
                             [37    65    139]

# The next commands illustrate the correct and the incorrect way to get
# the inverse of the product of twp matrices.
> evalm((M1&*M9)^(-1));

                       [46      181     1483 ]
                       [--     -----    -----]
                       [81     16848    50544]
                       [                     ]
                       [29     257      -361 ]
                       [--     ----     -----]
                       [81     4212     12636]
                       [                     ]
                       [-26    -539      331 ]
                       [---    -----    -----]
                       [81     16848    50544]

> evalm((M9^(-1))&*(M1^(-1)));

                       [46      181     1483 ]
                       [--     -----    -----]
                       [81     16848    50544]
                       [                     ]
                       [29     257      -361 ]
                       [--     ----     -----]
                       [81     4212     12636]
                       [                     ]
                       [-26    -539      331 ]
                       [---    -----    -----]
                       [81     16848    50544]

> evalm((M1^(-1))*(M9^(-1)));

                        [10    2867     -1939]
                        [--    -----    -----]
                        [27    50544    16848]
                        [                    ]
                        [11    1915     -731 ]
                        [--    -----    ---- ]
                        [27    12636    4212 ]
                        [                    ]
                        [-8    -4297    1913 ]
                        [--    -----    -----]
                        [27    50544    16848]

> evalm((M1^(-1))+(M9^(-1)));

                        [        41      -7 ]
                        [4/3     ---     -- ]
                        [        702     78 ]
                        [                   ]
                        [        194     -16]
                        [2/3     ---     ---]
                        [        351     39 ]
                        [                   ]
                        [        -157    23 ]
                        [-2/3    ----    -- ]
                        [        702     78 ]

> evalm((M1+M9)^(-1));

                         [        13     -43]
                         [1/3     ---    ---]
                         [        504    504]
                         [                  ]
                         [        31     -25]
                         [1/6     ---    ---]
                         [        252    252]
                         [                  ]
                         [        -4     19 ]
                         [-1/6    --     ---]
                         [        63     252]

# The last two commands have illustrated a dumb artihmetical mistake you
# should not make even with numbers.
