Aug 21, 2013

Bisection method


Tình cờ *học lại* môn nhập môn lập trình bên MIT, thực ra đó là khóa Introduciton to Computer Science. Thuật toán này không có gì mới lạ, nhưng ứng dụng của nó khá là hay ho.

Bài toán: Tìm căn bậc 2 của 1 số thực dương cho trước.

Phân tích một cách thông thường, ta sẽ đi tìm nghiệm của phương trình: x^2- a = 0. Với a nhập từ bàn phím. Bài toán quay về tìm nghiệm của phương trình trên.

Có 2 phương pháp chính ở đây được giới thiệu là bisection method và Newton method, cả 2 phương pháp này khá là dễ hiểu. Và nếu tìm hiểu trên wiki, ta sẽ thấy 2 giải thuật trên là nền móng cho kha khá cả giải thuật khác.Link

Bisection method

Tư tưởng của Bisection khá là dễ hiểu: cho một hàm $f$ liên tục, nếu ta biết được răng $f(a)f(b) < 0$ thì đâu đó trong đoạn [a, b] có một giá trị x để f(x) bằng 0. Dĩ nhiên, mấu chốt vẫn là $f$ phải liên tục. Ta tiếp tục thu nhỏ khoảng $[a,b]$ này lại cho đến khi nào giá trị tìm được của a hoặc b cho $f(x) = 0 $ hoặc khoảng cách a và b nhỏ hơn 1 epsilon nào đó. Để làm điều đó, ta xét giá trị c: $ c = \frac{a+b}/2 $. Nếu f(c) = 0, thì c là số cần tìm, còn nếu $f(c) < 0 $ thì chắc chắn giá trị cần tìm ở đâu đó trong đoạn [c, b]. Ngược lại nếu $f(c) > 0$ thì giá trị x trong khoảng [a, c]. Cứ tiếp tục vậy.

Đoạn code thử nghiệm

1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
def bisectionMethod(x, epsilon):
    """ Find square root of x with epsilon
    """
    assert x > 0
    assert epsilon > 0
    low = 0
    high = max(x, 1)
    value = (low + high) / 2.0
    loop = 0
    while (abs(value**2 - x) > epsilon ):
        if ((value**2-x) < -epsilon): low = value
        elif ((value**2-x) > epsilon): high = value
        value = (low + high) / 2.0
        loop += 1;
    print "Bisection method: Value:", value, ". Loop: ", loop

Newton method

Tư tưởng của phương pháp này dựa vào công thức tính đạo hàm của 1 hàm số:

\[ f'(x_1) = \frac{df(x_1)}{dx} \approx \frac{f(x_1) - f(x_2)}{x_1-x_2} \]

\[ x_2 = x_1 - \frac{f(x_1)}{f'(x_1)} - \frac{f(x_2)}{f'(x_1)} \]

Vì $f(x_2) = 0 $ nên :

\[ x_2 = x_1 - \frac{f(x_1)}{f'(x_1)} \]

Tổng quát hóa:

\[ x_n = x_{n-1} - \frac{f(x_{n-1})}{f'(x_{n-1})} \]

Đoạn code thử nghiệm:

1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
def newtonMethod(x, epsilon):
    """ Newton method
    """
    assert x > 0
    assert epsilon > 0
    value = x
    loop = 0
    while (abs(value**2 - x) > epsilon):
        value = value - (value**2-x)/(2.0*value)
        loop += 1
    print "Newton method. Value: ", value, ".Loop: ", loop

So sánh 2 phương pháp:

1
 2
 3
 4
 5
 6
 7
 8
 9
10
def test(x, epsilon):
    assert x > 0
    assert epsilon > 0
    bisectionMethod(x, epsilon)
    newtonMethod(x, epsilon)
    print
test(100, 0.2)
test(10000, 0.0000002)
test(123456789,0.0000002)
test(0.25, 0.00002)

Ta nhận thấy phương pháp Newton tốt hơn so với Bisection. Tuy nhiên, với phương pháp Newton, chọn giá trị khởi đầu như thế nào sẽ là một vấn đề quan trọng.