Line | |
---|
1 | #
|
---|
2 | # newton.rb
|
---|
3 | #
|
---|
4 | # Solves the nonlinear algebraic equation system f = 0 by Newton's method.
|
---|
5 | # This program is not dependent on BigDecimal.
|
---|
6 | #
|
---|
7 | # To call:
|
---|
8 | # n = nlsolve(f,x)
|
---|
9 | # where n is the number of iterations required,
|
---|
10 | # x is the initial value vector
|
---|
11 | # f is an Object which is used to compute the values of the equations to be solved.
|
---|
12 | # It must provide the following methods:
|
---|
13 | #
|
---|
14 | # f.values(x):: returns the values of all functions at x
|
---|
15 | #
|
---|
16 | # f.zero:: returns 0.0
|
---|
17 | # f.one:: returns 1.0
|
---|
18 | # f.two:: returns 1.0
|
---|
19 | # f.ten:: returns 10.0
|
---|
20 | #
|
---|
21 | # f.eps:: returns the convergence criterion (epsilon value) used to determine whether two values are considered equal. If |a-b| < epsilon, the two values are considered equal.
|
---|
22 | #
|
---|
23 | # On exit, x is the solution vector.
|
---|
24 | #
|
---|
25 | require "bigdecimal/ludcmp"
|
---|
26 | require "bigdecimal/jacobian"
|
---|
27 |
|
---|
28 | module Newton
|
---|
29 | include LUSolve
|
---|
30 | include Jacobian
|
---|
31 |
|
---|
32 | def norm(fv,zero=0.0)
|
---|
33 | s = zero
|
---|
34 | n = fv.size
|
---|
35 | for i in 0...n do
|
---|
36 | s += fv[i]*fv[i]
|
---|
37 | end
|
---|
38 | s
|
---|
39 | end
|
---|
40 |
|
---|
41 | def nlsolve(f,x)
|
---|
42 | nRetry = 0
|
---|
43 | n = x.size
|
---|
44 |
|
---|
45 | f0 = f.values(x)
|
---|
46 | zero = f.zero
|
---|
47 | one = f.one
|
---|
48 | two = f.two
|
---|
49 | p5 = one/two
|
---|
50 | d = norm(f0,zero)
|
---|
51 | minfact = f.ten*f.ten*f.ten
|
---|
52 | minfact = one/minfact
|
---|
53 | e = f.eps
|
---|
54 | while d >= e do
|
---|
55 | nRetry += 1
|
---|
56 | # Not yet converged. => Compute Jacobian matrix
|
---|
57 | dfdx = jacobian(f,f0,x)
|
---|
58 | # Solve dfdx*dx = -f0 to estimate dx
|
---|
59 | dx = lusolve(dfdx,f0,ludecomp(dfdx,n,zero,one),zero)
|
---|
60 | fact = two
|
---|
61 | xs = x.dup
|
---|
62 | begin
|
---|
63 | fact *= p5
|
---|
64 | if fact < minfact then
|
---|
65 | raise "Failed to reduce function values."
|
---|
66 | end
|
---|
67 | for i in 0...n do
|
---|
68 | x[i] = xs[i] - dx[i]*fact
|
---|
69 | end
|
---|
70 | f0 = f.values(x)
|
---|
71 | dn = norm(f0,zero)
|
---|
72 | end while(dn>=d)
|
---|
73 | d = dn
|
---|
74 | end
|
---|
75 | nRetry
|
---|
76 | end
|
---|
77 | end
|
---|
Note:
See
TracBrowser
for help on using the repository browser.