수치해석 (Numerical Analysis)

Presentation on theme: "수치해석 (Numerical Analysis)"— Presentation transcript:

1 수치해석 (Numerical Analysis)
다변수 방정식과 함수 (Part 1)

2 In this chapter … 다변수 방정식과 함수 변수가 두 개 이상인 함수, 예를 들어, 의 해(f(x,y,z)=0으로 하는 (x,y,z)의 값)를 수치 해석적으로 구하는 방법을 다룬다. 2차 함수인 f(x,y)를 집중적으로 다룬다. (3차 이상 확장 가능) We will cover … 이차원 이분 격자법 (Bisection Grid) 영점 곡선 추적 (Zero-Curve Tracking) 더욱 세밀한 이분 격자법 다차원 극값을 구하기 위한 경사도 탐색법 가파른 경사법 해(근)을 구하는 방법 극값을 구하는 방법

3 We are now … 이차원 이분 격자(bisection grid)법 영점 곡선 추적 (Zero-Curve Tracking)
더욱 세밀한 이차원 이분 격자법 다차원 극값을 구하기 위한 경사도 탐색 (Gradient Search) 가파른 경사법 (Steepest Descent)

4 이변수 방정식의 의미 이변수 방정식 f(x,y)=0의 해는 x-y 평면의 궤적이다.
Bisection Grid 이변수 방정식 f(x,y)=0의 해는 x-y 평면의 궤적이다. 예를 들어, 형태인 선형 방정식의 해는 x-y 평면에서 다음과 같은 직선이 된다. 또한, 비선형 방정식 의 해는 타원의 궤적에 해당한다. 중심이 원점이고, x축 길이가 a, y축 길이가 b인 타원 방정식

5 Recall: 일차원(일변수) 이분법 구간 분할: 중간 값을 취하는 방법을 사용한다.
Bisection Grid 구간 분할: 중간 값을 취하는 방법을 사용한다. 두 값 xl 과 xh 사이에 근이 존재할 때, 중간 값 xm은 다음과 같이 구한다. f(x) Xl Xh Xm Xl’ Xh’ Xm’ x f(xm)f(xh)와 f(xm)f(xl)을 조사하여 음수 값을 갖는 경우를 다음 구간으로 사용한다.

6 이차원 이분 격자법 개념 (1/2) Bisection Grid 일차원 이분법을 이차원으로 확장한 방법이다. 1) 일정한 크기의 격자로 나누고, 2) 해당 격자에서 x축 및 y축에 대해 이분법을 적용하여 범위를 축소시키면서 에러 범위 내의 (x, y) 해를 찾는다. (일차원) 이분법 이차원 이분 격자법

7 이차원 이분 격자법 개념 (2/2) Bisection Grid 격자 내의 y축 검사 격자 내의 x축 검사

8 Recall: 일차원 이분법 알고리즘 procedure bisection(xl, xh, e: real numbers)
Bisection Grid procedure bisection(xl, xh, e: real numbers) { xl is a left bound value of the range having a root.} { xh is a right bound value of the range having a root.} { e is an allowable error value.} while (xh − xl) > e begin xm := (xh + xl) / 2; {get a medium value} if f(xm)f(xh) = 0 then return xm; else if f(xm)f(xl) < 0 then xh := xm; else if f(xm)f(xh) < 0 then xl := xm; else break; { something wrong  cannot find the root.} end return xm;

9 이차원 이분 격자법 알고리즘 (1/3) Bisection Grid procedure bisection-grid(xl, xh, yl, yh, s, e: real numbers) { [xl, xh] is a domain of x.} { [yl, yh] is a domain of y.} { s is a sliding factor (or an interval factor) of a grid.} { e is an allowable error value.} x := xl; while (x  xh) begin y := yl; while (y  yh) bisx(x, y, y+s, e); { find a root on x where y is in (y, y+s) } bisy(y, x, x+s, e); { find a root on y where x is in (x, x+s) } y := y + s; end x := x + s; (xh, yh) (xl, yl)

10 이차원 이분 격자법 알고리즘 (2/3) bisx(): 함수 f(x,y)에서 x 값을 상수로 보고, y에 대한 근을 찾는다.
Bisection Grid bisx(): 함수 f(x,y)에서 x 값을 상수로 보고, y에 대한 근을 찾는다. procedure bisx(x, yl, yh, e: real numbers) if f(x,yl)f(x,yh)  0 return; {no root, or cannot find the root} while (yh − yl) > e begin ym := (yh + yl) / 2; {get a medium value} if f(x,ym)f(x,yh) = 0 then break; else if f(x,ym)f(x,yl) < 0 then yh := ym; else if f(x,ym)f(x,yh) < 0 then yl := ym; else return; { something wrong  cannot find the root.} end Insert (x, ym) into the root set;

11 이차원 이분 격자법 알고리즘 (3/3) bisy(): 함수 f(x,y)에서 y 값을 상수로 보고, x에 대한 근을 찾는다.
Bisection Grid bisy(): 함수 f(x,y)에서 y 값을 상수로 보고, x에 대한 근을 찾는다. procedure bisy(y, xl, xh, e: real numbers) if f(xl,y)f(xh,y)  0 return; {no root, or cannot find the root} while (xh − xl) > e begin xm := (xh + xl) / 2; {get a medium value} if f(xm,y)f(xh,y) = 0 then break; else if f(xm,y)f(xl,y) < 0 then xh := xm; else if f(xm,y)f(xh,y) < 0 then xl := xm; else return; { something wrong  cannot find the root.} end Insert (xm, y) into the root set;

12 이차원 이분 격자법 프로그램 (1/4) 대상 함수: #include <stdio.h>
Bisection Grid 대상 함수: #include <stdio.h> #include <stdlib.h> #include <math.h> float f(float, float); void bisx(float, float, float, float); void bisy(float, float, float, float); main(int argc, char *argv[]) { int i = 1; float x, xl, xh, y, yl, yh, s, e; if(argc < 7) { printf("Usage: %s xl xh yl yh s e\n", argv[0]); exit(0); } xl = (float)atof(argv[1]); xh = (float)atof(argv[2]); yl = (float)atof(argv[3]); yh = (float)atof(argv[4]); s = (float)atof(argv[5]); e = (float)atof(argv[6]);

13 이차원 이분 격자법 프로그램 (2/4) printf("(xl, xh) = (%.8f, %.8f)\n", xl, xh);
Bisection Grid printf("(xl, xh) = (%.8f, %.8f)\n", xl, xh); printf("(yl, yh) = (%.8f, %.8f)\n", yl, yh); printf("s = %.8f\n", s); printf("e = %.8f\n", e); printf(" x\t\t y\t\t f(x,y)\n"); for(x = xl;x <= xh;x += s) { for(y = yl;y <= yh;y += s) { bisx(x, y, y+s, e); bisy(y, x, x+s, e); } float f(float x, float y) { return ( 3.0*sin(3.0*x) + 4.0*cos(3.0*y) );

14 이차원 이분 격자법 프로그램 (3/4) void bisx(float x, float yl, float yh, float e)
Bisection Grid void bisx(float x, float yl, float yh, float e) { float ym; if((f(x,yh)*f(x,yl)) >= 0) return; while((yh-yl) > e) { ym = (yh + yl) / 2.0; if((f(x,ym)*f(x,yh)) == (float)0) break; else if((f(x,ym)*f(x,yl)) < (float)0) yh = ym; else if((f(x,ym)*f(x,yh)) < (float)0) yl = ym; else return; } printf("%.8f\t%.8f\t%.8f\n", x, ym, f(x,ym));

15 이차원 이분 격자법 프로그램 (4/4) void bisy(float y, float xl, float xh, float e)
Bisection Grid void bisy(float y, float xl, float xh, float e) { float xm; if((f(xh,y)*f(xl,y)) >= 0) return; while((xh-xl) > e) { xm = (xh + xl) / 2.0; if((f(xm,y)*f(xh,y)) == (float)0) break; else if((f(xm,y)*f(xl,y)) < (float)0) xh = xm; else if((f(xm,y)*f(xh,y)) < (float)0) xl = xm; else return; } printf("%.8f\t%.8f\t%.8f\n", xm, y, f(xm,y));

16 프로그램 실행 결과 (1/2) Bisection Grid

17 프로그램 실행 결과 (2/2) Bisection Grid

18 We are now … 이차원 이분 격자(bisection grid)법 영점 곡선 추적 (Zero-Curve Tracking)
더욱 세밀한 이차원 이분 격자법 다차원 극값을 구하기 위한 경사도 탐색 (Gradient Search) 가파른 경사법 (Steepest Descent)

19 영점-곡선 추적의 동기(motivation)
Zero-Curve Tracking 이분 격자법은 Domain 내의 모든 구간에 대해서 해를 구하는 시도를 해야 하므로, 불필요한 공간 탐색이 많이 이루어진다. 영점-곡선 추척에서는 1) (이분 격자법 등을 사용하여) 한 점(정확히는 두 점)을 먼저 찾아낸 후, 2) 찾아낸 점을 사용하여 다음 점을 찾아내는 방법을 사용한다. 이분 격자법에 비해서 검색 공간을 줄일 수 있다는 장점이 있다. But, 계산 과정이 비교적 복잡한 단점이 있다.

20 영점-곡선 추적법의 개념 1) (이분 격자법 등을 사용하여) 첫 번째 점 (x1, y1)을 찾아낸다.
Zero-Curve Tracking 1) (이분 격자법 등을 사용하여) 첫 번째 점 (x1, y1)을 찾아낸다. 2) 첫 번째 점 (x1, y1)에서 y축으로 w만큼 떨어진 곳에서, (이분 격자법 등을 사용하여) 두 번째 점 (x2, y2)를 찾아낸다. 3) 두 점 (x1, y1)과 (x2, y2)를 연결한 직선 상에서, (x2, y2)와 직교하는 직선을 구하고, 이를 w만큼 평행 이동한 직선과 곡선이 만나는 점을 세 번째 점 (x3, y3)로 삼는다. 4) 두 번째 및 세 번째 점을 사용하여 상기 2) ~ 3)의 과정을 반복한다. w (x3,y3) w (x2,y2) (x1,y1)

21 영점-곡선 추적에서 세 번째 점 계산 (1/6) 첫 번째 점을 (u,v), 두 번째 점을 (x,y)라 하자.
Zero-Curve Tracking 첫 번째 점을 (u,v), 두 번째 점을 (x,y)라 하자. 다음 그림은 이 두 점을 잇는 직선과 수직인 탐색선과의 관계를 나타낸다. 이때, 탐색선의 길이는 2w라 하고, 탐색선의 양 끝점을 각각 (a,s)와 (b,t)라 하자. (a,s) w w (x,y) w (b,t) (u,v)

22 영점-곡선 추적에서 세 번째 점 계산 (2/6) Zero-Curve Tracking 두 점 (u,v)와 (x,y), 그리고 두 점을 잇는 직선과 탐색선과의 교점이 이루는 삼각형에서 다음 관계가 성립한다. (a,s) w w (x,y) w d  (b,t) (u,v)

23 영점-곡선 추적에서 세 번째 점 계산 (3/6) 교점의 좌표 (i, j)는 다음과 같이 구할 수 있다. (a,s) w  w
Zero-Curve Tracking 교점의 좌표 (i, j)는 다음과 같이 구할 수 있다. (a,s) w  w (i,j) (x,y)   d w  (b,t) (u,v)

24 영점-곡선 추적에서 세 번째 점 계산 (4/6) 탐색선의 시작 및 끝 좌표는 다음과 같이 구할 수 있다. (a,s) w  w
Zero-Curve Tracking 탐색선의 시작 및 끝 좌표는 다음과 같이 구할 수 있다. (a,s) w  (i,j)  w (x,y) w    d  (b,t) (u,v)

25 영점-곡선 추적에서 세 번째 점 계산 (5/6) Zero-Curve Tracking 곡선과의 교점을 구하기 위하여, 탐색선의 시작점 (a, s)에서 c 만큼씩 이동하면서 이동 전의 점과 이동 후의 점에 대한 함수 값의 부호가 변화하는지 확인한다. (부호가 변하면, 그 중간에 해가 존재하기 때문이다.) (a,s) c  (g,h) c  이 두 점을 대상으로 이분법을 수행하여 원하는 교점을 찾아낸다. c 

26 영점-곡선 추적에서 세 번째 점 계산 (6/6) Zero-Curve Tracking 부호가 변하는 두 점을 찾았으면, c 를 절반으로 줄인 후, 다시금 이동하면서 이동 전의 점과 이동 후의 점에 대한 함수 값의 부호가 변화하는지 확인한다. (부호가 변하면, 그 중간에 해가 존재하기 때문이다.) c  c  c 원하는 정확도의 교점을 찾을 때까지 상기 과정을 반복한다.

27 영점-곡선 추적법 알고리즘 (1/2) Zero-Curve Tracking procedure zero-curve-tracking(xi, yi, w, c, e: real numbers) { xi and yi are starting points.} { w is a distance factor.} { c is a an interval factor.} { e is an allowable error value.} (u, v) = root(xi, yi, c, 0, e); // find the first root (x, y) = root(xi, yi+w, c, 0, e); // find the second root while (true) begin d := ; A := (x – u) / d; B := (y – v) / d; // A = cos , B = sin  u := x; v := y; x := x + w(A-B); y := y + w(A+B); (x, y) = root(x, y, cB, -cA, e); // p = cB, q = -cA end

28 영점-곡선 추적법 알고리즘 (2/2) procedure root(x, y, p, q, e: real numbers)
Zero-Curve Tracking procedure root(x, y, p, q, e: real numbers) while (true) begin xn := x + p; yn := y + q; if f(x,y)f(xn,yn)  0 then if (|p| > e) || (|q| > e) then p := p/2; q := q/2; end else return (xn, yn); else begin x := xn; y := yn; end

29 영점-곡선 추적법 프로그램 (1/4) Zero-Curve Tracking 대상 함수:

30 영점-곡선 추적법 프로그램 (2/4) Zero-Curve Tracking

31 영점-곡선 추적법 프로그램 (3/4) Zero-Curve Tracking

32 영점-곡선 추적법 프로그램 (4/4) Zero-Curve Tracking

33 영점-곡선 추적법 프로그램 실행 결과 Zero-Curve Tracking

34 영점-곡선 추적법 다른 예제 (1/2) Zero-Curve Tracking 대상 함수:

35 영점-곡선 추적법 다른 예제 (2/2) Zero-Curve Tracking 대상 함수:

36 We are now … 이차원 이분 격자(bisection grid)법 영점 곡선 추적 (Zero-Curve Tracking)
Zero-Crossing Squares 이차원 이분 격자(bisection grid)법 영점 곡선 추적 (Zero-Curve Tracking) 더욱 세밀한 이차원 이분 격자법 (영점-교차 사각형) 다차원 극값을 구하기 위한 경사도 탐색 (Gradient Search) 가파른 경사법 (Steepest Descent)

37 영점-교차 사각형의 동기(motivation) (1/2)
Zero-Crossing Squares 이분 격자법은 Domain 내의 많은 구간에 대해서 해를 구하는 시도를 해야 하므로, 불필요한 공간 탐색이 많이 이루어진다. 또한, 탐색 공간을 줄이기 위하여, 구간을 넓게 할 경우 정확한 해를 찾기가 어렵다. 영점-교차 사각형 방법에서는 1) 비교적 큰 사각형으로 구간을 분할한 후, 2) (이분 격자법 등을 사용하여) 해당 사각형들이 영점을 포함하는지 여부를 확인하여, 3) 영점을 포함하는 사각형에 대해서는 보다 세밀한 사각형을 구성하여 영점을 확인하는 방법을 사용한다. In CS, we call this technique as “filtering & refining.”

38 영점-교차 사각형의 동기(motivation) (2/2)
Zero-Crossing Squares 이차원 이분 격자법 영점-교차 사각형 방법

39 영점-교차 사각형의 개념 (1/5) Zero-Crossing Squares (정확히 이야기하면) 근이 있는 구간을 알아낸 후, 그 구간에 대해서 더욱 자세한 근을 구할 때 사용한다. 구간 내의 직선과 곡선의 교점을 구하는 방식이 아니라, 곡선이 지나가는 구간 자체를 파악하는 방식을 사용한다. (구간의 좌하점(혹은 우상점)을 근사해로 사용한다.) 구간의 크기가 작아지면, 결과적으로 정밀한 해를 찾을 수 있기 때문이다.

40 영점-교차 사각형의 개념 (2/5) 곡선이 지나가는 구간을 파악하는 방법: 두 개의 이차원 배열을 사용한다.
Zero-Crossing Squares 곡선이 지나가는 구간을 파악하는 방법: 두 개의 이차원 배열을 사용한다. ai,j: 구간을 나누는 직선의 교점이 가지는 함수 값을 나타낸다. li,j: 직선이 해당 구간을 지나는지의 여부를 나타낸다. a1,5 a2,5 a3,5 a4,5 a5,5 l1,4 l2,4 l3,4 l4,4 a1,4 a2,4 a3,4 a4,4 a5,4 l1,3 l2,3 l3,3 l4,3 a1,3 a2,3 a3,3 a4,3 a5,3 l1,2 l2,2 l3,2 l4,2 a1,2 a2,2 a3,2 a4,2 a5,2 l1,1 l2,1 l3,1 l4,1 a1,1 a2,1 a3,1 a4,1 a5,1

41 영점-교차 사각형의 개념 (3/5) 이차원 배열 ai,j의 구성 방법
Zero-Crossing Squares 이차원 배열 ai,j의 구성 방법 ai,j는 다음과 같이 구해지는 (x, y) 좌표 값의 함수 값을 가진다. 결국, ai,j는 시작 점 (u, v)에서 x 축으로 (i-1), y 축으로 (j-1)을 c 만큼씩 이동한 점의 함수 값이다.

42 영점-교차 사각형의 개념 (4/5) 이차원 배열 li,j의 구성 방법 (1/2) (Note: 모든 li,j의 초기값은 0임)
Zero-Crossing Squares 이차원 배열 li,j의 구성 방법 (1/2) (Note: 모든 li,j의 초기값은 0임) ai,j = 0 인 경우: 주변의 네 구역 모두 해로 포함시킨다. ai-1,j-1 ai-1,j ai-1,j+1 ai,j-1 ai,j ai,j+1 ai+1,j-1 ai+1,j ai+1,j+1 li-1,j-1=1 li-1,j=1 li,j-1=1 li,j=1

43 영점-교차 사각형의 개념 (5/5) 이차원 배열 li,j의 구성 방법 (2/2)
Zero-Crossing Squares 이차원 배열 li,j의 구성 방법 (2/2) (i,j)와 (i+1,j)에 교점이 있는 경우 (즉, ai,j ai+1,j < 0인 경우) : 상하 두 구간을 해에 포함시킨다. (i,j)와 (i,j+1)에 교점이 있는 경우 (즉, ai,j ai,j+1 < 0인 경우) : 좌우 두 구간을 해에 포함시킨다. li,j-1=1 li,j=1 ai,j-1 ai,j ai,j+1 ai+1,j-1 ai+1,j ai+1,j+1 ai-1,j ai-1,j+1 ai,j ai,j+1 ai+1,j ai+1,j+1 li-1,j=1 li,j=1

44 영점-교차 사각형 알고리즘 (1/2) Zero-Crossing Squares procedure zero-crossing-square(x, y, g, d: real numbers) { x and y are starting points.} { g is the number of grid divisions (한 축에 대한 division 수)} { d is the dimension (한 축의 길이)} c := d / g; for each i for each j begin ai,j := f(x + c(i-1), y + c(j-1)); li,j := 0; end (1) Calculate f(x,y) and assign it to ai,j. (2) Initialize li,j to 0.

45 영점-교차 사각형 알고리즘 (2/2) for each i for each j begin if ai,j = 0 then
Zero-Crossing Squares for each i for each j begin if ai,j = 0 then li,j := li-1,j := li,j-1 := li-1,j-1 := 1; else if ai,j  ai+1,j < 0 then li,j := li,j-1 := 1; else if ai,j  ai,j+1 < 0 then li,j := li-1,j := 1; end return (x+c(i-1), y+c(j-1)) as a root if li,j = 1; Calculate li,j by using ai,j.

46 영점-교차 사각형 프로그램 (1/2) Zero-Crossing Squares 대상 함수:

47 영점-교차 사각형 프로그램 (2/2) Zero-Crossing Squares

48 영점-교차 사각형 실행 결과 (1/2) Zero-Crossing Squares

49 영점-교차 사각형 실행 결과 (2/2) Zero-Crossing Squares