复旦大学：《Matlab》（ppt实验英文版）Chapter 6 Linear systems

Introduction to scientific Computing A Matrix Vector Approach Using Matlab Written by Charles FVan Loan 陈文斌 Wbchen(fudan. edu. cn 复日大学

Introduction to Scientific Computing -- A Matrix Vector Approach Using Matlab Written by Charles F.Van Loan 陈 文 斌 Wbchen@fudan.edu.cn 复旦大学

Chapter 6 Linear systems T trangular Problems Banded problems Full problems A nalySIS

Chapter 6 Linear systems • Triangular Problems • Banded Problems • Full Problems • Analysis

Triangular problems 7100x1「b 22 0‖x 31132 x2=(b2-l21x1)/l2 23 Forward substitution X ∑

          =                     3 2 1 3 2 1 31 32 33 21 22 11 0 0 0 b b b x x x l l l l l l 3 3 3 1 1 3 2 2 2 3 2 2 2 1 1 2 2 1 1 1 1 ( )/ ( )/ / x b l x l x l x b l x l x b l = − − = − = Triangular Problems ii i j i i ij j x b l x /l 1 1         = − − = Forward substitution

for i=1:n x(1)=b(1) rj=1:i-1 x(1)=x(1)-L(Jj)*x( enc x(1=x(i/L(ii floops end x(1)=b(1)/L(1,1); for i=2: n Row oriented x()-(b(i)-L(i,1;-1)*x(1:i-1)/L(i) end

for i=1:n x(i)=b(i); for j=1:i-1 x(i)=x(i)-L(i,j)*x(j); end x(i)=x(i)/L(i,i); end x(1)=b(1)/L(1,1); for i=2:n x(i)=(b(i)-L(i,1:i-1)*x(1:i-1))/L(i,i); end Row oriented floops

Column-oriented 200x1|6 150 798 50 2 98 16 -1/5 8x3=-16-9(-1/5)國x3=-71/40

Column-oriented           =                     5 2 6 7 9 8 1 5 0 2 0 0 3 2 1 x x x       − − =      −      =            16 1 7 1 3 5 2 9 8 5 0 3 2 x x x1 = 3 x2 = −1/5 8 16 9( 1/ 5) x3 = − − − x3 = −71/ 40

22 0 X 21 32 0 x131 b(2:n)-x1L(2:n,1) n3 bn -x,. nn Saxpy b(j+1: n)<b(j+l: n)-, l(i+l: n,j)

(2 : ) (2 : ,1) 0 0 0 1 1 1 3 1 3 1 2 1 2 1 3 2 2 3 3 2 3 3 2 2 b n x L n b x l b x l b x l x x x l l l l l l n n n n n n n = −             − − − =                                  1 1 11 x = b /l j j jj x = b / l b( j 1: n) b( j 1: n) x L( j 1: n, j) +  + − j + saxpy

function x- LTriSol(L, b) n= len ngth(b) X=zeros(n forj=l: n-1 x=b(/Lg,j) b(+1n)=b(j+1n)-L(+1:nj)*x(j end x(n) =b(n/L(n,n) The version involves n2 flops, just like the row-oriented dot product version developed eariler

function x = LTriSol(L,b) n = length(b); x = zeros(n,1); for j=1:n-1 x(j) = b(j)/L(j,j); b(j+1:n) = b(j+1:n) - L(j+1:n,j)*x(j); end x(n) = b(n)/L(n,n); The version involves n 2 flops, just like the row-oriented, dot product version developed eariler

Backward Substitution 13 b 0 22 23 2 00 b 33 x2=(b2-l23x3)/2 3=(b1-4

Backward Substitution           =                     3 2 1 3 2 1 3 3 2 2 2 3 1 1 1 2 1 3 0 0 0 b b b x x x u u u u u u 3 1 1 2 2 1 3 3 1 1 2 2 2 3 3 2 2 3 3 3 3 ( )/ ( )/ / x b u x u x u x b u x u x b u = − − = − =

x(n)=b(n)/U(n,n) fori=n-1:-1:1 x()-(b()-U(i,i+1:n)*x(i+1n)/U(i,1) ent Backward Substitution function x= UTriSol(U, b) n= length(b) ⅹ= zeros for j=n: 1: 2 Column-oriented Saxpy x()=bj)U〔j); b(1」-1)=b(1」-1)-X(j)*U(1j-1j) end x(1)=b(1)/U(1,1)

x(n)=b(n)/U(n,n); for i=n-1:-1:1 x(i)=(b(i)-U(i,i+1:n)*x(i+1:n))/U(i,i) end function x = UTriSol(U,b) n = length(b); x = zeros(n,1); for j=n:-1:2 x(j) = b(j)/U(j,j); b(1:j-1) = b(1:j-1) - x(j)*U(1:j-1,j); end x(1) = b(1)/U(1,1); Column-oriented saxpy Backward Substitution

Multiple right-Hand sides lX= B LX( k)=B(k) X=zeros(nr) for kl. r ⅹ= zeros( X( k)-LTriSol(L, B( k);For k-1: r end for j=l: n-1 X(, k)=B0, kLG j BG+l: n, k)=BG+l n, k)-LO+l: n,j)*XGk) ene X(n, k)=B(n, kL(n,n) End

Multiple Right-Hand Sides LX = B LX (:, k) = B(:, k) X=zeros(n,r); for k=1:r X(:,k)=LTriSol(L,B(:,k)); end x = zeros(n,1); For k=1:r for j=1:n-1 X(j,k) = B(j,k)/L(j,j); B(j+1:n,k) = B(j+1:n,k) - L(j+1:n,j)*X(j,k); end X(n,k) =B(n,k)/L(n,n); End