# 1. a)
> u:=<1,1,1,1,1>;
> v:=<1,-1,2,-2,3>;
> w:=<1,1,4,4,9>;

                                    [1]
                                    [ ]
                                    [1]
                                    [ ]
                               u := [1]
                                    [ ]
                                    [1]
                                    [ ]
                                    [1]


                                   [ 1]
                                   [  ]
                                   [-1]
                                   [  ]
                              v := [ 2]
                                   [  ]
                                   [-2]
                                   [  ]
                                   [ 3]


                                    [1]
                                    [ ]
                                    [1]
                                    [ ]
                               w := [4]
                                    [ ]
                                    [4]
                                    [ ]
                                    [9]

> with(LinearAlgebra);

  [Add, Adjoint, BackwardSubstitute, BandMatrix, Basis, BezoutMatrix,

        BidiagonalForm, BilinearForm, CharacteristicMatrix,

        CharacteristicPolynomial, Column, ColumnDimension,

        ColumnOperation, ColumnSpace, CompanionMatrix,

        ConditionNumber, ConstantMatrix, ConstantVector,

        CreatePermutation, CrossProduct, DeleteColumn, DeleteRow,

        Determinant, DiagonalMatrix, Dimension, Dimensions,

        DotProduct, Eigenvalues, Eigenvectors, Equal,

        ForwardSubstitute, FrobeniusForm, GaussianElimination,

        GenerateEquations, GenerateMatrix, GetResultDataType,

        GetResultShape, GivensRotationMatrix, GramSchmidt,

        HankelMatrix, HermiteForm, HermitianTranspose, HessenbergForm,

        HilbertMatrix, HouseholderMatrix, IdentityMatrix,

        IntersectionBasis, IsDefinite, IsOrthogonal, IsSimilar,

        IsUnitary, JordanBlockMatrix, JordanForm, LA_Main,

        LUDecomposition, LeastSquares, LinearSolve, Map, Map2,

        MatrixAdd, MatrixInverse, MatrixMatrixMultiply, MatrixNorm,

        MatrixScalarMultiply, MatrixVectorMultiply, MinimalPolynomial,

        Minor, Multiply, NoUserValue, Norm, Normalize, NullSpace,

        OuterProductMatrix, Permanent, Pivot, QRDecomposition,

        RandomMatrix, RandomVector, Rank, ReducedRowEchelonForm, Row,

        RowDimension, RowOperation, RowSpace, ScalarMatrix,

        ScalarMultiply, ScalarVector, SchurForm, SingularValues,

        SmithForm, SubMatrix, SubVector, SumBasis, SylvesterMatrix,

        ToeplitzMatrix, Trace, Transpose, TridiagonalForm, UnitVector,

        VandermondeMatrix, VectorAdd, VectorAngle,

        VectorMatrixMultiply, VectorNorm, VectorScalarMultiply,

        ZeroMatrix, ZeroVector, Zip]

> Transpose(u).v;

                                  3

> Transpose(u).w;

                                  19

> Transpose(v).w;

                                  27

# 1. b)
> Norm(u,2);

                               sqrt(5)

> Norm(v,2);

                               sqrt(19)

> Norm(u+v,2);

                               sqrt(30)

# 1. c)
# 
> evalf(arccos((Transpose(v).w)/(Norm(v,2)*Norm(w,2))));

                             .9549930364

> evalf(%*180/Pi);

                             54.71707045

# The first answer is in radians and the second is in degrees.
# 1. d)
> ((Transpose(u).v)/(Transpose(u).u))*u;

                                [3/5]
                                [   ]
                                [3/5]
                                [   ]
                                [3/5]
                                [   ]
                                [3/5]
                                [   ]
                                [3/5]

# 1. e)
> A := <u|v|w>;
> orthobasis:=NullSpace(Transpose(A));

                              [1     1    1]
                              [            ]
                              [1    -1    1]
                              [            ]
                         A := [1     2    4]
                              [            ]
                              [1    -2    4]
                              [            ]
                              [1     3    9]


                                    [-10]  [-6]
                                    [   ]  [  ]
                                    [  0]  [ 1]
                                    [   ]  [  ]
                     orthobasis := {[ 15], [ 8]}
                                    [   ]  [  ]
                                    [  1]  [ 0]
                                    [   ]  [  ]
                                    [ -6]  [-3]

# 2. a)
> b:=<-3,-2,6,-4,-8>;
# 

                                   [-3]
                                   [  ]
                                   [-2]
                                   [  ]
                              b := [ 6]
                                   [  ]
                                   [-4]
                                   [  ]
                                   [-8]

> xhat := MatrixInverse(Transpose(A).A).Transpose(A).b;

                                     [ 1]
                                     [  ]
                             xhat := [ 1]
                                     [  ]
                                     [-1]

> bhat := A.xhat;

                                     [ 1]
                                     [  ]
                                     [-1]
                                     [  ]
                             bhat := [-1]
                                     [  ]
                                     [-5]
                                     [  ]
                                     [-5]

# bhat = 1*u+1*v-1*w
# 2. b)
> Transpose(u).(b-bhat);

                                  0

> Transpose(v).(b-bhat);

                                  0

> Transpose(w).(b-bhat);

                                  0

> Ortho := <orthobasis[1]|orthobasis[2]>;

                                  [-10    -6]
                                  [         ]
                                  [  0     1]
                                  [         ]
                         Ortho := [ 15     8]
                                  [         ]
                                  [  1     0]
                                  [         ]
                                  [ -6    -3]

> LinearSolve(Ortho,b-bhat);

                                 [ 1]
                                 [  ]
                                 [-1]

# b-bhat= 1*orthobasis[1] - 1*orthobasis[2]
# 2. c)
> Norm(b-bhat,2);
> Norm(b-bhat,infinity);
> Norm(b-bhat,1);

                              2 sqrt(19)


                                  7


                                  16

# 3. a) The values found are, a=1, b=1, c=-1.
# 3. b)
> LS:=<u|v|<1,-1,8,-8,27>>;
> xhat2:=MatrixInverse(Transpose(LS).LS).Transpose(LS).b;

                              [1     1     1]
                              [             ]
                              [1    -1    -1]
                              [             ]
                        LS := [1     2     8]
                              [             ]
                              [1    -2    -8]
                              [             ]
                              [1     3    27]


                                    [-95 ]
                                    [--- ]
                                    [87  ]
                                    [    ]
                                    [1835]
                           xhat2 := [----]
                                    [522 ]
                                    [    ]
                                    [-311]
                                    [----]
                                    [522 ]

> bhat2 := LS.xhat2;

                                    [ 53 ]
                                    [ -- ]
                                    [ 29 ]
                                    [    ]
                                    [-349]
                                    [----]
                                    [ 87 ]
                                    [    ]
                                    [ 34 ]
                           bhat2 := [ -- ]
                                    [ 29 ]
                                    [    ]
                                    [-292]
                                    [----]
                                    [ 87 ]
                                    [    ]
                                    [-577]
                                    [----]
                                    [ 87 ]

> evalf(Norm(b-bhat2,2));

                             7.276161676

> evalf(Norm(b-bhat,2));

                             8.717797888

# So the equation for this least squares fit is
# y=-95/87+(1835/522)x-(311/522)x^3 and this fit is closer than the
# quadratic fit.
# 4. a)
> u:=<1,0,-1,0,3>;
> v:=<-5,6,5,6,7>;
> w:=<4,-7,-4,1,1>;
> z:=<-3,11,1,-5,5>;
> b:=<-5,19,3,-13,-1>;

                                   [ 1]
                                   [  ]
                                   [ 0]
                                   [  ]
                              u := [-1]
                                   [  ]
                                   [ 0]
                                   [  ]
                                   [ 3]


                                   [-5]
                                   [  ]
                                   [ 6]
                                   [  ]
                              v := [ 5]
                                   [  ]
                                   [ 6]
                                   [  ]
                                   [ 7]


                                   [ 4]
                                   [  ]
                                   [-7]
                                   [  ]
                              w := [-4]
                                   [  ]
                                   [ 1]
                                   [  ]
                                   [ 1]


                                   [-3]
                                   [  ]
                                   [11]
                                   [  ]
                              z := [ 1]
                                   [  ]
                                   [-5]
                                   [  ]
                                   [ 5]


                                   [ -5]
                                   [   ]
                                   [ 19]
                                   [   ]
                              b := [  3]
                                   [   ]
                                   [-13]
                                   [   ]
                                   [ -1]

> gs1 := u;

                                    [ 1]
                                    [  ]
                                    [ 0]
                                    [  ]
                             gs1 := [-1]
                                    [  ]
                                    [ 0]
                                    [  ]
                                    [ 3]

> r12 := (Transpose(gs1).v)/(Transpose(gs1).gs1);
> gs2 := v- r12*gs1;

                               r12 := 1


                                    [-6]
                                    [  ]
                                    [ 6]
                                    [  ]
                             gs2 := [ 6]
                                    [  ]
                                    [ 6]
                                    [  ]
                                    [ 4]

> r13 := (Transpose(gs1).w)/(Transpose(gs1).gs1);
> r23 := (Transpose(gs2).w)/(Transpose(gs2).gs2);
> gs3 := w- r13*gs1 - r23*gs2;

                               r13 := 1


                             r23 := -1/2


                                    [ 0]
                                    [  ]
                                    [-4]
                                    [  ]
                             gs3 := [ 0]
                                    [  ]
                                    [ 4]
                                    [  ]
                                    [ 0]

> r14 := (Transpose(gs1).z)/(Transpose(gs1).gs1);
> r24 := (Transpose(gs2).z)/(Transpose(gs2).gs2);
> r34 := (Transpose(gs3).z)/(Transpose(gs3).gs3);
> gs4 := z- r14*gs1 - r24*gs2 - r34*gs3;

                               r14 := 1


                              r24 := 1/2


                              r34 := -2


                                    [-1]
                                    [  ]
                                    [ 0]
                                    [  ]
                             gs4 := [-1]
                                    [  ]
                                    [ 0]
                                    [  ]
                                    [ 0]

# So gs1, gs2, gs3 and gs4 form a orthogonal basis for the span of
# u,v,w,z.
# 4. b)
> Q := <gs1|gs2|gs3|gs4>;
> R := <<1,0,0,0>|<r12,1,0,0>|<r13,r23,1,0>|<r14,r24,r34,1>>;

                          [ 1    -6     0    -1]
                          [                    ]
                          [ 0     6    -4     0]
                          [                    ]
                     Q := [-1     6     0    -1]
                          [                    ]
                          [ 0     6     4     0]
                          [                    ]
                          [ 3     4     0     0]


                          [1    1     1       1 ]
                          [                     ]
                          [0    1    -1/2    1/2]
                     R := [                     ]
                          [0    0     1      -2 ]
                          [                     ]
                          [0    0     0       1 ]

> A := <u|v|w|z>;

                          [ 1    -5     4    -3]
                          [                    ]
                          [ 0     6    -7    11]
                          [                    ]
                     A := [-1     5    -4     1]
                          [                    ]
                          [ 0     6     1    -5]
                          [                    ]
                          [ 3     7     1     5]

> Q.R;

                        [ 1    -5     4    -3]
                        [                    ]
                        [ 0     6    -7    11]
                        [                    ]
                        [-1     5    -4     1]
                        [                    ]
                        [ 0     6     1    -5]
                        [                    ]
                        [ 3     7     1     5]

# 4. c)
> x := MatrixInverse(R).MatrixInverse(Transpose(Q).Q).Transpose(Q).b;

                                   [ 1]
                                   [  ]
                                   [-1]
                              x := [  ]
                                   [-2]
                                   [  ]
                                   [ 1]

>