수치해석 Mini Project : 구의 중심 추정하기

기계시스템디자인공학과 김성환 교수

1. 구(sphere)의 방정식

공간상의 구의 방정식은 다음과 같다.

[eqn: (x-a) ^{2} +(y-b) ^{2} +(z-c) ^{2} =r ^{2}]

이때, 구의 중심은 [eqn: (a,`b,`c)], 구의 반지름은 [eqn: r] 이다.

2. 4개의 점을 지나는 구(sphere)를 구하기

구의 표면에 근사적으로 존재하는 점을 수천개 알고 있다. 이 구의 중심과 반지름을 Least-Square 방법으로 얻어내고자 한다. 이를 위해서는 아래 (식 1)의 n개의 식의 제곱의 합이 최소가 되도록 해야 한다.

[eqn: {cases{f _{1} (a,`b,`c,`r)`=`(x _{1} -a) ^{2} +(y _{1} -b) ^{2} +(z _{1} -c) ^{2} -r ^{2} =0&#f _{2} (a,`b,`c,`r)`=`(x _{2} -a) ^{2} +(y _{2} -b) ^{2} +(z _{2} -c) ^{2} -r ^{2} =0&#VDOTS &#f _{n} (a,`b,`c,`r)`=`(x _{n} -a) ^{2} +(y _{n} -b) ^{2} +(z _{n} -c) ^{2} -r ^{2} =0&}}]         (식 1)

여기서 [eqn: (x _{i} ,`y _{i} ,`z _{i} )], [eqn: i=1,`2,`` LDOTS `,`n] 는 주어진 점 좌표이고 [eqn: (a,`b,`c,`r)]은 구해야 할 방정식의 해이다.

3. Gauss-Newton 반복법에 의한 비선형 다변수 문제의 Least-Square 해

[eqn: (a,`b,`c,`r)]의 초기값을 정한다.

반복횟수를 0으로 한다.

다음과 같은 선형 연립방정식을 만든다.

[eqn: {bmatrix{{Partial f _{1}} over {Partial a}&{Partial f _{1}} over {Partial b}&{Partial f _{1}} over {Partial c}&{Partial f _{1}} over {Partial r}#{Partial f _{2}} over {Partial a}&{Partial f _{2}} over {Partial b}&{Partial f _{2}} over {Partial c}&{Partial f _{2}} over {Partial r}#VDOTS &VDOTS &VDOTS &VDOTS #{Partial f _{n}} over {Partial a}&{Partial f _{n}} over {Partial b}&{Partial f _{n}} over {Partial c}&{Partial f _{n}} over {Partial r}}}][eqn: LEFT { {matrix{DELTA a#DELTA b#DELTA c#DELTA r}} RIGHT }][eqn: = LEFT { {matrix{-f _{1}#-f _{2}#VDOTS #-f _{n}}} RIGHT }]

위 식을 다음과 같이 썼을 때

[eqn: LEFT [ A RIGHT ] LEFT { DELTA q RIGHT } = LEFT { -f RIGHT }]

이의 양변에 [eqn: LEFT [ A RIGHT ] ^{T}]를 곱해 다음과 같은 선형 시스템을 얻는다.

[eqn: LEFT [ A RIGHT ] ^{T} LEFT [ A RIGHT ] LEFT { DELTA q RIGHT } = LEFT [ A RIGHT ] ^{T} LEFT { -f RIGHT }]

이를 풀어 [eqn: LEFT { DELTA q RIGHT }] 를 구한다.

풀어진 해에서 모든 [eqn: LEFT | DELTA q RIGHT |] 가 [eqn: epsilon _{s}] 보다 작아지면 반복을 종료하고 결과를 출력한다. [eqn: epsilon _{s} =10 ^{-6}] 으로 한다.

반복횟수를 하나 증가시킨다. 이것이 최대 허용 반복횟수보다 커지면 종료하고 오류를 출력한다. 최대 반복횟수는 20으로 한다.

다음과 같이 변수를 갱신하고 단계 3으로 되돌아 간다.

[eqn: LEFT { q RIGHT } ^{j+1} = LEFT { q RIGHT } ^{j} + LEFT { DELTA q RIGHT }]

4. 도함수로 이루어진 Jacobian

편미분을 수행하여 Jacobian을 구하면 다음과 같다.

[eqn: LEFT [ A RIGHT ] = {bmatrix{{Partial f _{1}} over {Partial a}&{Partial f _{1}} over {Partial b}&{Partial f _{1}} over {Partial c}&{Partial f _{1}} over {Partial r}#{Partial f _{2}} over {Partial a}&{Partial f _{2}} over {Partial b}&{Partial f _{2}} over {Partial c}&{Partial f _{2}} over {Partial r}#VDOTS &VDOTS &VDOTS &VDOTS #{Partial f _{n}} over {Partial a}&{Partial f _{n}} over {Partial b}&{Partial f _{n}} over {Partial c}&{Partial f _{n}} over {Partial r}}} = {bmatrix{-2(x _{1} -a)&````-2(y _{1} -b)&````-2(z _{1} -c)&````-2r#-2(x _{2} -a)&````-2(y _{2} -b)&````-2(z _{2} -c)&````-2r#VDOTS &VDOTS &VDOTS &VDOTS #-2(x _{n} -a)&````-2(y _{n} -b)&````-2(z _{n} -c)&````-2r}}]

5. 문제

다음과 같은 근사점들이 주어져 있다. 개수에 따라 3개의 파일로 주어진다.

(1) 100개의 데이터 점  : “points_100.csv.txt”    파일로 주어짐

(2) 1000개의 데이터 점 : “points_1000.csv.txt”   파일로 주어짐

(3) 5000개의 데이터 점 : “points_5000.csv.txt”   파일로 주어짐

이들 점은 (a, b, c, r) = (15, 25, 20, 12)의 구 표면 근처에 있는 점들이다.

[문제 1] 초기값 (a, b, c, r)을 다음과 같이 변화시키며 실행시켜보고 무슨 일이 일어 났는지 써라.

[문제 2] Gauss 소거법의 연산횟수는 [eqn: {2n ^{3}} over {3} +O(n ^{2} )] 이다. 계산시간과 이를 비교하여 논의하라.

번호	초기값	구해진 해	반복횟수	실행시간
(1) 100개 데이터	Qini=[10 30 15 7]
(2) 1000개 데이터
(3) 5000개 데이터

5. MATLAB 프로그램

Sphere.m

function [Q, iter, sec] = SphereApproximation()

% Least-Square solution of sphere

% inital value of [a b c r]

Qini=[10 30 15 7];

% 구면 점 데이터 파일 얻기

[file,path] = uigetfile('*.txt');

if isequal(file,0)

disp('User selected Cancel');

return;

end

% 파일 읽기

fname = fullfile(path,file);

xyz = csvread(fname);

xp = xyz(:,1);

yp = xyz(:,2);

zp = xyz(:,3);

% 계산

[Q, iter, sec] = gauss_newton(xp, yp, zp, Qini);

end

%==========================================================

%   Gauss-Newton method

%==========================================================

function [Q, iter, sec] = gauss_newton(x, y, z, Q)

tic;

nfunc = size(x, 1);

iter = 0;

while(1)

a = Q(1);

b = Q(2);

c = Q(3);

r = Q(4);

% A matrix (Jacobian)

A = zeros(nfunc, 4);

for i = 1:nfunc

A(i, 1:4) = [             채울것            ];

end

% function value

f = zeros(nfunc, 1);

for i = 1:nfunc

f(i) =           채울것                 ;

end

% calculate Pseudo-Inverse solution {dq}

dq =           채울것                 ;

%update

Q = Q + dq';

iter = iter + 1;

% check convergence

if(max(abs(dq)) < 1.e-6), break, end

if(iter > 20), break, end

end

sec = toc;

end

(이상)