sp-14. ニュートン法

sp-14. ニュートン法 （Scheme プログラミング） URL: https://www.kkaneko.jp/pro/scheme/index.html 金子邦彦 1

• https://www.kkaneko.jp/pro/scheme/index.html

アウトライン 14-1 ニュートン法 14-2 パソコン演習 14-3 課題 2

14-1 ニュートン法 3

本日の内容 • ニュートン法による非線形方程式の求解 • x0, x1, x2 ... の収束の様子を観察し，ニュートン法 の理解を深める • ニュートン法で，解が求まらない場合がありえるこ とを理解する 4

8 x3 – 6x2+11x –6 6 4 解 2 解 解 x 0 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 -2 -4 -6 -8 5

8 6 4 x 2 0 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 -2 -4 -6 -8 関数上の点 関数上の点を， １つ選ぶ 6

8 6 接線 4 x 2 0 0 0.5 1 -2 -4 関数上の点 1.5 2 2.5 3 3.5 4 4.5 接線を引く -6 -8 7

8 6 接線 4 接線と 2 x 軸の交点 x 解 0 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 -2 -4 -6 -8 関数上の点 接線と x 軸の交点は， 解の１つに近い（と 望むことができる） 8

8 6 接線 y = f '(a)(x - a) + f(a) 4 接線と x 軸の交点 (a - 2f(a)/f '(a), 0) x 0 0 0.5 1 1.5 2.5 3 関数上の点 -2 -4 2 3.5 4 4.5 (a, f (a)) から 関数上の点 (a, f(a)) -6 接線と x 軸の交点 -8 が求まる (a − f ' (a) / f (a),0) 9

ニュートン法 • 初期近似値 x0 から出発 • 反復公式 xi +1 f ( xi ) = xi − f ' ( xi ) を繰り返す • x1, x2, x3 ... は，だんだんと解に近づいていく （と望むことができる） 10

ニュートン法 • 初期近似値 x0 から出発 • 反復公式 xi +1 f ( xi ) = xi − f ' ( xi ) を繰り返す • x1, x2, x3 ... は，だんだんと解に近づいていく （と望むことができる） これは，点（x i, f(xi) ）における 接線と， x軸 (y=0) との交点 (xi+1, 0) を求めている 11

8 接線と x 軸の交点 (0.54545454..., 0) x14 = 0.54545454... 6 x 2 0 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 -2 -4 -6 -8 関数上の点 (0, -6) x0 = 0 12

8 6 接線と x 軸の交点 (0.84895321..., 0) x4 2 = 0.84895321... x 2 0 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 -2 -4 -6 関数上の点 (0.54545454..., -1.62283996...) x1 = 0.54545454... -8 13

8 接線と x 軸の交点 (0.97467407..., 0) x = 0.97467407... 4 3 6 x 2 0 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 -2 -4 -6 関数上の点 (0.84895321..., -0.37398512...) x2 = 0.84895321... -8 14

8 接線と x 軸の交点 (0.99909154..., 0) x = 0.999909154... 4 3 6 x 2 0 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 -2 -4 -6 関数上の点 (0.97460407..., -0.052592310...) x3 = 0.97467407... -8 15

ニュートン法の例 f(x) = x3 – 6x2+11x –6, x0 = 0 では： x0 = 0 f(x0) = -6, f '(x0) = 11 ⇒ x1 = x0 - f(x0)/f '(x0) = 0.54545454... x1 = 0.54545454... f(x1) = -1.62283996..., f '(x1) = 5.3471074... ⇒ x2 = x1 - f(x1)/f '(x1) = 0.84895321... x2 = 0.84895321... f(x2) = 0.37398512..., f '(x2) = 2.9747261... ⇒ x3 = x2 - f(x2)/f '(x2) = 0.97460407... 16

どうやって計算を終了するか 「反復公式」を繰り返すのだが ⇒ いつ止めたらいいのか？ ・ 決定打は無いのだが，今日の授業では次のやり方で 行ってみる ある小さな正の数δに対して f ( xn )   となった時点で計算を終了し xn を解とする 17

Scheme での多項式の書き方 • 「+」では，足したいものを３つ以上並べても よい 例） 2x2 - 3x - 1 (+ (* 2 x x) (* -3 x) -1)) 18

14-2 パソコン演習 19

パソコン演習の進め方 • 資料を見ながら，「例題」を行ってみる • 各自，「課題」に挑戦する • 自分のペースで先に進んで構いません 20

DrScheme の使用 • DrScheme の起動 プログラム • → PLT Scheme → DrScheme 今日の演習では「Intermediate Student」 に設定 Language → Choose Language → Intermediate Student → Execute ボタン 21

例題１．ニュートン法による 非線形方程式の解 • 非線型方程式 f(x) = 0 をニュート ン法で解く関数 newton を作り， 実行する • 方程式: f • 初期近似値： x0 22

「例題１．ニュートン法による非線形方程式の解」の手順 (1/2) 1. 次を「定義用ウインドウ」で，実行しなさい • 入力した後に，Execute ボタンを押す ;; d/dx: (number->number) number number -> number ;; inclination of the tangent (define (d/dx f x h) (/ (- (f (+ x h)) (f (- x h))) (* 2 h))) (define (is-good? f guess delta) (< (abs (f guess)) delta)) (define (improve f guess) (- guess (/ (f guess) (d/dx f guess 0.0001)))) ;; newton: (number->number) number number number -> number (define (newton f guess delta number) (cond [(or (is-good? f guess delta) (< number 0)) guess] [else (newton f (improve f guess) delta (- number 1))])) (define (f2 x) (+ (* x x x) (* -6 x x) (* 11 x) -6)) 23

「例題１．ニュートン法による非線形方程式の解」の手順 (2/2) 2. その後，次を「実行用ウインドウ」で実行しなさい (newton f2 #i0 0.00001 10000) ☆ 次は，課題に進んでください 24

まず，関数 newton などを定義している 25

次に関数 f2 を定義している (define (f2 x) (+ (* x x x) (* -6 x x) (* 11 x) -6)) 26

これは， (newton f2 #i0 0.00001 10000) と書いて，guess の値を #i0 に, delta の値を 0.0001 に， number の値を 10000 に f を f2 に設定しての実行 実行結果である 「#i0.9999987646860998」 が表示される 27

newton の入力と出力 f2 #i1 0.00001 10000 f2 = 0 の近似解 の１つ newton 入力 入力は１つの関数と， ３つの数値 出力 出力は数値 newton は，関数を入力とするような関数 （つまり高階関数） 28

;; d/dx: (number->number) number number -> number ;; inclination of the tangent (define (d/dx f x h) (/ (- (f (+ x h)) (f (- x h))) (* 2 h))) (define (is-good? f guess delta) (< (abs (f guess)) delta)) (define (improve f guess) (- guess (/ (f guess) (d/dx f guess 0.0001)))) ;; newton: (number->number) number number number -> number (define (newton f guess delta number) (cond [(or (is-good? f guess delta) (< number 0)) guess] [else (newton f (improve f guess) delta (- number 1))])) 29

ニュートン法の注意点 • 初期近似値の決め方 • 初期近似値によって，求まる解が変わってくる • 求まる解は，あくまでも「近似解」 例：この例題では #i0.9999987646860998 → #i 付きの近似計算で十分 • 虚数解は求まらない 30

(newton f2 #i0 0.00001 10000)の結果 #i0.9999987646860998 (newton f2 #i10 0.00001 10000)の結果 #i3.0000000168443197 （違う値が得られた） 31

(newton f2 0 0.00001 10000)の結果 ⇒ 結果は「有理数」で得られる （画面には表示しきれない） 32

「ニュートン法のプログラム」 の理解のポイント • 繰り返しの終了条件は２つ • 収束の条件を満足した f • 繰り返しの上限回数を超えた ( xn )   • 入力は４つ • 求めるべき関数: f • 初期近似値: guess • 収束条件を決める値: delta • 繰り返し回数の上限: number 33

ニュートン法の繰り返し処理 • x0 から始める • x1 を求める • x2 を求める … • 「収束の条件を満足する」か「繰り 返しの上限回数を超える」まで続け る 34

ニュートン法の繰り返し処理 1. 「収束した」あるいは「繰り返しの上限回数を超え た」 ならば： → 終了条件 現在の xn の値 → 自明な解 2. そうで無ければ： – 「接線と x 軸の交点の計算」を行う 求まった交点の x 座標の値は xn+1 – 結局，「収束する」か「繰り返しの上限 回数を超える」まで，接線と x 軸の交点 の計算を繰り返す 35

ニュートン法の繰り返し処理 (or (is-good? f guess delta) (< number 0)) No Yes 現在の guess の値が解 (newton f (improve guess) delta (- number 1)) 36

(newton f2 #i0 0.00001 10000) から結果が得られる過程 (newton f2 #i0 0.00001 10000) 最初の式 = ... = (newton (improve f2 #i0) 0.00001 9999) = ... = (newton #i0.5454545449588037 0.00001 9999) = ... = (newton (improve f2 #i0.5454545449588037) 0.00001 9998) = ... = (newton #i0.8489532098093157 0.00001 9998) = ... = (newton (improve f2 #i0.8489532098093157 ) 0.00001 9997) コンピュータ内部での計算 = ... = #i0.9999987646860998 実行結果 37

(newton f2 #i0 0.00001 10000) から結果が得られる過程 (newton f2 #i0 0.00001 10000) = ... = (newton (improve f2 #i0) 0.00001 9999) = ... =これは， (newton #i0.5454545449588037 0.00001 9999) = ... (define (newton f guess delta number) = (newton (improve f2 #i0.5454545449588037) 0.00001 9998) (cond = ... [(or (is-good? f guess delta) = (newton #i0.8489532098093157 (< number 0)) guess] 0.00001 9998) = ... [else (newton f (improve f guess) delta (- number 1))])) =の(newton (improve ) 0.00001 9997) f を f2 で，guess を f2 #i0#i0.8489532098093157 で，delta を 0.00001で, number を 10000 で置き換 =えたもの ... = #i0.9999987646860998 38

ニュートン法での判定基準 • ニュートン法では，現在の xn の誤差（どれだけ真 の値に近いか）は，正確には分からない 39

ニュートン法での収束条件 • is-good? が収束条件を判定 (define (is-good? f guess delta) (< (abs (f guess)) delta)) abs は「絶対値」 を求める f ( xn )   が成立するとき，出力は true. （さもなければ false） 40

ニュートン法での繰り返し処理 number: カウンタ（繰り返しの残り回数） guess: xn の値（計算の途中結果） number ← number － 1 guess ← (improve guess) を繰り返す 41

ニュートン法での繰り返し処理 (define (newton f guess delta number) (cond [(or (is-good? f guess delta) (< number 0)) guess] [else (newton f (improve f guess) delta (- number 1))])) guess ← (improve f guess) number ← number － 1 42

f(x) = x2 - 5 での x1, x2 ... の収束の様子 方程式 f(x) = x2 - 5 (newton f #i1 0.00001 10000) f(x)の導関数 f '(x) = 2x =… = (newton f #i3.0 0.00001 10000) =… = (newton f #i2.3333333333333335 0.00001 9999) =… = (newton f #i2.238095238095238 0.00001 9998) =… = (newton f #i2.2360688956433634 0.00001 9997) =… = (newton f #i2.236067977499978 0.00001 9996) =… guess delta number 43

ニュートン法の注意点 • 初期近似値の決め方 • 初期近似値の設定の際、あまりに解と掛け離れた 値を与えると、収束するのに時間がかかったり、 収束しなかったりする • delta の決め方 • どの程度の精度で計算するかを決定していないと， εが決まらない 44

ニュートン法の能力と限界 ニュートン法は, 出発点と する十分近い解を見付け ることができれば, 収束が 早い. f(x) ３ １ ２ 0 初期近似値の選び方次第 では，収束がしないこと がありえる 対策 x 関数ｆ(x)が単調でなくて変曲点を持つ （つまりｆ’(x)の符号が変わる）とき 例） 上の図では： ２： f'(x) = 0 ３： y 軸との交点が求まらない （負の無限大に発散） → 収束しない 繰り返しの上限回数 number を設定 45

14-3 課題 46

課題１ • 立方根を求めるプログラムを，Newton 法のプロ グラムを利用して作成せよ • f(x) = x3 - a = 0 を解く • a が負のときにも正しく負の立方根を求めることができ ることを確認せよ • delta の値を変えて実行を繰り返し，得られた解の値の 変化も報告しなさい 47

課題２．ニュートン法での収束 1. f(x) = x2 - 5 のグラフを書け（手書き） 2. 1. で作成したグラフに，ニュートン法での，x0, x1, x2, x3 の値を書き加えなさい 但し， • 関数 方程式 f(x) = x2 - 5 • 入力 初期値 x0 = 1 収束条件を決める値 delta = 0.00001 繰り返し回数の上限 number = 10000 x1, x2, x3 の値は，実際にニュートン法のプログラムを実 行して得ること 48

課題３ • (x-1)4(x-2) = 0 を newton 法で解け • 初期近似値を変えて実行を繰り返し，得られた解の値 の変化も報告しなさい 49

演習４ • Newton 法により f(x) = tan-1 x + 0.3x を解け • 初期近似値の選び方によっては，正しく解を求めるこ とができない．その理由についても考え，グラフを書 いて説明しなさい． 50

さらに勉強したい人への 補足説明事項 区間二分法 51

区間二分法 • 非線形方程式の求解の一手法 52

区間二分法 (half-interval method) の考え方 • 連続関数 f の性質 「f(a)<0<f(b)ならば，[a, b]の間に f(x) = 0 の解 がある f(x) f(b) a 0 f(a) 解 b x 53

区間二分法 (half-interval method) の処理手順 • 「区間」を半分ずつ縮小するような処理手順 1. ２点 a, b を定める 2. ２点a,bの中点(a+b)/2について： f((a+b)/2)>0ならば → 解は，[a, (a+b)/2]にある f((a+b)/2)<0ならば → 解は，[(a+b)/2, b]にある 3. 2. を繰り返す．「区間」の幅が小さくなる →解 の近似値を得られる f(b) a f(a) a+b 2 b 54

例題２．区間二分法による非線形方程式の解 • f(x)=x2 –2 を区間二分法で解くプ ログラム half-interval を書く • f(x) = x2 - 2 解は √２と－√２ 55

「例題２．区間二分法法による非線形方程式の解」の手順(1/2) 1. 次を「定義用ウインドウ」で，実行しなさい • 入力した後に，Execute ボタンを押す (define (f x) (- (* x x) 2)) (define (good-enough? a b) (< (- b a) 0.000001)) (define (middle a b) (/ (+ a b) 2)) (define (half-interval a b) (cond [(good-enough? a b) a] [else (cond [(< (f (middle a b)) 0) (half-interval (middle a b) b)] [(= (f (middle a b)) 0) (middle a b)] [(> (f (middle a b)) 0) (half-interval a (middle a b))])])) 56

「例題２．区間二分法による非線形方程式の解」の手順 (2/2) 2. その後，次を「実行用ウインドウ」で実行 しなさい • 正しい解が得られている場合と，得られていな い場合があることを確認しなさい (half-interval 0 2) (half-interval -2 0) (half-interval 2 4) 57

「区間二分法による非線形方程式の解」の実行結果 > (half-interval 0 2) 1.4142131805419921875 → 正しい > (half-interval -2 0) -0.00000095367431640625 → 正しくない > (half-interval 2 4) 2 → 正しくない 58

入力と出力 a の値: 0 b の値: 2 1.4142131805419921875 half-interval 入力 出力 入力は２つの数値 出力は数値 59

区間二分法のプログラム • 関数 • 方程式 f(x) • 入力 • 初期値 a, b • 出力 • 解 60

区間二分法の処理手順 1. 初期値 a, b の決定 2. 区間を半分に縮小 • ｆ((a+b)/2) の値が 正： b ← ((a+b)/2 ０： (a+b)/2 が解である 負： a ← ((a+b)/2 3. 2. を繰り返す 61

区間二分法の処理手順 f(x) = x2-2, a = 0, b = 2 の場合 a = 0, b = 2 f((a+b)/2) = f(1) = -1 ⇒ a を「1」に a = 1, b = 2 f((a+b)/2) = f(1.5) = 0.25 ⇒ b を「1.5」に a = 1, b = 1.5 f((a+b)/2) = f(1.25) = -0.4375 ⇒ a を「1.25」に 62

区間二分法の繰り返し処理 • half-interval の内部に half-interval が登場 (define (half-interval a b) (cond [(good-enough? a b) a] [else (cond [(< (f (middle a b)) 0) (half-interval (middle a b) b)] [(= (f (middle a b)) 0) (middle a b)] [(> (f (middle a b)) 0) (half-interval a (middle a b))])])) • half-interval の実行が繰り返される 63

区間二分法での収束条件 • good-enough? が収束条件を判定 (define (good-enough? a b) (< (- b a) 0.000001)) b - a < 0.000001 が成立するとき，出力は true. （さもなければ false） * 「0.000001」 は適当に決めた値 64

「区間二分法のプログラム」 の理解のポイント • 繰り返しの終了条件は２つ • b – a の値が，「しきい値」を下回った（収束した） • 「ｆ((a+b)/2) = 0」となった • 入力は２つ • 初期値 a, b • 繰り返し処理： 反復的プロセス 65

区間二分法での注意点 • f(a) ≧0, f(b)≦0 でないと，解が得られな い f(x) f(b) a 0 f(a) 解 b x 66

例題３．区間二分法での繰り返し処理 • 例題２の「区間二分法」のプログラムについて， (half-interval 0 2) から(half-interval (middle 1 3/2) 3/2) まで の過程を確認する (half-interval 0 2) =… = (half-interval (middle 0 2) 2) =… = (half-interval 1 2) =… = (half-interval 1 (middle 1 2)) =… = (half-interval 1 3/2) =… = (half-interval (middle 1 3/2) 3/2) 67

「例題３．区間二分法での繰り返し処理」の手順 1. 次を「定義用ウインドウ」で，実行しなさい • • Intermediate Student で実行すること 入力した後に，Execute ボタンを押す (define (f x) (- (* x x) 2)) (define (good-enough? a b) (< (- b a) 0.000001)) (define (middle a b) (/ (+ a b) 2)) (define (half-interval a b) 例題２に (cond １行書き加える [(good-enough? a b) a] [else (cond [(< (f (middle a b)) 0) (half-interval (middle a b) b)] [(= (f (middle a b)) 0) (middle a b)] [(> (f (middle a b)) 0) (half-interval a (middle a b))])])) (half-interval 0 2) 2. DrScheme を使って，ステップ実行の様子を 確認しなさい （Step ボタン，Next ボタンを使用） • 理解しながら進むこと 68