2.4
f(λ) = -1 + λ + λ2 (1)
f'(λ) = 1 + 2λ (2)
式(1)のグラフをFig. 2.4に示す。
a graph of eq.(1) is shown in Fig. 2.4.
Hand calculation
λ0 = 1 ( >0 )
f(λ0) = -1 + 1 + 12 = 1
f'(λ0) = 1 + 2(1) = 1 + 2 = 3
λ1 = λ0 - f(λ0) / f'(λ0) = 1 - 1 / 3 = 2/3
f(λ1) = -1 + 2/3 + (2/3)2 = -1 + 2/3 + 4/9 = 1/9
f'(λ1) = 1 + 2(2/3) = 1 + 4/3 = 7/3
λ2 = λ1 - f(λ1) / f'(λ1) = 2/3 - (1/9) / (7/3) = 13/21
f(λ2) = -1 + 13/21 + (13/21)2 = 1/212
f'(λ2) = 1 + 2(13/21) = 47/21
λ3 = λ2 - f(λ2) / f'(λ2) = 13/21 - (1/212) / (47/21) = 610/987
f(λ3) = -1 + 610/987 + (610/987)2 = 1/9872
f'(λ3) = 1 + 2(610/987) = 2207/987
λ4 = λ3 - f(λ3) / f'(λ3) = 610/987 - (1/9872)/(2207/987) = 1346269/2178309
∴λ0 = 1
λ1 = 2/3 (=0.667)
λ2 = 13/21 (=0.619)
λ3 = 610/987 (=0.618)
λ4 = 1346269/2178309 (=0.618)
Newton法の収束の様子をFig. 2.4.2に示す。
Maxima
Maximaによる解を次に示す。(added on 2016-11-05)
Maximaによる解は、Hand calculationによる解と一致した。
注
%i Maxima に対する手入力。
%o Maxima による出力。
%i1 関数 f(x) を定義(:=)する。
%i3 変数 x0 に 1 を代入(:)する。
Maxima 5.38.1 http://maxima.sourceforge.net
using Lisp SBCL 1.3.4
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.
(%i1) f(x) := -1 + x + x ^ 2;
2
(%o1) f(x) := (- 1) + x + x
(%i2) fd(x) := 1 + 2 * x;
(%o2) fd(x) := 1 + 2 x
(%i3) x0 : 1;
(%o3) 1
(%i4) x1 : x0 - f(x0) / fd(x0);
2
(%o4) -
3
(%i5) x2 : x1 - f(x1) / fd(x1);
13
(%o5) --
21
(%i6) x3 : x2 - f(x2) / fd(x2);
610
(%o6) ---
987
(%i7) x4 : x3 - f(x3) / fd(x3);
1346269
(%o7) -------
2178309
(%i8)
f(x)=0の根を「探索」で求める
added on 2016-11-13.
二分法(bisection method)で方程式の根を求める方法が、二分探索(binary search)
に似ているので、探索の手法で方程式の根を求めるプログラムを作成した。
f(x)を次式とする。
f(x) = x2 + x - 1
区間[0, 1]において、間隔0.001毎にxとf(x)を計算し、表を作成する。 間隔0.001は、方程式の根を小数点以下3桁まで求めるために設定した。
| i | 0 | 1 | … | i | … | 999 | 1000 |
| xi | 0 | 0.001 | … | xi | … | 0.999 | 1 |
| f(xi) | 0 | -0.998999 | … | f(xi) | … | 0.997001 | 1 |
逐次探索(sequential search)
上表のデータを次で作成する。(in Python 3)
x_table = list(map(lambda x: x / 1000, range(1001))) # 0.000 to 1.000 step 0.001 f_table = list(map(lambda x: f(x), x_table)) fd_table = list(map(lambda x: fd(x), x_table))
次の関数でリストtab、f(xi)のリスト、の個々の要素を順番に調べていって、
f(xi) が 0 の要素を探索する。
ただし、f(xi) が 0 である要素がない場合は、
f(xi)の絶対値が最小になる xi を方程式の根とする。
関数 lookup_a(tab) は、リストtab のインデックス i と探索した回数 n をリスト
[i, n]で返す。
# sequential search def lookup_a(tab): n = 0 if tab[0] < 0: for i in range(len(tab)): n = n + 1 if tab[i] >= 0: break if tab[0] > 0: for i in range(len(tab)): n = n + 1 if tab[i] <= 0: break if tab[0] == 0: return [0, n] if abs(tab[i-1]) < abs(tab[i]): return [i-1, n] else: return [i, n]
二分探索(binary search)
次の関数では、リストtab、f(xi)のリスト、
の最初の要素のインデックスを low とし、最後の要素のインデックスを high とし、
二分探索をリストtabに適用する。
探索中に low + 1 = high になった場合は、
f(xi) の絶対値の小さい方を根とする。
関数 lookup_b()は、リストのインデックス i と探索した回数 n をリスト[i, n]で返す。
# binary seaarch def lookup_b(tab): n = 0 low = 0 high = len(tab) - 1 while low <= high: n = n + 1 mid = int((low + high) / 2) cmp = tab[mid] print('{:2d}{:4d}{}{:}'.format(n, mid, ' ',cmp)) if cmp > 0: high = mid elif cmp < 0: low = mid else: return [mid, n] if low + 1 == high: if abs(tab[low]) < abs(tab[high]): return [low, n] else: return [high, n]
Newtonの方法
次のプログラムは、Newtonの方法を上の数表に適用したものである。 小数点以下3桁が変化しなくなることを収束条件とした。
# Newton's search def lookup_c(x_tab, f_tab, fd_tab): n = 0 if f_tab[len(f_tab)-1] > 0: i1 = len(f_tab) - 1 n = n + 1 x2 = x_tab[i1] - f_tab[i1] / fd_tab[i1] while 1: x2 = round(x2, 3) i2 = x_tab.index(x2) print('{:10d}{:10d}{:10d}{:15.3f}'.format(n,i1,i2,x2)) if i1 == i2: return [i, n] i1 = i2 n = n + 1 x2 = x_tab[i1] - f_tab[i1] / fd_tab[i1] return [i, n]
探索結果
探索結果を次に示す。
逐次探索は620回、二分探索は10回、 Newtonの方法は4回で小数点以下3桁までの根が探索された。
sequential search result
n i x_table[i] f_table[i]
--------------------------------------------------
620 618 0.618 -7.599999999996498e-05
--------------------------------------------------
binary search
n i f(x)
------------------------------
1 500 -0.25
2 750 0.3125
3 625 0.015625
4 562 -0.12215599999999993
5 593 -0.05535100000000004
6 609 -0.020118999999999998
7 617 -0.002310999999999952
8 621 0.006641000000000119
9 619 0.0021610000000000795
10 618 -7.599999999996498e-05
------------------------------
binary search result
n i x_table[i] f_table[i]
--------------------------------------------------
10 618 0.618 -7.599999999996498e-05
--------------------------------------------------
Newton's search
n i1 i2 x_table[i]
---------------------------------------------
1 1000 667 0.667
2 667 619 0.619
3 619 618 0.618
4 618 618 0.618
---------------------------------------------
Newton's search result
n i x_table[i] f_table[i]
--------------------------------------------------
4 618 0.618 -7.599999999996498e-05
--------------------------------------------------
探索プログラム
注、逐次探索と二分探索の部分は、次を見て作成した。
Brian W.Kernighan, Rob Pike, 福崎俊博 訳、プログラミング作法、株式会社アスキー、2000、pp.54-56 、ISBN4-7561-3649-4.
Python 3 code
Fortran code
revised on 201-03-26, 2025-02-01 open revised 2.4.f90
2.4.f90Gauche code
revised on 201-03-26, 2025-02-01 open revised 2.4.scm
2.4.scmJavaScript code
added on 2016-02-06, 2025-02-01 open revised 2.4.js
2.4.jsJulia code
コードを次のリンクに示す。
2.4.jl
備考 Fig. 2.4他について
Note, about Fig. 2.4 and so on
(added on 2005-2-14, revised on 2011-10-29.)
Fig.2.4は、次を使用して作成した。
Fig. 2.4 was made with the following.
GNUPLOT Version 4.4 patchlevel 3 last modified March 2011 System: MS-Windows 32 bit Copyright (C) 1986-1993, 1998, 2004, 2007-2010 Thomas Williams, Colin Kelley and many others
図化コマンドを次のリンクに示す。
history
revised on 2005-2-14, 2007-04-21, 2008-02-19, 2009-08-09, 2011-10-26,
2011-11-02, 2014-02-11, 2015-03-28, 2016-02-06, 2016-03-26, 2016-09-28,
2016-10-23, 2016-10-30.
2016-09-28 Fortarn, Gauche, Javascript codes moved to link.
2016-10-23 Fig. Newton's method added.
2016-10-30 Python code and gnuplot script moved to link.
2016-11-05 solution by Maxima added.
2016-11-13 solution in search way added.
2018-01-01 png changed to svg.
2018-01-03 XHTML1.0 changed to HTML5.
2018-06-03 charset shift-jis to utf-8.
2021-02-14 move history.
2024-08-28 revise 2.4.js.
2024-10-26 add 2.4.jl.
2025-02-01 open revised 2.4.f90, 2.4.scm, and 2.4.js.