ONLY DO WHAT ONLY YOU CAN DO

こけたら立ちなはれ 立ったら歩きなはれ

ジュリアに傷心 - Juliaで非線形方程式を解く (2分法)

非線形方程式の解法(2分法)を利用して2の平方根を求める

f:id:fornext1119:20140531160605p:plain
1. まず, 条件 a < b, f(a) < 0, f(b) > 0 を満たす点 a, b を考えると,
関数 f(x) の解は, 区間 (a,b) の中に存在する.
2. 次に, 区間 (a,b) の中点 c = (a + b) / 2 を考えると,
f(c) < 0 であれば, 解は区間 (c,b) の中に存在し,
同様に, f(c) > 0 であれば, 区間 (a,c) の中に存在する.
3. この作業を繰り返して, 区間を絞り込んで行くことで解を求める.

※ちなみに, 区間 (a,b) とは, 両端を含まない数の集合で, 開区間という
 (a, b) = \{x | a < x < b\}
また, 両端を含む数の集合は, 閉区間という
 [a, b] = \{x | a ≤ x ≤ b\}

function f(x)
    x * x - 2
end

function bisection(a::Float64, b::Float64)
    c = 0
    while (true)
        # 区間 (a, b) の中点 c = (a + b) / 2
        c = (a + b) / 2
        @printf("%12.10f\t%13.10f\n", c, c - sqrt(2))

        fc = f(c)
        if abs(fc) < 0.0000000001
            break
        end

        if (fc < 0)
            # f(c) < 0 であれば, 解は区間 (c, b) の中に存在
            a = c
        else
            # f(c) > 0 であれば, 解は区間 (a, c) の中に存在
            b = c
        end
    end
    c
end

a = 1.0
b = 2.0
@printf "%12.10f\n" bisection(a, b)
1.5000000000     0.0857864376
1.2500000000    -0.1642135624
1.3750000000    -0.0392135624
1.4375000000     0.0232864376
1.4062500000    -0.0079635624
1.4218750000     0.0076614376
1.4140625000    -0.0001510624
1.4179687500     0.0037551876
1.4160156250     0.0018020626
1.4150390625     0.0008255001
1.4145507813     0.0003372189
1.4143066406     0.0000930783
1.4141845703    -0.0000289921
1.4142456055     0.0000320431
1.4142150879     0.0000015255
1.4141998291    -0.0000137333
1.4142074585    -0.0000061039
1.4142112732    -0.0000022892
1.4142131805    -0.0000003818
1.4142141342     0.0000005718
1.4142136574     0.0000000950
1.4142134190    -0.0000001434
1.4142135382    -0.0000000242
1.4142135978     0.0000000354
1.4142135680     0.0000000056
1.4142135531    -0.0000000093
1.4142135605    -0.0000000019
1.4142135642     0.0000000019
1.4142135624     0.0000000000
1.4142135624
参考文献