Download Solving Systems of Linear Algebraic Equations and more Exercises Reinforced Concrete Design in PDF only on Docsity!
TECHNOLOGICAL INSTITUTE OF THE PHILIPINES
938 Aurora Blvd., Cubao, Quezon City CE 007 (Numerical Solutions to CE Problems) Machine Problem 2. Numerical Solutions of Equations and Systems of Equations Submitted by: Gaddi, Marken John B. 1811969 ACADEMIC INTEGRITY PLEDGE I swear on my honor that I did not use any inappropriate aid, nor give such to others, in accomplishing this coursework. I understand that cheating and/or plagiarism is a major offense, as stated in TIP Memorandum No. P-04, s. 2017-2018, and that I will be sanctioned appropriately once I have committed such acts. Marken John B. Gaddi 1811969
MACHINE EXERISE:
a. Make an Algorithm for solving systems of linear algebraic equations using the methods below:
- Gauss Elimination Step 1. Start Step 2. Declare the variables and read the order of the matrix n. Step 3. Take the coefficients of the linear equation as: Do for k=1 to n Do for j=1 to n+ Read a[k][j] End for j End for k Step 4. Do for k=1 to n- 1 Do for i=k+1 to n Do for j=k+1 to n+ a[i][j] = a[i][j] – a[i][k] /a[k][k] * a[k][j] End for j End for i End for k Step 5. Compute x[n] = a[n][n+1]/a[n][n] Step 6. Do for k=n-1 to 1 sum = 0 Do for j=k+1 to n sum = sum + a[k][j] * x[j] End for j x[k] = 1/a[k][k] * (a[k][n+1] – sum) End for k Step 7: Display the result x[k] Step 8. Stop
Step 7. Print Array X as the solution Step 8. Stop
2. Crout’s Decomposition Crout's method decomposes a nonsingular n × n matrix A into the product of an n×n lower triangular matrix L and an n×n unit upper triangular matrix U. A unit triangular matrix is a triangular matrix with 1's along the diagonal. Crout's algorithm proceeds as follows: Evaluate the following pair of expressions for k = 0,... , n-1; And Note that Lik = 0 for k > i, Uik = 0 for k < i, and Ukk = 1 for k = 0, …, n-1. The matrix U forms a unit upper triangular matrix, and the matrix L forms a lower triangular matrix. The matrix A = LU. After the LU decomposition of A is performed, the solution to the system of linear equations A x = L U x = B is solved by solving the system of linear equations L y = B by forward substitution for y, and then solving the system of linear equations U x = y by backward substitution for x. 3. Cholesky’s Decomposition
- Iterative Methods 1. Gauss-Jacobi 2. Gauss-Seidel
4. Conjugate Gradient The algorithm is detailed below for solving Ax = b where A is a real, symmetric, positive-definite matrix. The input vector x0 can be an approximate initial solution or 0. b. Write a program for solving the system Ax = b by
- Gaussian Elimination Algorithm #include #include using namespace std; int main() {int n,i,j,k; cout << setprecision(3)<<fixed; cout.setf(ios::fixed); cout<<"Enter the number of unknowns: "; cin>>n;
float a[n][n+1],x[n]; cout<<"Enter elements of augmented matrix (row-wise): "<< endl; for (i=0;i<n;i++) for (j=0;j<=n;j++) cin>>a[i][j]; for (i=0;i<n-1;i++) for (k=i+1;k<n;k++) {double t=a[k][i]/a[i][i]; for (j=0;j<=n;j++) a[k][j]=a[k][j]-ta[i][j];} cout<<"\n\nThe matrix after elementary row operations to conduct gauss- elimination is as follows:\n"; for (i=0;i<n;i++) {for (j=0;j<=n;j++) cout<<a[i][j]<<setw(16); cout<<"\n";} for (i=n-1;i>=0;i--) {x[i]=a[i][n]; for (j=i+1;j<n;j++) if (j!=i) x[i]=x[i]-a[i][j]x[j]; x[i]=x[i]/a[i][i];} cout<<"\nThe solution in the system:\n"; for (i=0;i<n;i++) cout<<x[i]<<endl; return 0; }
for (w=0; w<z; w++) for (q=0; q<z; q++) cin >> a[w][q]; LUdecomposition(a, l, u, z); cout << "L DECOMPOSITION..."<<endl; for (w=0; w<z; w++) { for (q=0; q<z; q++) { cout<<l[w][q]<<" ";} cout << endl} cout << "U DECOMPOSITION..."<<endl; for (w=0; w<z; w++) { for (q=0; q<z; q++) { cout<<u[w][q]<<" ";} cout << endl;} return 0;}
- Iterative Methods 2. Conjugate Gradient vec conjugateGradientSolver( const matrix &A, const vec &B ) { double TOLERANCE = 1.0e- 10 int n = A.size(); vec X( n, 0.0 ); vec R = B; vec P = R; int k = 0; while ( k < n ) {vec Rold = R vec AP = matrixTimesVector( A, P ); double alpha = innerProduct( R, R ) / max( innerProduct( P, AP ), NEARZERO ); X = vectorCombination( 1.0, X, alpha, P );
R = vectorCombination( 1.0, R, - alpha, AP ); if ( vectorNorm( R ) < TOLERANCE ) break; double beta = innerProduct( R, R ) / max( innerProduct( Rold, Rold ), NEARZERO ); P = vectorCombination( 1.0, R, beta, P ); k++;} return X;} MACHINE PROBLEM:
1. (Hilbert matrix) Suppose A=[hij],hij=1i+j−1,i,j=1,...,n. Also, let b=Ax, xi=1, 1, ... , n, n=5, 10, 20. - Solve for Ax' = b (given any value for x(0)), and compare the results for different values of n with x* = 1 (i.e., one vector). Discuss in detail the results obtained for each case. - Note: Hilbert matrices are notoriously ill-conditioned matrices. Because of this, you may get almost zero pivot during the process, and as a result, you may not be able to solve the linear system for some large n.
- Gauss-Jordan A= [1.44, - 0.36, 5.52, 0.00; - 0.36, 10.33, - 7.78, 0.00; 5.52, - 7.78, 28.40,9.00; 0.00, 0.00, 9.00, 61.00] b = [0.04; - 2.15; 0; 0.88] Ab = [A,b] X = rref(Ab) x = X(:,5)
X = (U^-1)*Y
2. Crout’s Decomposition A= [1.44, - 0.36, 5.52, 0.00; - 0.36, 10.33, - 7.78, 0.00; 5.52, - 7.78, 28.40,9.00; 0.00, 0.00, 9.00, 61.00] b = [0.04; - 2.15; 0; 0.88] L = chol(A,'lower') U = chol(A,'upper') Y = (L^-1)b X = (U^-1)Y
- Gauss-Jacobi** A=[1.44, - 0.36, 5.52, 0.00; - 0.36, 10.33, - 7.78, 0.00; 5.52, - 7.78, 28.40,9.00; 0.00, 0.00, 9.00, 61.00] b=[0.04; - 2.15; 0; 0.88]' x=[0 0 0 0]' n=size(x,1); normVal=Inf; tol=1e-5; itr=0; while normVal>tol xold=x; for i=1:n sigma=0; for j=1:n if j~=i sigma=sigma+A(i,j)x(j); end end x(i)=(1/A(i,i))(b(i)-sigma); end itr=itr+1; normVal=abs(xold-x); end
fprintf('Solution of the system is : \n%f\n%f\n%f\n%f in %d iterations',x,itr);
2. Gauss-Seidel A=[1.44, - 0.36, 5.52, 0.00; - 0.36, 10.33, - 7.78, 0.00; 5.52, - 7.78, 28.40,9.00; 0.00, 0.00, 9.00, 61.00] b=[0.04; - 2.15; 0; 0.88]' x=[0 0 0 0]' n=size(x,1); normVal=Inf; tol=1e-5; itr=0; while normVal>tol x_old=x;
