bird function의 해 구하기
오랜만에 스터디를 하다가 스터디 과제 하나를 들고 왔다. 다른 과제나 진도 따라갈 수업도 산더미인데 언제하지 ㅅㅂ. 하여튼 과제는 일단 간단하다. bird function f(x,y) 에 대해 f(x,y)=0을 만족하는 해의 쌍(-2pi<=x<=2pi, -2pi<=y<=2pi)을 구하는 것. bird function은 무슨 최적화에 쓰이는 듣보잡 함수라고 한다. 암튼 내가 모르면 듣보잡임
bird function
Mazhar Ansari Ardeh
그럼 일단 매우 단순하게, 저 주어진 범위를 엄청나게 많은 격자로 나눠서 해에 가까운 점을 찾아보는 건 어떨까? 물론 격자점들 사이에 있는 해는 잡아낼 수 없고, 실수 오차 때문에 정확한 해를 더욱더 찾을 수 없겠지만 일단 해에 가까운 점들이 어느 정도 있는지를 대충 예상해 볼 수 있을 것 같다.
#include <stdio.h>
#include <math.h>
double e = 2.71828;
double pi = 3.14159;
double bird(double x, double y) {
return (sin(x) * pow(e, pow(1 - cos(y), 2)) + cos(y) * pow(e, pow(1 - sin(x), 2)) + pow(x - y, 2));
}
int main() {
int cnt = 0;
for (double nx = -2 * pi; nx < 2 * pi; nx += 4 * pi / 10000) {
for (double ny = -2 * pi; ny < 2 * pi; ny += 4 * pi / 10000) {
if (-1e-5 < bird(nx, ny) && bird(nx, ny) < 1e-5) {
printf("x = %lf, y = %lf\n", nx, ny);
printf("bird function value : %lf\n", bird(nx, ny));
cnt++;
}
}
}
printf("number of numerical root : %d", cnt);
}
주어진 범위를 10000x10000, 즉 1억개의 격자로 나누고 격자점들 중에 어느 정도 0에 가까운 함숫값을 뱉는 점들을 찾아봤다.
총 7개 나온다. 참고로 값의 범위를 1e-5가 아니라 1e-4로 좀 더 널널하게 잡아줬을 때는 135개의 점이 나왔다.
그럼 여기서 더 정확한 해를 찾으려면 어떻게 할까? 물론 주어진 범위를 10만x10만이라든지, 1억x1억같이 더 촘촘한 칸으로 나누면 더 정확한 해를 구할 수 있을 것이다. 하지만 우리는 모든 점의 함숫값을 알고 싶은 것이 아니고 그냥 f(x,y)=0이 되는 (x,y)를 알고 싶을 뿐이다. 그러면 모든 점을 더 세부적으로 나눌 것이 아니라, 가능성이 있는 점들에서만 더 세분화해서 탐색하면 된다.
직관적으로, 만약 A점과 B점 사이에서 함숫값의 부호가 바뀐다면 그 근처 어딘가에 함숫값이 0이 되는 지점이 반드시 있다. 그리고 A점과 B점 사이의 거리가 매우 작다면 대충 그 중간 정도 값을 함숫값이 0이 되는 지점이라고 생각할 수 있을 것이다.
#include <stdio.h>
#include <math.h>
double e = 2.71828;
double pi = 3.14159;
double bird(double x, double y) {
return (sin(x) * pow(e, pow(1 - cos(y), 2)) + cos(y) * pow(e, pow(1 - sin(x), 2)) + pow(x - y, 2));
}
int main() {
int cnt = 0;
double dx = 4 * pi / 10000;
double dy = dx;
for (double nx = -2 * pi; nx < 2 * pi; nx += dx) {
for (double ny = -2 * pi; ny < 2 * pi; ny += dy) {
if (bird(nx, ny) * bird(nx + dx, ny + dy) < 0) {
printf("x = %lf, y = %lf\n", nx + dx / 2, ny + dy / 2);
printf("bird function value : %lf\n", bird(nx + dx / 2, ny + dy / 2));
cnt++;
}
}
}
printf("number of numerical root : %d", cnt);
}
아까와 같은 방식으로 함숫값의 부호가 달라지는 격자들을 찾았더니 무려 24158개가 나온다. 대부분 0.1 미만의 오차를 가지는 값들이지만 그래도 훨씬 더 세분화할 필요가 있다. 이제 이분탐색의 아이디어를 이용하자. 탐색하는 구간을 점점 좁혀나가면서 가장 오차가 작은 곳을 찾는 것이다.
#include <stdio.h>
#include <math.h>
double e = 2.718281828459045;
double pi = 3.141592653589793;
double bird(double x, double y) {
return (sin(x) * pow(e, pow(1 - cos(y), 2)) + cos(y) * pow(e, pow(1 - sin(x), 2)) + pow(x - y, 2));
}
int main() {
int cnt = 0;
//해의 개수를 세줌
double dx = 4 * pi / 1000;
double dy = dx;
for (double nx = -2 * pi; nx < 2 * pi; nx += dx) {
for (double ny = -2 * pi; ny < 2 * pi; ny += dy) {
double lowx = nx, lowy = ny;
double highx = nx + dx, highy = ny + dy;
int iter = 0;
//몇 번의 이분탐색을 통해 해를 찾아냈는지를 세줌
while (lowx <= highx && lowy <= highy) {
double midx = (lowx + highx) / 2, midy = (lowy + highy) / 2;
if (bird(lowx, lowy) * bird(highx, highy) > 0) { break; }
//구간에서 함숫값의 부호가 달라지지 않는다고 추측될 때에는
//이분탐색을 진행할 필요조차 없음
printf("x = %.10lf, y = %.10lf\n", midx, midy);
printf("current function value : %.15lf\n", bird(midx, midy));
//이분탐색 과정을 출력해줌
if (-1e-13 < bird(midx, midy) && bird(midx, midy) < 1e-13) {
//오차가 1e-13 미만인 해를 찾았을 때
printf("we got the root, x = %.10lf, y = %.10lf\n", midx, midy);
printf("bird function value : %.15lf\n", bird(midx, midy));
printf("number of iteration : %d\n\n", iter);
cnt++;
break;
}
iter++;
if (bird(lowx, lowy) * bird(midx, midy) < 0) { highx = midx; highy = midy; }
if (bird(midx, midy) * bird(highx, highy) < 0) { lowx = midx; lowy = midy; }
}
}
}
printf("number of numerical root : %d\n", cnt);
}
위와 같은 코드를 보면 일단 해가 있을 법한 구간을 찾고 이진탐색을 이용해서 '해가 있을 만한 구간'을 좁혀나간다. 이 코드를 사용하면 평균적으로 40~45번의 탐색으로 오차가 1e-13 미만인 해를 찾을 수 있다. 약 2400개의 근사해를 순식간에 찾아준다. 와!
결국 이진탐색의 아이디어로 적당히 조져서 해 찾는다는 소리다. 사실 구간 나누면서 찾는다는 것만 이진탐색의 아이디어지 뭐 정렬된 구간이라든지 그런 건 없다. bisection method를 찾아보라..