gistfile1.py

#!/usr/bin/python
# -*- coding: utf-8 -*-

'''
    非線形方程式の根を求める
        再帰的に近似を求め関数の根を求める

    根を求める関数
        f(x) = x^3 - x + 1

    使用アルゴリズム一覧
        ・ 二分法
        ・ ニュートン・ラフソン法
        ・ 割線法
        ・ 挟み撃ち法
        ・ Ridder法
'''

from random import random
import math

fn = lambda x: x ** 3 - x + 1
fnd = lambda x: 3 * (x ** 2) - 1
EPS = 1e-8


def sign(x):
    return 1 if x >= 0 else -1


def findvar(fn, match):
    tmp = 0
    while match(fn(tmp)) == False:
        tmp = random() * 100 - 50
    return tmp


def ridders(fn):
    '''
        Ridder法
            挟み撃ち法を改良したもの

        段取り
            正負をまたいだ二数{xl, xh}を用意する
            xl, xh の中点 を xm と定める
            その三点を次の式に代入して次の点を得る
            new = xm + (xm - xl) *
                        (sign[f(xl) - f(xh)] * f(xh)) /
                        sqrt(f(xm)^2 - f(xl) * f(xh))
            この新しい点は同符号の xl, xh と差し替える。
            上記を繰り返し十分な精度まで達したら終える
    '''
    xl = findvar(fn, lambda x: x <= 0)
    xh = findvar(fn, lambda x: x >= 0)
    newx = 0.0
    fnewx = fn(newx)

    if fn(xl) == 0: return xl
    elif fn(xh) == 0: return xh
    elif fn(newx) == 0: return newx

    while fnewx > EPS:
        xm = (xl + xh) * 0.5
        fl, fm, fh = fn(xl), fn(xm), fn(xh)
        s = math.sqrt(fm ** 2.0 - fl * fh)
        newx = xm + (xm - xl) * (sign(fl - fh) * fm) / s
        fnewx = fn(newx)
        if fnewx == 0:
            break
        if fnewx * fl >= 0:
            xl = newx
        elif fnewx * fh >= 0:
            xh = newx
    return newx


def falseposition(f):
    '''
        挟み撃ち法
            割線法に二分法を加えた物
    '''
    x, dx = 0, EPS + 1
    x0 = findvar(fn, lambda x: x < 0)
    x1 = findvar(fn, lambda x: x >= 0)
    fx, fx0, fx1 = f(x), f(x0), f(x1)

    if fx == 0: return x
    elif fx0 == 0: return x0
    elif fx1 == 0: return x1

    while dx >= EPS:
        # x = x1 - fx1 * (x1 - x0) / (fx1 - fx0)
        x = (fx1 * x0 - fx0 * x1) / (fx1 - fx0)
        fx = fn(x)
        if fx == 0:
            break
        if fx * fx0 >= 0:
            dx, x0, fx0 = abs(x - x0), x, fx
        elif fx * fx1 >= 0:
            dx, x1, fx1 = abs(x - x1), x, fx
    return x


def secant(f):
    '''
        割線法
            {x0, x1} を通る直線と x 軸の交点を再帰的に求める
    '''
    x, x0, x1 = 0, math.pi, math.e
    f0, f1 = f(x0), f(x1)

    if f0 < f1:
        x0, x1 = x1, x0
        f0, f1 = f0, f1

    while abs(f0 - f1) > EPS:
        x = x1 - f1 * (x1 - x0) / (f1 - f0)
        x0, x1 = x1, x
        f0, f1 = f(x0), f(x1)
    return x


def newton(f, fd):
    '''
        ニュートン・ラフソン法
        ニュートン法によって fn(x) = 0 を求める
        fnd = fn の部分関数
    '''
    def nextx(f, fd, x):
        # Reccurrence
        return x - f(x) / fd(x)

    x0 = math.pi
    x1 = nextx(f, fd, x0)
    f0, f1 = f(x0), f(x1)
    while abs(f0 - f1) > EPS:
        x0 = x1
        x1 = nextx(f, fd, x0)
        f0, f1 = f(x0), f(x1)
    return x1


def bisection(f):
    '''
        二分法
            二分法で fn(x) = 0 を求める
    '''
    EPS = 1e-8
    xl = findvar(f, lambda x: x <= 0)
    xh = findvar(f, lambda x: x >= 0)
    xm = 0.0

    if f(xl) == 0: return xl
    elif f(xm) == 0: return xm
    elif f(xh) == 0: return xh

    while f(xm) > EPS:
        xm = (xl + xh) * 0.5
        fl, fm, fh = f(xl), f(xm), f(xh)
        if fm == 0:
            break
        if fl * fm >= 0:
            xl, fl = xm, fm
        elif fm * fh >= 0:
            xh, fh = xm, fm
    return xm


def main():
    print 'Bisection:', bisection(fn)
    print 'Newton-Raphson:', newton(fn, fnd)
    print 'Secant method:', secant(fn)
    print 'False Position method:', falseposition(fn)
    print 'Ridders\' low:', ridders(fn)


if __name__ == '__main__':
    main()