> read("/afs/eos.ncsu.edu/users/k/kaltofen/www/courses/LinAlgebra/Maple/ > initlib.mpl"); libname := /afs/eos.ncsu.edu/users/k/kaltofen/www/courses/LinAlg\ ebra/Maple, "/afs/bp.ncsu.edu/dist/maple551/update", "/afs/bp.ncsu.edu/dist/maple551/lib" > with(refpkg); [E_I, E_II, E_III, mydet, myinverse, mysolve, ref, xref] > with(linalg); [BlockDiagonal, GramSchmidt, JordanBlock, LUdecomp, QRdecomp, Wronskian, addcol, addrow, adj, adjoint, angle, augment, backsub, band, basis, bezout, blockmatrix, charmat, charpoly, cholesky, col, coldim, colspace, colspan, companion, concat, cond, copyinto, crossprod, curl, definite, delcols, delrows, det, diag, diverge, dotprod, eigenvals, eigenvalues, eigenvectors, eigenvects, entermatrix, equal, exponential, extend, ffgausselim, fibonacci, forwardsub, frobenius, gausselim, gaussjord, geneqns, genmatrix, grad, hadamard, hermite, hessian, hilbert, htranspose, ihermite, indexfunc, innerprod, intbasis, inverse, ismith, issimilar, iszero, jacobian, jordan, kernel, laplacian, leastsqrs, linsolve, matadd, 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, stackmatrix, submatrix, subvector, sumbasis, swapcol, swaprow, sylvester, toeplitz, trace, transpose, vandermonde, vecpotent, vectdim, vector, wronskian] > A := matrix(4,4); A := array(1 .. 4, 1 .. 4, []) > rnd:=rand(-100..100); rnd := proc() local t; global _seed; _seed := irem(427419669081*_seed, 999999999989); t := _seed; irem(t, 201) - 100 end > for i from 1 to 4 do for j from 1 to 4 do A[i,j] := rnd(); > od: od: > print(A); [-34 -68 -85 -38] [ ] [ 90 52 43 91] [ ] [ 30 90 -24 -52] [ ] [ 41 58 82 21] > b1 := vector([1,2,3,4]); > b1 := [1, 2, 3, 4] > b2 := vector([a,b,c,d]); b2 := [a, b, c, d] > Ainv := inverse(A); [391407 86003 11471 363985 ] [------- ------- ------- -------] [9699947 9699947 9699947 9699947] [ ] [-266519 -11971 146925 -248491] [------- ------- -------- -------] [9699947 9699947 19399894 9699947] Ainv := [ ] [ 60233 -47556 -50104 191002 ] [------- ------- ------- -------] [9699947 9699947 9699947 9699947] [ ] [-263271 50847 -29648 -308244] [------- ------- ------- -------] [9699947 9699947 9699947 9699947] > x1:=evalm(Ainv &* b1); [2053766 -2128075 578817 -1483497] x1 := [-------, --------, -------, --------] [9699947 19399894 9699947 9699947 ] > evalm(A &* x1); [1, 2, 3, 4] > x2 := evalm(Ainv &* b2); [391407 86003 11471 363985 x2 := [------- a + ------- b + ------- c + ------- d, [9699947 9699947 9699947 9699947 266519 11971 146925 248491 - ------- a - ------- b + -------- c - ------- d, 9699947 9699947 19399894 9699947 60233 47556 50104 191002 ------- a - ------- b - ------- c + ------- d, 9699947 9699947 9699947 9699947 263271 50847 29648 308244 ] - ------- a + ------- b - ------- c - ------- d] 9699947 9699947 9699947 9699947 ] > evalm(A &* x2); [a, b, c, d] > U := xref(A, 4, 'T', 'L'); [-34 -68 -85 -38 ] [ ] [ -163 ] [ 0 -128 -182 ---- ] [ 17 ] [ ] U := [ -4533 -95501 ] [ 0 0 ----- ------ ] [ 32 1088 ] [ ] [ -9699947] [ 0 0 0 --------] [ 308244 ] > print(L); [ 1 0 0 0] [ ] [-45 ] [--- 1 0 0] [17 ] [ ] [-15 -15 ] [--- --- 1 0] [17 64 ] [ ] [-41 -436 ] [--- 3/16 ---- 1] [34 4533 ] > evalm(L &* U); [-34 -68 -85 -38] [ ] [ 90 52 43 91] [ ] [ 30 90 -24 -52] [ ] [ 41 58 82 21] > print(A); [-34 -68 -85 -38] [ ] [ 90 52 43 91] [ ] [ 30 90 -24 -52] [ ] [ 41 58 82 21] > y1 := linsolve(L,b1); [ 79 5409 494499] y1 := [1, --, ----, ------] [ 17 1088 102748] > xx1 := linsolve(U,y1); [2053766 -2128075 578817 -1483497] xx1 := [-------, --------, -------, --------] [9699947 19399894 9699947 9699947 ] > print(x1); [2053766 -2128075 578817 -1483497] [-------, --------, -------, --------] [9699947 19399894 9699947 9699947 ] > y2 := linsolve(L, b2); [ 45 1635 15 y2 := [a, -- a + b, ---- a + c + -- b, [ 17 1088 64 87757 997 436 ] ------ a + d - ---- b + ---- c] 102748 6044 4533 ] > xx2 := linsolve(U,y2); [391407 86003 11471 363985 xx2 := [------- a + ------- b + ------- c + ------- d, [9699947 9699947 9699947 9699947 266519 11971 146925 248491 - ------- a - ------- b + -------- c - ------- d, 9699947 9699947 19399894 9699947 60233 47556 50104 191002 ------- a - ------- b - ------- c + ------- d, 9699947 9699947 9699947 9699947 263271 50847 29648 308244 ] - ------- a + ------- b - ------- c - ------- d] 9699947 9699947 9699947 9699947 ] > print(x2); [391407 86003 11471 363985 [------- a + ------- b + ------- c + ------- d, [9699947 9699947 9699947 9699947 266519 11971 146925 248491 - ------- a - ------- b + -------- c - ------- d, 9699947 9699947 19399894 9699947 60233 47556 50104 191002 ------- a - ------- b - ------- c + ------- d, 9699947 9699947 9699947 9699947 263271 50847 29648 308244 ] - ------- a + ------- b - ------- c - ------- d] 9699947 9699947 9699947 9699947 ] > print(A); [-34 -68 -85 -38] [ ] [ 90 52 43 91] [ ] [ 30 90 -24 -52] [ ] [ 41 58 82 21] > ?LUdecomp > Lf := map(evalf, L); [ 1. 0 0 0 ] [ ] [-2.647058824 1. 0 0 ] Lf := [ ] [-.8823529412 -.2343750000 1. 0 ] [ ] [-1.205882353 .1875000000 -.09618354291 1.] > Uf := map(evalf, U); [-34. -68. -85. -38. ] [ ] [ 0 -128. -182. -9.588235294] Uf := [ ] [ 0 0 -141.6562500 -87.77665441] [ ] [ 0 0 0 -31.46840490] > Ainvf := map(evalf, Ainv); Ainvf := [.04035145759 , .008866337105 , .001182583781 , .03752443183] [-.02747633570 , -.001234130455 , .007573494989 , -.02561776884] [.006209621558 , -.004902707200 , -.005165389048 , .01969103543] [-.02714148851 , .005241987405 , -.003056511546 , -.03177790559] > x1f := evalm(Ainvf &* b1); x1f := [.2117296104, -.1096951870, .05967218174, -.1529386707] > y1f := linsolve(Lf, b1); y1f := [1., 4.647058824, 4.971507353, 4.812736015] > xx1f := linsolve(Uf, y1f); xx1f := [.2117296106, -.1096951870, .05967218169, -.1529386707] > map(evalf, x1); [.2117296105, -.1096951870, .05967218171, -.1529386707] >