Blame view

chap2/newton.py 1.6 KB
e5dbe416   Chunk   staged.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
__author__ = 'chunk'

import numpy as np
import math


def deriv(fn):
    dx = 0.00001
    return lambda x: (fn(x + dx) - fn(x)) / dx


def dampnm(fun, x0, lmbda, epsilon, gfun=None):
    assert lmbda < 1 and lmbda > 0

    if gfun == None:
        gfun = deriv(fun)
    x_new = x0
    x_old = None
    # while np.absolute(fun(x_new)) > epsilon:
    while x_old == None or np.absolute(fun(x_new) - fun(x_old)) > epsilon:
        x_old = x_new
        s = fun(x_old) / gfun(x_old)
        x_new = x_old - s
        while np.absolute(fun(x_new)) >= np.absolute(fun(x_old)):
            x_new = x_old - lmbda * s
            lmbda /= 2.0
        print lmbda, x_new
    return x_new


def test():
    fn = lambda x: math.pow(x, 4) - x - 2
    gfn = lambda x: 4 * math.pow(x, 3) - 1
    print dampnm(fn, 1.5, 0.99, 0.00001, gfn)

    fn = lambda x: math.pow(x, 3) - x - 1
    gfn = lambda x: 3 * math.pow(x, 2) - 1
    print dampnm(fn, 0.6, 0.99, 0.00001, gfn)

    fn = lambda x: -math.pow(x, 3) + 5 * x
    gfn = lambda x: -3 * math.pow(x, 2) + 5
    print dampnm(fn, 1.2, 0.99, 0.00001, gfn)

    # 0.99 1.375
    # 0.99 1.35377701579
    # 0.99 1.35321036029
    # 0.99 1.3532099642
    # 1.3532099642
    # 0.01546875 1.13521875
    # 0.01546875 1.36976137194
    # 0.01546875 1.32649943078
    # 0.01546875 1.32472090757
    # 0.01546875 1.32471795725
    # 0.01546875 1.32471795724
    # 1.32471795724
    # 0.2475 -1.90976470588
    # 0.2475 -2.34458457506
    # 0.2475 -2.24316474291
    # 0.2475 -2.23610151441
    # 0.2475 -2.23606797825
    # 0.2475 -2.2360679775
    # -2.2360679775

if __name__ == '__main__':
    test()

    pass