A=[4 0 -8; 1 2 3; -2 4 -1] A = 4 0 -8 1 2 3 -2 4 -1 % initialize U and L U=A U = 4 0 -8 1 2 3 -2 4 -1 L=eye(size(A)) L = 1 0 0 0 1 0 0 0 1 % Gauss elim on column 1 L(2,1)=U(2,1)/U(1,1) L = 1.0000 0 0 0.2500 1.0000 0 0 0 1.0000 U(2,:)=U(2,:)-L(2,1)*U(1,:) U = 4 0 -8 0 2 5 -2 4 -1 L(3,1)=U(3,1)/U(1,1) L = 1.0000 0 0 0.2500 1.0000 0 -0.5000 0 1.0000 U(3,:)=U(3,:)-L(3,1)*U(1,:) U = 4 0 -8 0 2 5 0 4 -5 L(3,2)=U(3,2)/U(2,2) L = 1.0000 0 0 0.2500 1.0000 0 -0.5000 2.0000 1.0000 U(3,:)=U(3,:)-L(3,2)*U(2,:) U = 4 0 -8 0 2 5 0 0 -15 L*U ans = 4 0 -8 1 2 3 -2 4 -1 b=[12; 2; 7] b = 12 2 7 y=forsub(L,b) y = 12 -1 15 x=backsub(U,y) x = 1 2 -1 clear A=[-10 7 6 5; 4 3 -1 -7; 1 2 6 8; -2 -4 5 3] A = -10 7 6 5 4 3 -1 -7 1 2 6 8 -2 -4 5 3 [L,U]=lu(A) L = 1.0000 0 0 0 -0.4000 1.0000 0 0 -0.1000 0.4655 1.0000 0 0.2000 -0.9310 0.8580 1.0000 U = -10.0000 7.0000 6.0000 5.0000 0 5.8000 1.4000 -5.0000 0 0 5.9483 10.8276 0 0 0 -11.9449 L*U ans = -10.0000 7.0000 6.0000 5.0000 4.0000 3.0000 -1.0000 -7.0000 1.0000 2.0000 6.0000 8.0000 -2.0000 -4.0000 5.0000 3.0000 clear A=[4 0 -8; 1 2 3; -2 4 -1] A = 4 0 -8 1 2 3 -2 4 -1 P=[0 1 0; 1 0 0; 0 0 1] P = 0 1 0 1 0 0 0 0 1 P*A ans = 1 2 3 4 0 -8 -2 4 -1 A([1 2],:)=A([2 1],:) A = 1 2 3 4 0 -8 -2 4 -1 clear A=[1 2 3 -2; 2 4 1 0; 3 3 2 5; -1 6 2 1] A = 1 2 3 -2 2 4 1 0 3 3 2 5 -1 6 2 1 % initialize U and L U=A U = 1 2 3 -2 2 4 1 0 3 3 2 5 -1 6 2 1 L=eye(size(A)) L = 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 % Gauss elim on column 1 L(2,1)=U(2,1)/U(1,1) L = 1 0 0 0 2 1 0 0 0 0 1 0 0 0 0 1 U(2,:)=U(2,:)-L(2,1)*U(1,:) U = 1 2 3 -2 0 0 -5 4 3 3 2 5 -1 6 2 1 L(3,1)=U(3,1)/U(1,1) L = 1 0 0 0 2 1 0 0 3 0 1 0 0 0 0 1 U(3,:)=U(3,:)-L(3,1)*U(1,:) U = 1 2 3 -2 0 0 -5 4 0 -3 -7 11 -1 6 2 1 L(4,1)=U(4,1)/U(1,1) L = 1 0 0 0 2 1 0 0 3 0 1 0 -1 0 0 1 U(4,:)=U(4,:)-L(4,1)*U(1,:) U = 1 2 3 -2 0 0 -5 4 0 -3 -7 11 0 8 5 -1 % switch rows 2 and 3 to avoid zero pivot U([2 3],:)=U([3 2],:) U = 1 2 3 -2 0 -3 -7 11 0 0 -5 4 0 8 5 -1 L([2 3],1)=L([3 2],1) L = 1 0 0 0 3 1 0 0 2 0 1 0 -1 0 0 1 % permutation matrix P=eye(size(A)); P([2 3],:)=P([3 2],:) P = 1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 1 % Gauss elim on new column 2 L(4,2)=U(4,2)/U(2,2) L = 1.0000 0 0 0 3.0000 1.0000 0 0 2.0000 0 1.0000 0 -1.0000 -2.6667 0 1.0000 U(4,:)=U(4,:)-L(4,2)*U(2,:) U = 1.0000 2.0000 3.0000 -2.0000 0 -3.0000 -7.0000 11.0000 0 0 -5.0000 4.0000 0 0 -13.6667 28.3333 % Gauss elim on column 3 L(4,3)=U(4,3)/U(3,3) L = 1.0000 0 0 0 3.0000 1.0000 0 0 2.0000 0 1.0000 0 -1.0000 -2.6667 2.7333 1.0000 U(4,:)=U(4,:)-L(4,3)*U(3,:) U = 1.0000 2.0000 3.0000 -2.0000 0 -3.0000 -7.0000 11.0000 0 0 -5.0000 4.0000 0 0 0 17.4000 L*U - P*A ans = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 clear A = [10 -7 0; -3 2.09999 6; 5 -1 5] A = 10.0000 -7.0000 0 -3.0000 2.1000 6.0000 5.0000 -1.0000 5.0000 U=A U = 10.0000 -7.0000 0 -3.0000 2.1000 6.0000 5.0000 -1.0000 5.0000 L=eye(size(A)) L = 1 0 0 0 1 0 0 0 1 % Gauss elim on column 1 L(2,1)=U(2,1)/U(1,1) L = 1.0000 0 0 -0.3000 1.0000 0 0 0 1.0000 U(2,:)=U(2,:)-L(2,1)*U(1,:) U = 10.0000 -7.0000 0 0 -0.0000 6.0000 5.0000 -1.0000 5.0000 L(3,1)=U(3,1)/U(1,1) L = 1.0000 0 0 -0.3000 1.0000 0 0.5000 0 1.0000 U(3,:)=U(3,:)-L(3,1)*U(1,:) U = 10.0000 -7.0000 0 0 -0.0000 6.0000 0 2.5000 5.0000 % Gauss elim on column 2 L(3,2)=U(3,2)/U(2,2) L = 1.0e+05 * 0.0000 0 0 -0.0000 0.0000 0 0.0000 -2.5000 0.0000 U(3,:)=U(3,:)-L(3,2)*U(2,:) U = 1.0e+06 * 0.0000 -0.0000 0 0 -0.0000 0.0000 0 0 1.5000 b=[7 3.90001 6] format long y=forsub(L,b) y = 1.0e+06 * 0.000007000000000 0.000006000010000 1.500004999990173 x=backsub(U,y) x = 0.000000000031086 -0.999999999955591 1.000000000000000 clear format long A = [10 -7 0; -3 2.09999 6; 5 -1 5]; U=A; L=eye(size(A)); % Gauss elim on column 1 L(2,1)=U(2,1)/U(1,1); U(2,:)=U(2,:)-L(2,1)*U(1,:); L(3,1)=U(3,1)/U(1,1); U(3,:)=U(3,:)-L(3,1)*U(1,:) U = 10.000000000000000 -7.000000000000000 0 0 -0.000010000000000 6.000000000000000 0 2.500000000000000 5.000000000000000 clear A = [10 -7 0; -3 2.09999 6; 5 -1 5] A = 10.000000000000000 -7.000000000000000 0 -3.000000000000000 2.099990000000000 6.000000000000000 5.000000000000000 -1.000000000000000 5.000000000000000 [L,U]=lu(A) L = 1.000000000000000 0 0 -0.300000000000000 -0.000004000000000 1.000000000000000 0.500000000000000 1.000000000000000 0 U = 10.000000000000000 -7.000000000000000 0 0 2.500000000000000 5.000000000000000 0 0 6.000020000000001 [L,U,P]=lu(A) L = 1.000000000000000 0 0 0.500000000000000 1.000000000000000 0 -0.300000000000000 -0.000004000000000 1.000000000000000 U = 10.000000000000000 -7.000000000000000 0 0 2.500000000000000 5.000000000000000 0 0 6.000020000000001 P = 1 0 0 0 0 1 0 1 0 b=[7; 3.90001;6]; y=forsub(L,P*b) y = 7.000000000000000 2.500000000000000 6.000020000000001 x=backsub(U,y) x = 0 -1 1