CH-01 (523167), страница 5
Текст из файла (страница 5)
Note that xN is already found equal to vN after the Nthnormalization.Program Gauss listed below in both QuickBASIC and FORTRAN languagesis developed for interactive specification of the number of unknowns, N, and thevalues of the elements of [C] and {V}. It proceeds to solve for {X} and prints outthe computed values. Sample applications of both languages are provided immediately following the listed programs.A subroutine Gauss.Sub is also made available for the frequent need in theother computer programs which require the solution of matrix equations.© 2001 by CRC Press LLCQUICKBASIC VERSIONSample Application© 2001 by CRC Press LLCFORTRAN VERSION© 2001 by CRC Press LLCSample ApplicationGAUSS-JORDAN METHODOne slight modification of the elimination step will make the backward substitution steps completely unnecessary. That is, during the elimination of the xi termsfrom the linear algebraic equations except the ith one, Equations 25 and 26 shouldbe applied for k equal to 1 through N and excluding k = i.
For example, the x3 termsshould be eliminated from the first, second, fourth through Nth equations. In thismanner, after the Nth normalization, [C] becomes an identity matrix and {V} willhave the elements of the required solution {X}. This modified method is calledGauss-Jordan method.A subroutine called GauJor is made available based on the above argument. Inthis subroutine, a block of statements are also added to include the consideration ofthe pivoting technique which is required if ci,i = 0. The normalization steps,Equations 49 and 50, cannot be implemented if ci,i is equal to zero. For such asituation, a search for a nonzero ci,k is necessary for i = k + 1,k + 2,…,N.
That is,to find in the kth column of [C] and below the kth row a nonzero element. Oncethis nonzero ci,k is found, then we can then interchange the ith and kth rows of [C]and {V} to allow the normalization steps to be implemented; if no nonzero ci,k canbe found then [C] is singular because the determinant of [C] is equal to zero! Thiscan be explained by the fact that when ck,k = 0 and no pivoting is possible and thedeterminant D of [C] can be calculated by the formula:ND = c1,1c2,2 … c k ,k … c N,N =∏ck =1where indicates a product of all listed factors.© 2001 by CRC Press LLCk ,k(28)A subroutine has been written based on the Gauss-Jordan method and calledGauJor.Sub.
Both QuickBASIC and FORTRAN versions are made available andthey are listed below.QUICKBASIC VERSION© 2001 by CRC Press LLCFORTRAN VERSION© 2001 by CRC Press LLCSample ApplicationsThe same problem previously solved by the program Gauss has been used againbut solved by application of subroutine GauJor. The results obtained with the QuickBASIC and FORTRAN versions are listed, in that order, below:MATLAB APPLICATIONSFor solving the vector {X} from the matrix equation [C]{X} = {R} when boththe coefficient matrix [C] and the right-hand side vector {R} are specified, MATLABsimply requires [C] and {R} to be interactively inputted and then uses a statementX = C\R to obtain the solution vector {X} by multiplying the vector {R} on the leftof the inverse of [C] or dividing {R} on the left by [C]. More details are discussedin the program MatxAlgb.
Here, for providing more examples in MATLAB applications, a m file called GauJor.m is presented below as a companion of the FORTRAN and QuickBASIC versions:© 2001 by CRC Press LLCThis file GauJor.m should then be added into MATLAB. As an example ofinteractive application of this m file, the sample problem used in the FORTRANand QuickBASIC versions is again solved by specifying the coefficient matrix [C]and the right hand side vector {R} to obtain the resulting display as follows:The results of the vector {X} and determinant D for the coefficient matrix [C]are same as obtained before.MATHEMATICA APPLICATIONSFor solving a system of linear algebraic equations which has been arranged inmatrix form as [A]{X} = {R}, Mathematica’s function LinearSolve can be applied© 2001 by CRC Press LLCto solve for {X} when the coefficient matrix [A] and the right-hand side vector {R}are both provided.
The following is an example of interactive application:In[1]: = A = {{3,6,14},{6,14,36},{14,36,98}}Out[1]: ={{3, 6, 14}, {6, 14, 36}, {14, 36, 98}}In[2]: = MatrixForm[A]Out[2]//MatrixForm: =361461436143698In[3]: = R = {9,20,48}Out[3]: ={9, 20, 48}In[4]: = LinearSolve[A,R]Out[4]: ={–9,13,–3}Output[2] and Output[1] demonstrate the difference in display of matrix [A]when MatrixeForm is requested, or, not requested, respectively.
It shows that withoutrequesting of MatrixForm, some screen space saving can be gained. Output[4] givesthe solution {X} = [–9 13 –3]T for the matrix equation [A]{X} = {R} where thecoefficient matix [A] and vector {R} are provided by Input[1] and Input[3], respectively.1.5 MATRIX INVERSION, DETERMINANT,AND PROGRAM MatxInvDGiven a square matrix [C] of order N, its inverse as [C]–1 of the same order is definedby the equation:[C][C]−1 = [C]−1[C] = [I](1)where [I] is an identity matrix having elements equal to one along its main diagonaland equal to zero elsewhere.
That is:.1 00 1 0[I] = . . ..0 0© 2001 by CRC Press LLC............001(2)To find [C]–1, let cij and dij be the elements at the ith row and jth column of thematrices [C] and [C]–1, respectively. Both i and j range from 1 to N.
Furthermore,let {Dj} and {Ij} be the jth column of the matrices [C]–1 and [I], respectively. It iseasy to observe that {Ij} has elements all equal to zero except the one in the jth rowwhich is equal to unity. Also,{D } = [d djlj 2 j… d Nj]T(3)and[C]−1 = [D1D2 … DN ]T(4)Based on the rules of matrix multiplication, Equation 1 can be interpreted as[C]{D1} = {I1}, [C]{D2} = {I2}, …, and [C]{DN} = {IN}.
This indicates that programGauss can be successively employed N times by using the same coefficient matrix[C] and the vectors {Ii} to find the vectors {Di} for i = 1,2,…,N. Program MatxInvDis developed with this concept by modifying the program Gauss.
It is listed belowalong with a sample interactive run.QUICKBASIC VERSION© 2001 by CRC Press LLCSample ApplicationFORTRAN VERSION© 2001 by CRC Press LLC© 2001 by CRC Press LLCSample ApplicationsMATLAB APPLICATIONMATLAB offers very simple matrix operations. For example, matrix inversioncan be implemented as:To check if the obtained inversion indeed satisfies the equation [A}[A]–1 = [I]where [I] is the identity matrix, we enter:Once [A]–1 becomes available, we can solve the vector {X} in the matrix equation[A]{X} = {R} if {R} is prescribed, namely {X} = [A]–1{R}.
For example, may entera {R} vector and find {X} such as:© 2001 by CRC Press LLCMATHEMATICA APPLICATIONSMathematica has a function called Inverse for inversion of a matrix. Let usreuse the matrix A that we have entered in earlier examples and continue to demonstrate the application of Inverse:In[1]: = A = {{1,2},{3,4}}; MatrixForm[A]Out[1]//MatrixForm =1234In[2]: = B = {{5,6},{7,8}}; MatrixForm[B]Out[2]//MatrixForm =5678In[3]: = MatrixForm[A + B]Out[3]//MatrixForm =681012In[4]: = Dif = A-B; MatrixForm[Dif]Out[4]//MatrixForm =–4–4–4–4In[5]: = AT = Transpose[A]; MatrixForm[AT]Out[5]//MatrixForm =1324In[6]: = Ainv = Inverse[A]; MatrixForm[Ainv]Out[6]//MatrixForm =–213 1−2 2© 2001 by CRC Press LLCTo verify whether or not the inverse matrix Ainv obtained in Output[6] indeedsatisfies the equations [A][A]–1 = [I] which is the identity matrix, we apply Mathematica for matrix multiplication:In[7]: = Iden = A.Ainv; MatrixForm[Iden]Out[7]//MatrixForm =1001A dot is to separate the two matrices A and Ainv which is to be multiplied in thatorder.
Output[7] proves that the computed matrix, Ainv, is the inverse of A! It shouldbe noted that D and I are two reserved variables in Mathematica for the determinantof a matrix and the identity matrix. In their places, here Dif and Iden are adopted,respectively.
For further testing, we show that [A][A]T is a symmetric matrix:In[8]: = S = A.AT; MatrixForm[S]Out[8]//MatrixForm =5111125And, the unknown vector {X} in the matrix equation [A]{X} = {R} can besolved easily if {R} is given and [A]–1 are available:In[9]: = R = {13,31}; X = Ainv.ROut[9] = {5, 4}The solution of x1 = 5 and x2 = 4 do satisfy the equations x1 + 2x2 = 13 and 3x1+ 4x2 = 31.TRANSFORMATIONOFCOORDINATE SYSTEMS, ROTATION,ANDANIMATIONMatrix algebra can be effectively applied for transformation of coordinate systems. When the cartesian coordinate system, x-y-z, is rotated by an angle z aboutthe z-axis to arrive at the system x-y-z as shown in Figure 2, where z and z axescoincide and directed outward normal to the plane of paper, the new coordinates ofa typical point P whose coordinates are (xP,yP,zP) can be easily obtained as follows:x ′P = OP cos(θ P − θ z ) = (OP cos θ P ) cos θ z + (OP sin θ P ) sin θ z= x P cos θ z + y p sin θ zy ′P = OP sin(θ P − θ z ) = (OP sin θ P ) cos θ z − (OP cos θ P ) sin θ z= x p sin θ z + y p sin θ z© 2001 by CRC Press LLCFIGURE 2.
The cartesian coordinate system, x-y-z, is rotated by an angle z about the zaxis to arrive at the system x-y-z.andz′P = z PIn matrix notation, we may define {P} = [xP yP zP]T and {P'} = [xP' yP' zP']T andwrite the above equations as {P'} = [Tz]{P} where the transformation matrix for arotation of z-axis by z is: cos θz[Tz ] = − sin θz 0sin θzcos θz000 1 (5)In a similar manner, it can be shown that the transformation matrices for rotatingabout the x- and y-axes by angles x and y, respectively, are:1[Tx ] = 000cos θx− sin θx0 sin θx cos θx (6)− sin θ y 0 cos θ y (7)andcos θ yTy = 0 sin θy[ ]© 2001 by CRC Press LLC010FIGURE 3. Point P whose coordinates are (xP,yP,zP) is rotated to the point P' by a rotationof z.It is interesting to note that if a point P whose coordinates are (xP,yP,zP) is rotatedto the point P' by a rotation of z as shown in Figure 3, the coordinates of P' canbe easily obtained by the formula {P'} = [Rz]{P} where [Rz] = [Tz]T.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
