sp-14. ニュートン法

186 Views

January 26, 22

スライド概要

(Scheme プログラミング)
URL: https://www.kkaneko.jp/cc/scheme/index.html

profile-image

金子邦彦(かねこくにひこ) 福山大学・工学部・教授 ホームページ: https://www.kkaneko.jp/index.html 金子邦彦 YouTube チャンネル: https://youtube.com/user/kunihikokaneko

シェア

埋め込む »CMSなどでJSが使えない場合

各ページのテキスト
1.

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

2.

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

3.

14-1 ニュートン法 3

4.

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

5.

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

6.

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 関数上の点 関数上の点を, 1つ選ぶ 6

7.

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.

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 軸の交点は, 解の1つに近い(と 望むことができる) 8

9.

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 (a, f (a)) 4.5 から 関数上の点 (a, f(a)) -6 接線と x 軸の交点 -8 が求まる (a − f ' (a) / f (a),0) 9

10.

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

11.

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

12.

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

13.

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

14.

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

15.

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

16.

ニュートン法の例 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

17.

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

18.

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

19.

14-2 パソコン演習 19

20.

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

21.

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

22.

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

23.

「例題1.ニュートン法による非線形方程式の解」の手順 (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

24.

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

25.

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

26.

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

27.

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

28.

newton の入力と出力 f2 #i1 0.00001 10000 f2 = 0 の近似解 の1つ newton 入力 入力は1つの関数と, 3つの数値 出力 出力は数値 newton は,関数を入力とするような関数 (つまり高階関数) 28

29.

;; 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

30.

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

31.

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

32.

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

33.

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

34.

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

35.

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

36.

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

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) = ... = (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

38.

(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

39.

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

40.

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

41.

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

42.

ニュートン法での繰り返し処理 (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

43.

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

44.

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

45.

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

46.

14-3 課題 46

47.

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

48.

課題2.ニュートン法での収束 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

49.

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

50.

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

51.

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

52.

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

53.

区間二分法 (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

54.

区間二分法 (half-interval method) の処理手順 • 「区間」を半分ずつ縮小するような処理手順 1. 2点 a, b を定める 2. 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

55.

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

56.

「例題2.区間二分法法による非線形方程式の解」の手順(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

57.

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

58.

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

59.

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

60.

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

61.

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

62.

区間二分法の処理手順 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

63.

区間二分法の繰り返し処理 • 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

64.

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

65.

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

66.

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

67.

例題3.区間二分法での繰り返し処理 • 例題2の「区間二分法」のプログラムについて, (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

68.

「例題3.区間二分法での繰り返し処理」の手順 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) 例題2に (cond 1行書き加える [(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