1. 선형시스템을 MATLAB으로 해결하는 방법

Gauss 소거법의 MATLAB실습 : 서울과학기술대학교 기계시스템디자인공학과 / 김성환교수

>> K = [150 -100 0; -100 150 -50; 0 -50 50]

>> mg = [588.6; 686.7; 784.8]

>> x = A \ mg                   % 바람직 함

>> x = inv(A) * mg              % 바람직하지 않음

2. GaussPivot 함수를 이해하기 위한 기본기

(1) max 함수의 사용법

>> y = [1; 3; 5; 2; 4];

>> [ymax, imax] = max(y)    % 최대값과 최대위치가 얻어짐

ymax = 5

imax = 3

(2) 행과 행의 교환

① 잘못된 방법

>> K = [150 -100 0; -100 150 -50; 0 -50 50]

>> K(1, :)             % 1행만의 부분 행렬

>> K(1, :) = K(3, :)   % 덮어씀. 1행이 소실됨. 잘못 !!

② 임시 저장변수를 쓰는 방법

>> K = [150 -100 0; -100 150 -50; 0 -50 50]

>> tmp = K(1, :)       % 1행을 tmp에 보관

>> K(1, :) = K(3, :)   % 1행을 3행으로 덮어씀

>> K(3, :) = tmp

③ 보다 좋은 방법

>> K = [150 -100 0; -100 150 -50; 0 -50 50]

>> K([1 3], :)                 % 1, 3행만의 부분 행렬

>> K([1 3], :) = K([3 1], :)   % 1, 3행을 한꺼번에 덮어씀

(3) 행렬의 부분곱

>> A = [0 2 3; 4 6 7; 2 -3 6]

>> x = [8; -3; 5]

>> A(2,2:3) * x(2:3)

[eqn: {bmatrix{a _{11}&a _{12}&a _{13}#a _{21}&a _{22}&a _{23}#a _{31}&a _{32}&a _{33}}} LEFT { {matrix{x _{1}#x _{2}#x _{3}}} RIGHT }]

3. GaussPivot 함수

함수를 입력

function x = GaussPivot(A,b)

% GaussPivot (A,b): Gauss elimination with partial pivoting

% input:

%   A = coefficient matrix

%   b = right hand side vector

% output:

%   x = solution vector

[m,n] = size(A);

if m ~= n, error('Matrix A must be square'); end

nb = n+1;

Aug =[A b];

% forward elimination

for k = 1:n-1

% partial pivoting

[big, i] = max(abs(Aug(k:n,k)));

ipr = i + k - 1;

if ipr  ~= k

% pivot the row

Aug([k,ipr],:) = Aug([ipr,k],:);

end

for i = k+1:n

factor = Aug(i,k)/Aug(k,k);

Aug(i,k:nb) = Aug(i,k:nb)-factor*Aug(k,k:nb);

end

end

% back substitution

x = zeros(n,1);

x(n) = Aug(n,nb)/Aug(n,n);

for i = n-1:-1:1

x(i) = (Aug(i,nb)-Aug(i,i+1:n)*x(i+1:n))/Aug(i,i);

end

실행

>> A = [0 2 3; 4 6 7; 2 -3 6]

>> b = [8; -3; 5]

>> GaussPivot(A, b)

Break Point(F12)를 설정한 후 Step(F10) 하면서 행렬의 변화를 살펴본다.

혹은 함수 중간에 disp() 함수를 추가하여 Aug 행렬의 변화를 출력하여 본다.

4. 도전과제

>> A = [1 2 3; 4 5 6; 7 8 9]

>> b = [8; -3; 5]

>> A \ b

>> GaussPivot(A, b)

위를 수행하면 오류를 만나게 된다.

A 행렬이 행렬식이 0인 (0에 가까운) 특이(singular) 행렬이기 때문이다.

GaussPivot 함수를 실행하면 이상한 값이 나온다.

특이행렬이라는 것을 어디서 알아낼 수 있을까?

어떻게 처리해야 할까?

이상.