Exercise 3.4 Solutions

2.

A=(({{1, 1, 6, 7}, {-1, 2, 9, 2}, {1, -2, 3, 10}}))

A1 = rowadd[1, 1, 2, A]

A2 = rowadd[-1, 1, 3, A1]

A3 = rowadd[1, 2, 3, A2]

6.

Clear[A, B]

quadr[x_] := A + B x  +  C (x^2)

Solve[{quadr[1] == 4, quadr[2] == 7, 
        quadr[3] == 14}, {A, B, C}]

quadr[1]

quadr[2]

quadr[3]

am=(({{1, 1, 1}, {1, 2, 4}, {1, 3, 9}}))

am1 = rowadd[-1, 1, 2, am]

am2 = rowadd[-1, 1, 3, am1]

am3 = rowadd[-2, 2, 3, am2]

bmat=(({{4}, {7}, {14}}))

BackSubDetails[ am3, bmat, 3 ]

7.

Clear[A, B]

cubic[x_] := A + B x  +  C (x^2) + D (x^3)

cubic[0]

cubic[1]

cubic[2]

cubic[3]

am=(({{1, 0, 0, 0}, {1, 1, 1, 1}, {1, 2, 4, 8}, {2, 3, 9, 27}}))

am1 = rowadd[-1, 1, 2, am]

am2 = rowadd[-1, 1, 3, am1]

am3 = rowadd[-2, 1, 4, am2]

am4 = rowadd[-2, 2, 3, am3]

am5 = rowadd[-3, 2, 4, am4]

am6 = rowadd[-3, 3, 4, am5]

bmat=(({{0}, {1}, {2}, {2}}))

BackSubDetails[ am6, bmat, 4 ]

12.

augm=(({{1, 1, 0, 0, 7}, {2, 3, -1, 0, 9}, {0, 4, 2, 3, 10}, {0, 0, 2, 6, 2}}))

am1 = rowadd[-2, 1, 2, augm]

am2 = rowadd[-4, 2, 3, am1]

am3 = rowadd[-1/3, 3, 4, am2]

BackSubDetails[(({{1, 1, 0, 0}, {0, 1, -1, 0}, {0, 0, 6, 3}, {0, 0, 0, 5}})),(({{7}, {-5}, {30}, {-8}})),4]

15.

hmatrix = HilbertMatrix[4]

GaussJordan[HilbertMatrix[4], 4]

GaussJordan[N[HilbertMatrix[4], 5], 4]

Algorithms and Programs 1.

See Tri-diagonal program below

Solution 2 (a).

n = 6;   

    d = Table[2.0, {n}]

a = c=Table[1.0, {n - 1}]

b=Table[
        16.0 - Abs[i - 16.0], {i, 1, n}]

Z2=TriDiagonal[a, d, c, b]

We can compare this answer with the previous one in Example 1.

Mathematica is able to get a slightly better answer.

Solution 1 (a).

n = 6;  
B = 
        Table[16.0 - Abs[i - 16.0], {i, 1, 
            n}];   
Print["B = ", B ];

The column form for B is:

Print["B = ", MatrixForm[B] ];

The matrix A can be constructed with the command.

\(\(\(\(A = DiagonalMatrix[Table[2.0, {n}]];\)\ \ 

    \(For[i = 1, 
        i ≤ 
          n - 1, \(i++\), 
\(A\_\(\([[\)\(i, i + 1\)\(]]\)\) = \(A\_\(\([[\)\(i + 1, i\)\(]]\)\) = 1.0\);\)];\)\ \ 

    \(Print["A = ", MatrixForm[A]\ ];\)\)\(\ \ \)\)\)

We can print out the system to be solved if we wish.

Clear[x];

    vars  =  Table[(Subscript["x", i]), {i, 1, n}];

    Print["Solve the system"];

     Print[MatrixForm[A], MatrixForm[vars], " = ", 
        MatrixForm[B] ];

Solve the system Mathematica's built in  LinearSolve[A,B]  procedure.

Z1 = LinearSolve[A, B];  
\(\(Print[\*""\!\(X\_1\) = "", MatrixForm[vars], " = ", 
        MatrixForm[Z1]\ , " = ", 
        NumberForm[MatrixForm[Z1], 17]];\)\)

Notice that this linear system is producing some round off error near the 14th decimal place.

Verify the answer.

Print["Does A X = B ?"];

     Print[MatrixForm[A], MatrixForm[Chop[Z1]], " = ", 
        MatrixForm[B] ];    

     Print[MatrixForm[A . Z1], " = ", 
        MatrixForm[B] ];    

     Print[
        MatrixForm[A . Z1]  ==  
          MatrixForm[B] ];

Algorithms and Programs 2.

poly6th[x_] := 
      a1 + a2 x  + a3 (x^2) + a4 (x^3) + a5 (x^4) + a6 (x^5) + 
        a7 (x^6)

TableForm[Table[poly6th[k], {k, 0, 6}]]

mata = {{1, 0, 0, 0, 0, 0, 0, 1}, 
        Flatten[Append[Table[1^k, {k, 0, 6}], 3]], 
        Flatten[Append[Table[2^k, {k, 0, 6}], 2]], 
        Flatten[Append[Table[3^k, {k, 0, 6}], 1]], 
        Flatten[Append[Table[4^k, {k, 0, 6}], 3]], 
        Flatten[Append[Table[5^k, {k, 0, 6}], 2]], 
        Flatten[Append[Table[6^k, {k, 0, 6}], 1]]}

LandUmod[mata, 7]

N[Solve[
        mata . {x1, x2, x3, x4, x5, x6, x7} == {1, 3, 2, 1, 3, 2, 
            1}, {x1, x2, x3, x4, x5, x6, x7}]]

g1 = 
        Plot[1 - 1.8  x + 11.025  x^2 - 10.5625  x^3 + 3.9375  x^4 - 
            0.6375  x^5 + 0.0375  x^6, {x, 0, 8}, 
          DisplayFunction -> Identity];

g2 = 
        ListPlot[{{0, 1}, {1, 3}, {2, 2}, {3, 1}, {4, 3}, {5, 2}, {6, 1}}, 
          PlotStyle -> PointSize[0.02], 
          DisplayFunction -> Identity];

Show[g1, g2, DisplayFunction -> $DisplayFunction];

Algorithms and Programs 3.

A1 = Table[i^(j - 1), {i, 1, 3}, {j, 1, 3}]

B1 = 
      Transpose[{Flatten[
            Append[{3}, 
              Table[((i^3) - 1)/(i - 1), {i, 2, 
                  3}]]]}]

AppendRows[A1, B1]

LandUmod[AppendRows[A1, B1], 3]

solver[
        n_] := 
Module[{amat, 
          bmat}, 
amat = 
          Table[i^(j - 1), {i, 1, n}, {j, 1, 
              n}]; 
bmat = 
          Transpose[{Flatten[
                Append[{n}, 
                  Table[((i^n) - 1)/(i - 1), {i, 2, 
                      n}]]]}]; 
Print[
          MatrixForm[AppendRows[amat, bmat]]]; 
LandUmod[
          AppendRows[amat, bmat], n]
]

solver[3]

Functions Used in this notebook

<< LinearAlgebra`MatrixManipulation`

$Post=If[MatrixQ[#],MatrixForm[#] ,#]&;

rowadd[scalar_,r1_,r2_,mat_]:=
	Module[{result,mat1=mat},
  	mat1[[r2]]=(scalar mat1[[r1]])+mat1[[r2]];
	mat1
	]

rowmultiply[scalar_,r1_,mat_]:=
	Module[{mat1=mat},
	mat1[[r1]]=scalar mat1[[r1]];
	mat1
	]

rowexchange[r1_,r2_,mat_]:=
	Module[{mat1=mat,temp},
	temp=mat1[[r1]];
	mat1[[r1]]=mat[[r2]];
	mat1[[r2]]=temp;
	mat1
	]

GaussJordan[A0_, 
        n_]:=
Module[ {A = A0, i, 
          p},  
Print[ 
          MatrixForm[A] ];  
For[ p = 1,  
          p ≤ 
            n,  p++,  
(Subscript[A, [[p]]])  =  ((Subscript[A, [[p]]])/(Subscript[A, [[p, p]]]));  
For[ i = 1,  
            i ≤ 
              n,  i++,  
If[ 
                i ≠ 
                  p,  
(Subscript[A, [[i]]])  =  
                  (Subscript[A, [[i]]])  -  (Subscript[A, [[i, p]]]) 
                      (Subscript[A, [[p]]])  ]; ];  
Print[
            MatrixForm[A]];  ];  ]

\(\(\(TriDiagonal[a0_, d0_, c0_, 
          b0_]\  := 
Module[{\ a = a0, b = b0, c = c0, 
            d = d0, k, m, n = Length[b0], x, 
            showmat\ }, 
showmat = 
            DiagonalMatrix[d]; 
For[i = 1, 
            i ≤ n - 1, \(i++\), 
            showmat[\([i, i + 1]\)] = c[\([i]\)]; 
            showmat[\([i + 1, i]\)] = a[\([i]\)];]; 
          Print[MatrixForm[showmat]]; 
Print[
            MatrixForm[b]]; 
For[\ k = 2, \ 
            k ≤ 
              n, \ \(k++\), 
d\_\(\([[\)\(k\)\(]]\)\)\  = \ d\_\(\([[\)\(k\)\(]]\)\)\ \  - \ \(\(\ \)\(a\_\(\([[\)\(k - 1\)\(]]\)\)\)\)\/\(\(\ \)\(d\_\(\([[\)\(k - 1\)\(]]\)\)\)\)\ \ c\_\(\([[\)\(k - 1\)\(]]\)\); \ 
            Print["k=", 
              k, \*"" (Subscript[d, [[k]]])="", 
              d\_\(\([[\)\(k\)\(]]\)\)\ ]; 
b\_\(\([[\)\(k\)\(]]\)\)\  = \ 
              b\_\(\([[\)\(k\)\(]]\)\)\  - \ \(\(\ \)\(a\_\(\([[\)\(k - 1\)\(]]\)\)\)\)\/\(\(\ \)\(d\_\(\([[\)\(k - 1\)\(]]\)\)\)\)\ b\_\(\([[\)\(k - 1\)\(]]\)\); \ 
            Print[\*""(Subscript[b, [[k]]]) ="", 
              b\_\(\([[\)\(k\)\(]]\)\)\ ]; 
showmat[\([k, k]\)] = d[\([k]\)]; 
            showmat[\([k, k - 1]\)] = 0; Print[MatrixForm[showmat]]; 
            Print[MatrixForm[b]]]; \ 
x = 
            Table[0, {n}]; \ 
x\_\(\([[\)\(n\)\(]]\)\)\  = \(\(\ \)\(b\_\(\([[\)\(n\)\(]]\)\)\)\)\/d\_\(\([[\)\(n\)\(]]\)\); \ 
For[\ k = n - 1, \ 
            1 ≤ 
              k, \ \(k--\), 
\(x\_\(\([[\)\(k\)\(]]\)\)\  = \ \(b\_\(\([[\)\(k\)\(]]\)\)\  - \ c\_\(\([[\)\(k\)\(]]\)\)\ x\_\(\([[\)\(k + 1\)\(]]\)\)\)\/d\_\(\([[\)\(k\)\(]]\)\);\)\ ]; \ 
Return[x];\ ];\)\(\ \)\)\)

\(\(\(BackSubDetails[\ a0_, b0_, 
          n_\ ] := 
Module[\ {\ a = a0, b = b0, i, j, x, 
            X\ }, 
x\  = \ 
            Table["x"\_i, {i, 1, n}]; \ 
X\  = \ 
            Table["x"\_i, {i, 1, n}]; \ 
For[\ i = n, \ 
            1 ≤ i, \ \(i--\), 
Print[
              Flatten[Chop[
                      a\_\(\([[\)\(i\)\(]]\)\)] . X] == 
                b\_\(\([[\)\(i, 1\)\(]]\)\)]; 
If[i < n, 
              Print[Flatten[
                      Chop[a\_\(\([[\)\(i\)\(]]\)\)]] . Flatten[x] == 
                  b\_\(\([[\)\(i, 1\)\(]]\)\)]; 
Print[\(a\_\(\([[\)\(i, i\)\(]]\)\)\) 
                    x\_\(\([[\)\(i\)\(]]\)\) == \((b\_\(\([[\)\(i\)\(]]\)\) - ∑\+\(j = \ i + 1\)\%n\( a\_\(\([[\)\(i, j\)\(]]\)\)\) x\_\(\([[\)\(j\)\(]]\)\))\)\_\(\([[\)\(1\)\(]]\)\)];]; 
x\_\(\([[\)\(i\)\(]]\)\)\  = \ \(b\_\(\([[\)\(i\)\(]]\)\) - ∑\+\(j = \ i + 1\)\%n\( a\_\(\([[\)\(i, j\)\(]]\)\)\) x\_\(\([[\)\(j\)\(]]\)\)\)\/a\_\(\([[\)\(i, i\)\(]]\)\); \ 
Print["x"\_i, " = ", 
              x\_\(\([[\)\(i, 1\)\(]]\)\), "
"]]; \ 
Return[
            x];\ ];\)\(\ \)\)\)

ForeSub[ a0_, b0_, 
          n_ ] := 
Module[ { a = a0,  b = b0,  i, j , 
            x}, 
x  =  
            Table[0, {n}];  
For[ i = 1,  
            i ≤ 
              n,  i++, 
(Subscript[x, [[i]]])  =  (((Subscript[b, [[i]]]) - (∑_(j =  1)^(i - 1))(Subscript[a, [[i, j]]]) (Subscript[x, [[j]]]))/(Subscript[a, [[i, i]]])) ];  
Return[
            x]; ];

BackSub[ a0_, b0_, 
          n_ ] := 
Module[ { a = a0,  b = b0,  i, j, 
            x }, 
x  =  
            Table[0, {n}];  
For[ i = n,  
            1 ≤ 
              i,  i--, 
(Subscript[x, [[i]]])  =  (((Subscript[b, [[i]]]) - (∑_(j =  i + 1)^n) (Subscript[a, [[i, j]]]) (Subscript[x, [[j]]]))/(Subscript[a, [[i, i]]])); ];  
Return[
            x]; ];

LandUmod[A0_, 
        n_] := 
Module[ {A = A0, B = B0, i, 
          p},  
U = A0;  
L = 
          IdentityMatrix[n];  
For[ p = 1,  
          p ≤ 
            n - 1,  p++,  
For[ i = p + 1,  
              i ≤ n,  i++,  
m = 
                ((Subscript[A, [[i, p]]])/(Subscript[A, [[p, p]]])); 
(Subscript[L, [[i, p]]]) = 
                m;  
(Subscript[A, [[i]]])  =  
                (Subscript[A, [[i]]])  -  m (Subscript[A, [[p]]]); 
U  =  A;  
Print[
                MatrixForm[U]];  ]; ]; 
        Print[N[ 
            BackSub[U, Transpose[U][[n + 1]], n]]] ]

\(\(\(LandU[A0_, 
        n_]\)\(:=\)\(
\)\(Module[\ {A = A0, i, 
          p}, \ 
U = A0; \ 
L = 
          IdentityMatrix[n]; \ 
Print[MatrixForm[L], 
          MatrixForm[U], " = ", 
          MatrixForm[A0]]; \ 
For[\ p = 1, \ 
          p ≤ 
            n - 1, \ \(p++\), \ 
\(For[\ i = p + 1, \ 
              i ≤ n, \ \(i++\), \ 
m = 
                A\_\(\([[\)\(i, p\)\(]]\)\)\/A\_\(\([[\)\(p, p\)\(]]\)\); 
L\_\(\([[\)\(i, p\)\(]]\)\) = 
                m; \ 
A\_\(\([[\)\(i\)\(]]\)\)\  = \ 
                A\_\(\([[\)\(i\)\(]]\)\)\  - \ m\ A\_\(\([[\)\(p\)\(]]\)\); 
U\  = \ A; \ 
Print[MatrixForm[L], 
                MatrixForm[U], " = ", 
                MatrixForm[A0]];\ \ ];\)\ \ ];\ \ ]\)\(\ \)\)\)