수치해석 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)를 구하기

공간에 4개의 점이 있고 이들이 한 평면위에 놓여있지 않다면 이를 지나는 구가 존재한다. 이 구의 중심을 기하학적으로 구하는 것은 쉽지 않다.

이를 대수적으로 구할 수 있는데, 이를 위해서는 다음 4개의 비선형 연립 방정식을 풀어야 한다.

[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&#f _{3} (a,`b,`c,`r)`=`(x _{3} -a) ^{2} +(y _{3} -b) ^{2} +(z _{3} -c) ^{2} -r ^{2} =0&#f _{4} (a,`b,`c,`r)`=`(x _{4} -a) ^{2} +(y _{4} -b) ^{2} +(z _{4} -c) ^{2} -r ^{2} =0&}}]         (식 1)

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

3. Newton-Raphson 반복법에 의한 비선형연립방정식의 풀이

[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}#{Partial f _{3}} over {Partial a}&{Partial f _{3}} over {Partial b}&{Partial f _{3}} over {Partial c}&{Partial f _{3}} over {Partial r}#{Partial f _{4}} over {Partial a}&{Partial f _{4}} over {Partial b}&{Partial f _{4}} over {Partial c}&{Partial f _{4}} over {Partial r}}}][eqn: LEFT { {matrix{DELTA a#DELTA b#DELTA c#DELTA r}} RIGHT }][eqn: = LEFT { {matrix{-f _{1}#-f _{2}#-f _{3}#-f _{4}}} RIGHT }]

위 식을 다음과 같이 썼을 때

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

이를 풀어 [eqn: LEFT { DELTA q RIGHT }] 를 구한다. (만일 [eqn: LEFT [ J 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 [ J 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}#{Partial f _{3}} over {Partial a}&{Partial f _{3}} over {Partial b}&{Partial f _{3}} over {Partial c}&{Partial f _{3}} over {Partial r}#{Partial f _{4}} over {Partial a}&{Partial f _{4}} over {Partial b}&{Partial f _{4}} over {Partial c}&{Partial f _{4}} 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#-2(x _{3} -a)&````-2(y _{3} -b)&````-2(z _{3} -c)&````-2r#-2(x _{4} -a)&````-2(y _{4} -b)&````-2(z _{4} -c)&````-2r}}]

5. 문제

주어진 4개의 점이 다음과 같다고 하자.

P1(x, y, z) = (26.547002   26.547002   26.547002)

P2(x, y, z) = (21.318489   21.318489   21.318489)

P3(x, y, z) = (20.014806   20.014806   20.014806)

P4(x, y, z) = (15.039961   15.039961   15.039961)

이 4점을 지나는 구의 중심은 (a, b, c) = (15, 25, 20)이고 반지름 r = 12이다.

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

번호	초기값	수렴했는가?	반복횟수	결과 검토 의견
(1)	Qini=[10 30 15 7]
(2)	Qini=[0 0 0 0]
(3)	Qini=[0 0 0 1]
(4)	Qini=[10 10 10 20]

5. MATLAB 프로그램

Sphere.m

function [Q, iter] = Sphere(Qini)

% Find sphere center and radius

% Qini : [a b c r]의 초기값

% 구면 위에 있어야 하는 주어진 점

xp = [26.547002 21.318489 20.014806 15.039961];

yp = [28.094010 29.424254 34.842110 25.205581];

zp = [21.045869 29.192533 24.688774 31.998172];

% 계산

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

end

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

%   Newton-Raphson iteration

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

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

iter = 0;

while(1)

% Jacobian 계산

% function value 계산

% calculate {dq}

% update Q & iter 증가

% check convergence

end

end

6. 다른 해법

지금까지는 비선형연립방정식을 Newton-Rapshon 방법으로 해결하는 과정을 논의하였다. 그런데 이 문제는 최적화 기법을 이용하여 해결될 수도 있다. (교재 318페이지 하단 참조)

즉 (식 1)의 [eqn: f _{1} ,`f _{2} ,`f _{3} ,`f _{4}]에 대해 [eqn: F= sum _{i} ^{4} LEFT [ f _{i} (a,`b,`c,`r) RIGHT ] ^{2}] 을 최소하는  [eqn: a,`b,`c,`r]를 구하면 된다.

최적화 해를 얻기 위해서는 MATLAB의 기본 내장함수 fminsearch()를 사용할 수 있다. (교재 210페이지 참조)

fminsearch( )를 사용하는 프로그램을 작성하여 실행시켜 보라.

(이상)