1 / 10

Lecture 13 Sparse matrix computation

Lecture 13 Sparse matrix computation. Matlab function: spdiags. Download spdiags01.m And spdiags02.m. n=5; e = ones(n,1); A = spdiags([e -2*e 3*e], -1:1, n, n) full(A) >> spdiags01 A = (1,1) -2 (2,1) 1 (1,2) 3 (2,2) -2 (3,2) 1

gregoryfigueroa
Download Presentation

Lecture 13 Sparse matrix computation

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. Lecture 13 Sparse matrix computation Matlab function: spdiags Download spdiags01.m And spdiags02.m

  2. n=5; e = ones(n,1); A = spdiags([e -2*e 3*e], -1:1, n, n) full(A) >> spdiags01 A = (1,1) -2 (2,1) 1 (1,2) 3 (2,2) -2 (3,2) 1 (2,3) 3 (3,3) -2 (4,3) 1 (3,4) 3 (4,4) -2 (5,4) 1 (4,5) 3 (5,5) -2 ans = -2 3 0 0 0 1 -2 3 0 0 0 1 -2 3 0 0 0 1 -2 3 0 0 0 1 -2

  3. n=5; e1 = [1; 2; 3; 4; 5]; e2=[-1; -2; -3; -4; -5]; e3=[11; 12; 13; 14; 15]; A = spdiags([e1 e2 e3], -1:1, n, n); full(A) >> spdiags02 … ans = -1 12 0 0 0 1 -2 13 0 0 0 2 -3 14 0 0 0 3 -4 15 0 0 0 4 -5 full(A) – display matrix in usual format

  4. Define diagonals –2, 0 and 1 A = spdiags([e1 e2 e3], [-2 0 1], n, n) >> full(A) ans = -1 12 0 0 0 0 -2 13 0 0 1 0 -3 14 0 0 2 0 -4 15 0 0 3 0 -5

  5. Convert from sparse matrix to reduced form Taking into account nonzero diagonals only [B,d] = spdiags(A) extracts all nonzero diagonals from the m-by-n matrix A. B is a min(m,n)-by-p matrix whose columns are the p nonzero diagonals of A. d is a vector of length p whose integer components specify the diagonals in A. Example: >> [B d]=spdiags(A) B = 1 -1 0 2 -2 12 3 -3 13 0 -4 14 0 -5 15 d = -2 0 1

  6. The same in standard representation of matrix >> A0=full(A) A0 = -1 12 0 0 0 0 -2 13 0 0 1 0 -3 14 0 0 2 0 -4 15 0 0 3 0 -5 >> [B1,d1]=spdiags(A0) B1 = 1 -1 0 2 -2 12 3 -3 13 0 -4 14 0 -5 15 d1 = -2 0 1

  7. If a column of B is longer than the diagonal it's replacing, spdiags takes elements of super-diagonals from the lower part of the column of B, and elements of sub-diagonals from the upper part of the column of B. >> full(A) ans = -1 12 0 0 0 0 -2 13 0 0 1 0 -3 14 0 0 2 0 -4 15 0 0 3 0 -5 B = 1 -1 0 2 -2 12 3 -3 13 0 -4 14 0 -5 15

  8. Inclass 1 Modify Use spdiags to create the following sparse matrix: 1 2 0 0 0 0 0 0 0 0 0 1 2 0 0 0 0 0 0 0 -1 0 1 2 0 0 0 0 0 0 0 -1 0 1 2 0 0 0 0 0 0 0 -1 0 1 2 0 0 0 0 0 0 0 -1 0 1 2 0 0 0 0 0 0 0 -1 0 1 2 0 0 0 0 0 0 0 -1 0 1 2 0 0 0 0 0 0 0 -1 0 1 2 0 0 0 0 0 0 0 -1 0 1

  9. Newton’s method for system of two nonlinear equations Download newton2d01.m Function newton2d01 Set initial (x,y) point, maximum number of iterations, and convergence tolerance. for i= 1:maxiter Perform Newton update; check for conver-gence; print iterates end Define functions f1(x,y), and f2(x,y) and JacobianJacob(x,y).

  10. function newton2d01 %solution of system of two nonlienar equations % f1(x,y) = exp(x)*y = 0 % f2(x,y) = x*cos(y) = 0 % using a Newton iteration. Calls function f1 and f2 and % Jacobian matrix Jacob. % % Setup x = 0.2; y = .5; maxiter = 25; tol = 1.e-6; err = inf; iter = 0; fprintf('iter = %3.0f, x = %6.4f, y = %6.4f, \n', 0, x, y); % Newton update. while err > tol & iter <= maxiter J = Jacob(x,y); rhs = -[f1(x,y);f2(x,y)]; dxy = J\rhs; x = x + dxy(1); y = y + dxy(2); err = norm([f1(x,y);f2(x,y)],Inf); iter = iter+1; disp(['Iteration number = ',num2str(iter),' x=',num2str(x),' y=',num2str(y),' Backward error=',num2str(err)]); end function z = f1(x,y) z = exp(x)*y; function z = f2(x,y) z = x*cos(y); function xx = Jacob(x,y) xx = [exp(x)*y, exp(x); cos(y), -x*sin(y)];

More Related