수치해석 () 다변수방정식과함수 (Part 1) (Multi-Variable Equations and Functions Part 1) 2005 년가을학기 문양세컴퓨터과학과강원대학교자연과학대학 In this chapter 다변수방정식과함수 변수가두개이상인함수, 예를들어, f ( x, y, z) = log( x+ y) + sin( x+ z) 의해 (f(x,y,z)=0 으로하는 (x,y,z) 의값 ) 를수치해석적으로구하는방법을다룬다. 2 차함수인 f(x,y) 를집중적으로다룬다. (3 차이상확장가능 ) We ill cover 이차원이분격자법 () 영점곡선추적 () 더욱세밀한이분격자법 다차원극값을구하기위한경사도탐색법 가파른경사법 해 ( 근 ) 을구하는방법 극값을구하는방법 Page 2 1
We are no 이차원이분격자 (bisection grid) 법영점곡선추적 () 더욱세밀한이차원이분격자법다차원극값을구하기위한경사도탐색 (Gradient Search) 가파른경사법 (Steepest Descent) Page 3 이변수방정식의의미 이변수방정식 f(x,y)=0 의해는 x-y 평면의궤적이다. 예를들어, f ( x, y) = ax + by + c 형태인선형방정식의해는x-y 평면에서다음과같은직선이된다. a c y = x b b 또한, 비선형방정식한다. 2 2 y f( x, y) = x + 1 의해는타원의궤적에해당 4 Page 4 2
Recall: 일차원 ( 일변수 ) 이분법 구간분할 : 중간값을취하는방법을사용한다. 두값 x l 과 x h 사이에근이존재할때, 중간값 x m 은다음과같이구한다. f(x) x m = x l + x 2 h X l X m X h x X l X m X h f(x m ) f(x h ) 와 f(x m ) f(x l ) 을조사하여음수값을갖는경우를다음구간으로사용한다. Page 5 이차원이분격자법개념 (1/2) 일차원이분법을이차원으로확장한방법이다. 1) 일정한크기의격자로나누고, 2) 해당격자에서 x 축및 y 축에대해이분법을적용하여범위를축소시키면서에러범위내의 (x, y) 해를찾는다. ( 일차원 ) 이분법 이차원이분격자법 Page 6 3
이차원이분격자법개념 (2/2) 격자내의 y 축검사 격자내의 x 축검사 Page 7 Recall: 일차원이분법알고리즘 procedure bisection(x l, x h, e: real numbers) { x l is a left bound value of the range having a root.} { x h is a left bound value of the range having a root.} { e is an alloable error value.} hile (x h x l ) > e begin x m := (x h + x l ) / 2; {get a medium value} if f(x m ) f(x h ) = 0 then return x m ; else if f(x m ) f(x l ) < 0 then x h := x m ; else if f(x m ) f(x h ) < 0 then x l := x m ; else break; { something rong cannot find the root.} end return x m ; Page 8 4
Notice ( 있지도않지만 ) 교재의알고리즘을참고하지마세요. 교재의프로그램은더더욱참고하지마세요. Computer Scientist에의해작성되지않아서, 알고리즘 ( 프로그램 ) 이그다지 Clear하지않습니다. 강의노트의알고리즘과프로그램을참고하기바랍니다. Page 9 이차원이분격자법알고리즘 (1/3) procedure bisection-grid(x l, x h, y l, y h, s, e: real numbers) { [x l, x h ] is a domain of x.} { [y l, y h ] is a domain of y.} { s is a sliding factor (or an interval factor) of a grid.} { e is an alloable error value.} x := x l ; hile (x x h ) begin y := y l ; hile (y y h ) begin bisx(x, y, y+s, e); { find a root on x here y is in (y, y+s) } bisy(y, x, x+s, e); { find a root on y here x is in (x, x+s) } y := y + s; end x := x + s; end Page 10 5
이차원이분격자법알고리즘 (2/3) procedure bisx(x, y l, y h, e: real numbers) if f(x,y l ) f(x,y h ) 0 return; {no root, or cannot find the root} hile (y h y l ) > e begin y m := (y h + y l ) / 2; {get a medium value} if f(x,y m ) f(x,y h ) = 0 then break; else if f(x,y m ) f(x,y l ) < 0 then y h := y m ; else if f(x,y m ) f(x,y h ) < 0 then y l := y m ; else return; { something rong cannot find the root.} end Insert (x, y m ) into the root set; Page 11 이차원이분격자법알고리즘 (3/3) procedure bisy(y, x l, x h, e: real numbers) if f(x l,y) f(x h,y) 0 return; {no root, or cannot find the root} hile (x h x l ) > e begin x m := (x h + x l ) / 2; {get a medium value} if f(x m,y) f(x h,y) = 0 then break; else if f(x m,y) f(x l,y) < 0 then x h := x m ; else if f(x m,y) f(x h,y) < 0 then x l := x m ; else return; { something rong cannot find the root.} end Insert (x m, y) into the root set; Page 12 6
이차원이분격자법프로그램 (1/4) 대상함수 : #include <stdio.h> #include <stdlib.h> #include <math.h> f ( x, y) = 3 sin(3 x) + 4 cos(3 y) 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]); Page 13 이차원이분격자법프로그램 (2/4) 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) ); } Page 14 7
이차원이분격자법프로그램 (3/4) void bisx(float x, float yl, float yh, float e) { float ym; if((f(x,yh)*f(x,yl)) >= 0) return; } hile((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)); Page 15 이차원이분격자법프로그램 (4/4) void bisy(float y, float xl, float xh, float e) { float xm; if((f(xh,y)*f(xl,y)) >= 0) return; } hile((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)); Page 16 8
프로그램실행결과 (1/2) 1.00 0.90 0.80 0.70 0.60 0.50 0.40 0.30 0.20 0.10 0.00 0.00 0.20 0.40 0.60 0.80 1.00 Page 17 프로그램실행결과 (2/2) 0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00 0.00 0.20 0.40 0.60 0.80 1.00 Page 18 9
We are no 이차원이분격자 (bisection grid) 법영점곡선추적 () 더욱세밀한이차원이분격자법다차원극값을구하기위한경사도탐색 (Gradient Search) 가파른경사법 (Steepest Descent) Page 19 영점 -곡선추적의동기 (motivation) 이분격자법은 Domain 내의모든구간에대해서해를구하는시도를해야하므로, 불필요한공간탐색이많이이루어진다. 영점 - 곡선추척에서는 1) ( 이분격자법등을사용하여 ) 한점 ( 정확히는두점 ) 을먼저찾아낸후, 2) 찾아낸점을사용하여다음점을찾아내는방법을사용한다. 이분격자법에비해서검색공간을줄일수있다는장점이있다. But, 계산과정이비교적복잡한단점이있다. Page 20 10
영점 -곡선추적법의개념 1) ( 이분격자법등을사용하여 ) 첫번째점 (x 1, y 1 ) 을찾아낸다. 2) 첫번째점 (x 1, y 1 ) 에서 y 축으로 만큼떨어진곳에서, ( 이분격자법등을사용하여 ) 두번째점 (x 2, y 2 ) 를찾아낸다. 3) 두점 (x 1, y 1 ) 및 (x 2, y 2 ) 를연결한직선상에서, (x 2, y 2 ) 와직교하는직선을구하고, 이직선과곡선이만나는점을세번째점 (x 3, y 3 ) 로삼는다. 4) 두번째및세번째점을사용하여상기 2) ~ 3) 의과정을반복한다. (x 3,y 3 ) (x 2,y 2 ) (x 1,y 1 ) Page 21 영점 -곡선추적에서세번째점계산 (1/6) 첫번째점을 (u,v), 두번째점을 (x,y) 라하자. 다음그림은이두점을잇는직선과수직인탐색선과의관계를나타낸다. 이때, 탐색선의길이는 2 라하고, 탐색선의양끝점을각각 (a,s) 와 (b,t) 라하자. (a,s) (x,y) (u,v) (b,t) Page 22 11
영점 -곡선추적에서세번째점계산 (2/6) 두점 (u,v) 와 (x,y), 그리고두점을잇는직선과탐색선과의교점이이루는삼각형에서다음관계가성립한다. 2 2 x u y v d= ( x u) + ( y v), cos θ=, sinθ= d d (a,s) (x,y) d (u,v) θ (b,t) Page 23 영점 -곡선추적에서세번째점계산 (3/6) 교점의좌표 (i, j) 는다음과같이구할수있다. ( i, j) = ( x+ cos θ, y+ sinθ) (a,s) (u,v) d θ (x,y) θ (i,j) (b,t) α β =α sin θ θ γ =α cosθ Page 24 12
영점 -곡선추적에서세번째점계산 (4/6) 탐색선의시작및끝좌표는다음과같이구할수있다. ( as, ) = ( x+ cosθ sin θ, y+ sin θ+ cosθ) ( bt, ) = ( x+ cosθ+ sin θ, y+ sin θ cosθ) (a,s) θ (u,v) d θ (x,y) θ (i,j) θ (b,t) α β =α sin θ θ γ =α cosθ Page 25 영점 -곡선추적에서세번째점계산 (5/6) 곡선과의교점을구하기위하여, 탐색선의시작점 (a, s) 에서 c 만큼씩이동하면서이동전의점과이동후의점에대한함수값의부호가변화하는지확인한다. ( 부호가변하면, 그중간에해가존재하기때문이다.) p= c sin θ q = c cosθ (a,s) θ c (i,j) θ c θ c ( i, j) = ( a+ c sin θ, s c cosθ) 이두점을대상으로이분법을수행하여원하는교점을찾아낸다. Page 26 13
영점 -곡선추적에서세번째점계산 (6/6) 부호가변하는두점을찾았으면, c 를절반으로줄인후, 다시금이동하면서이동전의점과이동후의점에대한함수값의부호가변화하는지확인한다. ( 부호가변하면, 그중간에해가존재하기때문이다.) θc c θ c 원하는정확도의교점을찾을때까지상기과정을반복한다. Page 27 영점 -곡선추적법알고리즘 (1/2) procedure zero-curve-tracking(x i, y i,, c, e: real numbers) { x i and y i are starting points.} { is a distance factor.} { c is a an interval factor.} { e is an alloable error value.} (u, v) = root(x i, y i, c, 0, e); (x, y) = root(x i, y i +, c, 0, e); hile (true) begin 2 2 d := ( x u) + ( y v) ; A := (x u) / d; B := (y v) / d; u := x; v := y; x := x + (A-B); y := y + (A+B); (x, y) = root(x, y, c B, -c A, e); end Page 28 14
영점 -곡선추적법알고리즘 (2/2) procedure root(x, y, p, q, e: real numbers) hile (true) begin x n := x + p; y n := y + q; if f(x,y) f(x n,y n ) 0 then begin if (p > e) (q > e) then begin p := p/2; q := q/2; end else return (x n, y n ); end else begin x := x n ; y := y n ; end end Page 29 영점 -곡선추적법프로그램 (1/4) 2 2 대상함수 : f( x, y) = x + y 4 Page 30 15
영점 -곡선추적법프로그램 (2/4) Page 31 영점 -곡선추적법프로그램 (3/4) Page 32 16
영점 -곡선추적법프로그램 (4/4) Page 33 영점 -곡선추적법프로그램실행결과 2.50 2.00 1.50 1.00 0.50 0.00-2.50-1.50-0.50 0.50 1.50 2.50-0.50-1.00-1.50-2.00-2.50 Page 34 17
We are no 이차원이분격자 (bisection grid) 법영점곡선추적 () 더욱세밀한이차원이분격자법 ( 영점-교차사각형 ) 다차원극값을구하기위한경사도탐색 (Gradient Search) 가파른경사법 (Steepest Descent) Page 35 영점 -교차사각형의동기 (motivation) (1/2) 이분격자법은 Domain 내의많은구간에대해서해를구하는시도를해야하므로, 불필요한공간탐색이많이이루어진다. 또한, 탐색공간을줄이기위하여, 구간이넓게할경우정확한해를찾기가어렵다. 영점 - 교차사각형방법에서는 1) 비교적큰사각형으로구간을분할한후, 2) ( 이분격자법등을사용하여 ) 해당사각형들이영점을포함하는지여부를확인하여, 3) 영점을포함하는사각형에대해서는보다세밀한사각형을구성하여영점을확인하는방법을사용한다. In CS, e call this technique as filtering & refining. Page 36 18
영점 -교차사각형의동기 (motivation) (2/2) 이차원이분격자법 영점 - 교차사각형방법 Page 37 영점 -교차사각형의개념 (1/5) 구간내의직선과곡선의교점을구하는방식이아니라, 곡선이지나가는구간자체를파악하는방식을사용한다. ( 구간의좌하점 ( 혹은우상점 ) 을근사해로사용한다.) 구간의크기가작아지면, 결과적으로정밀한해를찾을수있기때문이다. Page 38 19
영점 -교차사각형의개념 (2/5) 곡선이지나가는구간을파악하는방법 : 두개의이차원배열을사용한다. a i,j : 구간을나누는직선의교점의함수값을나타낸다. l i,j : 직선이해당구간을지나는지의여부를나타낸다. a 1,5 a 2,5 a 3,5 a 4,5 a 5,5 l 1,4 l 2,4 l 3,4 l 4,4 a 1,5 a 2,4 a 3,4 a 4,4 a 5,4 l 1,3 l 2,3 l 3,3 l 4,3 a 1,3 a 2,3 a 3,3 a 4,3 a 5,3 l 1,2 l 2,2 l 3,2 l 4,2 a 1,2 a 2,2 a 3,2 a 4,2 a 5,2 l 1,1 l 2,1 l 3,1 l 4,1 a 1,1 a 2,1 a 3,1 a 4,1 a 5,1 Page 39 영점 -교차사각형의개념 (3/5) 이차원배열 a i,j 의구성방법 a i,j 는다음과같이구해지는 (x, y) 좌표값의함수값을가진다. ai, j= f( x, y) x= u+ ( i 1) c y = v + ( j 1) c 결국, a i,j 는시작점 (u, v) 에서 x 축으로 (i-1), y 축으로 (j-1) 을 c 만큼씩이동한좌표이다. Page 40 20
영점 -교차사각형의개념 (4/5) 이차원배열 l i,j 의구성방법 (1/2) (Note: 모든 l i,j 의초기값은 0 임 ) a i,j = 0 인경우 : 주변의네구역모두해로포함시킨다. a i-1,j+1 a i,j+1 a i+1,j+1 l i-1,j =1 l i,j =1 a i-1,j a i,j a i+1,j l i-1,j-1 =1 l i,j-1 =1 a i-1,j-1 a i,j-1 a i+1,j-1 Page 41 영점 -교차사각형의개념 (5/5) 이차원배열 l i,j 의구성방법 (2/2) (i,j) 와 (i+1,j) 에교점이있는경우 ( 즉, a i,j a i+1,j < 0인경우 ) : 상하두구간을해에포함시킨다. (i,j) 와 (i,j+1) 에교점이있는경우 ( 즉, a i,j a i,j+1 < 0인경우 ) : 좌우두구간을해에포함시킨다. a i,j+1 a i+1,j+1 a i,j+1 l i,j =1 a i-1,j+1 a i+1,j+1 a i,j a i+1,j l i-1,j =1 l i,j =1 l i,j-1 =1 a i-1,j a i,j a i+1,j a i,j-1 a i+1,j-1 Page 42 21
영점 -교차사각형알고리즘 (1/2) procedure zero-crossing-square(x, y, g, d: real numbers) { x and y are starting points.} { g is the number of grid divisions.} { d is the dimension.} c := d / g; for each i for each j begin a i,j := f(x + c (i-1), y + c (i-1)); l i,j := 0; end Page 43 영점 -교차사각형알고리즘 (2/2) for each i for each j begin if a i,j = 0 then l i,j := l i-1,j := l i,j-1 := l i-1,j-1 := 0; else if a i,j a i+1,j < 0 then l i,j := l i,j-1 := 0; else if a i,j a i+1,j < 0 then l i,j := l i-1,j := 0; end return (x+c (i-1), y+c (i-1)) as a root if l i,j = 0; Page 44 22
영점 -교차사각형프로그램 (1/2) 2 ( ) ( ) 대상함수 : f( x, y) = x 1 + y 1 1 2 Page 45 영점 -교차사각형프로그램 (2/2) Page 46 23
영점 -교차사각형실행결과 (1/2) 2.50 2.00 1.50 1.00 0.50 0.50 1.00 1.50 2.00 2.50 Page 47 영점 -교차사각형실행결과 (2/2) 2.50 2.00 1.50 1.00 0.50 0.50 1.00 1.50 2.00 2.50 Page 48 24