sp-13. 数値微分と数値積分

>100 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-13. 数値微分と数値積分 (Scheme プログラミング) URL: https://www.kkaneko.jp/cc/scheme/index.html 金子邦彦 1

2.

アウトライン 13-1 数値微分と数値積分 13-2 パソコン演習 13-3 課題 2

3.

13-1 数値微分と数値積分 3

4.

内容 1. 接線の傾き d/dx 2. 台形則による数値積分 trapezoid 4

5.

接線の傾き f(x) 接線 f ( x + h) − f ( x − h) 2h f ( x + h) f ( x − h) 0 傾きは x−h x+h x h を 0 に近づけることで、接線の傾き (=f '(x) )の「近似値」が求まる 5

6.

定積分 f(x) a 定積分: I = b  b a x f ( x ) dx 区間 [a, b] で,連続関数 f, x軸, x=a, x=b で囲まれた面積 6

7.

区間 [a, b] の小区間への分割 f(x) x1x2 x0 = a xn-1 x xn = b • n 個の等間隔な小区間に分割 • 幅: h = (b-a) / n • 小区間: [x0, x1], [x1, x2], ..., [xn-1, xn] 但し,x0 = a, xi = x0 + i × h 7

8.

小長方形 f(x) f(xi) x0 = a xi xi+1 小長方形 x xn = b • 小長方形の面積は f ( xi )h 8

9.

小台形 f(x) f(xi) x0 = a f(xi+ 1) xi xi+1 小台形 x xn = b • 小台形の面積は h ( f ( x i ) + f ( x i + 1 )) 2 9

10.

台形則 trapezoidal rule f(x) x1x2x3 xn-1 xn-2 x0 = a x xn = b • 小台形の面積の和は Sn = n −1  i=0 h ( f ( x i ) + f ( x i + 1 )) 2 • 定積分 I を,この和 Sn で近似 ⇒ 台形則という 10

11.

台形則による数値積分 • 区間[a,b]を n 等分 (1区間の幅h=(b-a)/n) • n 個の台形を考え, その面積の和 Sn で,定積分 I を近 似 • f(x) が連続関数のときは,n を無限大に近づければ,Sn は I に近づく I = lim n→  n −1  i=0 h ( f ( x i ) + f ( x i +1 )) = lim S n n→  2 • 但し,単純に「n を大きくすればよい」とは言えない • n を大きくすると ⇒ 計算時間の問題,丸め誤差の問題が 発生 11

12.

数値積分の意味 • 式で書いてある関数の積分 ⇒ ごく限られた関数しか定積分できない • 「数値積分」が重要になる 12

13.

台形則 • 両端 x0 = a と xn = b を除いて,f(xi) は2度現 れる Sn = n −1  i=0 h = 2 h ( f ( x i ) + f ( x i +1 ) ) 2 ( f ( x 0 ) + 2 ( f ( x1 ) +  + 1 = h   f ( x0 ) + 2 1 = h   f (a ) + 2 n −1  i =1 n −1  f ( x n − 1 )) + f ( x n ) ) 2回現れる部分 1  f ( xi ) + f ( xn )  2  f ( a + ih ) + i =1 但し 1  f (b )  2  h = (b-a) / n 13

14.

13-2 パソコン演習 14

15.

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

16.

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

17.

例題1.接線の傾き • 接線の傾きを求める関数 d/dx を作り, 実行する • 数値 x, h と関数 f から,x における f の 傾き(= f' (x))の近似値を求める • d/dx は高階関数 17

18.

「例題1.接線の傾き」の手順 1. 次を「定義用ウインドウ」で,実行しなさい • 入力した後に,Execute ボタンを押す ;; d/dx: (number->number) number number -> number ;; inclination of the tangent ;; Example: (d/dx f 3 0.0001) ;; = ((- (f 3.0001) (f 2.9999)) 0.0002) (define (d/dx f x h) (/ (- (f (+ x h)) (f (- x h))) (* 2 h))) (define (f2 x) (- (* x x) 2)) 2. その後,次を「実行用ウインドウ」で実行しなさい (d/dx f2 3 0.0001) ☆ 次は,例題2に進んでください 18

19.

まず,関数 d/dx を定義している 19

20.

次に関数 f2 を定義している (define (f2 x) (- (* x x) 2)) 20

21.

これは, (d/dx f2 3 0.0001) と書いて,x の値を 3 に, h の値を 0.0001 に, f を f2 に設定しての実行 実行結果である「6」が 表示される 21

22.

入力と出力 f2 3 0.0001 f2'(3) の近似値 d/dx 入力 入力は2つの数値と 関数 出力 出力は数値 d/dx は,関数を入力とするような関数 (つまり高階関数) 22

23.

d/dx 関数 ;; d/dx: (number->number) number number -> number ;; inclination of the tangent ;; Example: (d/dx f 3 0.0001) ;; = ((- (f 3.0001) (f 2.9999)) 0.0002) (define (d/dx f x h) (/ (- (f (+ x h)) (f (- x h))) (* 2 h))) 関数が, 「d/dx」の入力になっている 23

24.

(dx/d f2 3 0.0001) から 6 が得られる過程の概略 (dx/d f2 3 0.0001) 最初の式 = (/ (- (f2 (+ 3 0.0001)) (f2 (- 3 0.0001))) (* 2 0.0001)) = (/ (- (- (* (+ 3 0.0001) (+ 3 0.0001)) 2) (- (* (- 3 0.0001) (- 3 0.0001)) 2) (* 2 0.0001))) コンピュータ内部での計算 =… =6 実行結果 但し,f2 は (define (f2 x) (- (* x x) 2)) 24

25.

(dx/d f2 3 0.0001) から 6 が得られる過程の概略 (dx/d f2 3 0.0001) = (/ (- (f2 (+ 3 0.0001)) (f2 (- 3 0.0001))) (* 2 0.0001)) = (/ (- (- (* (+ 3 0.0001) (+ 3 0.0001)) 2) これは, (- (* (- 3 0.0001) (- 3 0.0001)) 2) (define (d/dx f x h) (*(/2 0.0001))) (- (f (+ x h)) (f (- x h))) (* 2 h))) =… の x を 3 で,h を 0.0001 で, f を f2 で置き換 =6 えたもの 25

26.

(dx/d f2 3 0.0001) から 6 が得られる過程の概略 (dx/d f2 3 0.0001) = (/ (- (f2 (+ 3 0.0001)) (f2 (- 3 0.0001))) (* 2 0.0001)) = (/ (- (- (* (+ 3 0.0001) (+ 3 0.0001)) 2) (- (* (- 3 0.0001) (- 3 0.0001)) 2) (* 2 0.0001))) =… これは, = 6(define (f2 x) (- (* x x) 2)) の x を (+ 3 0.0001) や (- 3 0.0001) で置き換えたも の 26

27.

例題2.台形則による数値積分 • 台形則による数値積分の関数 trapezoid を作り,実行する • 但し, 積分したい関数 f(x) = e 積分区間: [a, b] 区間数: n -x2 27

28.

「例題2.台形則による数値積分」の手順 1. 次を「定義用ウインドウ」で,実行しなさい • 入力した後に,Execute ボタンを押す (define (trapezoid-iter f a h k) (cond [(= k 1) (f (+ a h))] [else (+ (f (+ a (* h k))) (trapezoid-iter f a h (- k 1)))])) ;; trapezoid: number number number -> number ;; to compute the area under the graph of f between a and b (define (trapezoid f a b n) (* (/ (- b a) n) (+ (trapezoid-iter f a (/ (- b a) n) (- n 1)) (/ (f a) 2) (/ (f b) 2)))) (define (f2 x) (exp (- (* x x)))) 2. その後,次を「実行用ウインドウ」で実行しなさい (trapezoid f2 0 1 1000) ☆ 次は,課題に進んでください 28

29.

まず,関数 trapezoid-iter, trapezoid を定義している 29

30.

次に関数 f2 を定義している (define (f2 x) (exp (- (* x x)))) これは,数値積分したい関数 30

31.

これは, (trapezoid f2 0 1 1000) と書いて,a の値を 0 に, b の値を 1 に,h の値を 1000 に, f を f2 に設定しての実行 実行結果である 「#i0.7468240714991842」 が表示される 31

32.

(define (trapezoid-iter f a h k) (cond [(= k 1) (f (+ a h))] [else (+ (f (+ a (* h k))) (trapezoid-iter f a h (- k 1)))])) ;; trapezoid: number number number -> number ;; to compute the area under the graph of f between a and b (define (trapezoid f a b n) (* (/ (- b a) n) (+ (trapezoid-iter f a (/ (- b a) n) (- n 1)) (/ (f a) 2) (/ (f b) 2)))) 32

33.

例題3.台形則の計算の精度 • 例題2での台形則の計算の精度が,分割幅を小さ くするほど高精度になることを確かめる. • 但し, 2 -x 積分したい関数 f(x) = e 積分区間 分割数 0 から 1 10, 20, 40, 80, 160, 320, 640, 1280 の8通り 33

34.

例題3.「台形則の計算の精度」の手順 (1/2) 1. 次を「定義用ウインドウ」で,実行しなさい • 入力した後に,Execute ボタンを押す (define (trapezoid-iter f a h k) (cond [(= k 1) (f (+ a h))] [else (+ (f (+ a (* h k))) (trapezoid-iter f a h (- k 1)))])) これは,例題2と同じ ;; trapezoid: number number number -> number ;; to compute the area under the graph of f between a and b (define (trapezoid f a b n) (* (/ (- b a) n) (+ (trapezoid-iter f a (/ (- b a) n) (- n 1)) (/ (f a) 2) (/ (f b) 2)))) (define (f2 x) (exp (- (* x x)))) 34

35.

例題3.「台形則の計算の精度」の手順(2/2) 2. その後,次を「実行用ウインドウ」で実行 しなさい (trapezoid (trapezoid (trapezoid (trapezoid (trapezoid (trapezoid (trapezoid (trapezoid 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 10) 20) 40) 80) 160) 320) 640) 1280) ☆ 次は,課題に進んでください 35

36.

数値積分の精度 • 分割幅を小さくするほど高精度 36

37.

おわりに • 近似値を求める手法 • 数値微分,数値積分 • 十分に役に立つ • 「必ず誤差が生じる」ことを意識しながら使うことが 必要 37

38.

13-3 課題 38

39.

課題1 • f(x) = x cos x について,f'(0), f'(0.2), f'(0.4), f'(0.6) の値を報告しなさい • 関数 d/dx (授業の例題1)を使いなさい • h = 0.1 としなさい 39

40.

演習2 • f(x) = x cos x について,h = 0.1, 0.01, 0.001 と変えて関数 d/dx を実行し,結 果を比べなさい • h と誤差の関係が分かるグラフを作成せ よ 40

41.

演習3 • 台形則を使って,次を計算せよ log x 2 =  dx 0 1+ x 1 e • 計算結果を,以下と比較せよ DrScheme での (log 2) の実行結果 41

42.

演習4 • 演習3について,区間数 n を,n = 4, 8, 16, 32, 64, 128 と変えて計算を行い,刻み幅と誤 差の関係を調べよ • 区間数 n と誤差の関係を書いたグラフを作成せよ • この場合,おおよそ次の関係が成り立っているこ とを確認せよ c Sn − I  2 n (c は定数) 42

43.

演習4.関数のグラフ • 関数 d/dx (例題1)を使って,関数 のグラフをグラフィックスで表示させ なさい • プログラムは次ページ以降にある.この プログラムを実行させ,実行結果を報告 しなさい • 但し,実行結果の報告では,必ず 関数 f2 を別のものに書き換えて実行しなさい 43

44.

;; d/dx: (number->number) number number -> number ;; inclination of the tangent ;; Example: (d/dx f 3 0.0001) ;; = ((- (f 3.0001) (f 2.9999)) 0.0002) (define (d/dx f x h) (/ (- (f (+ x h)) (f (- x h))) (* 2 h))) ;; samples: (number->number) number number number -> list of 'posn' structure (define (samples f a h k) (cond [(= k 1) (cons (make-posn (+ a h) (f (+ a h))) empty)] [else (cons (make-posn (+ a (* h k)) (f (+ a (* h k)))) (samples f a h (- k 1)))])) 44

45.

; window size (define WX 500) (define WY 500) ; grid color, axis color and (define GRID_COLOR 'blue) (define AXIS_COLOR 'red) (define GRAPH_COLOR 'red) ; draw-a-line (define (sessen x px py d) (+ (* d (- x px)) py)) (define (draw-a-sessen from to px py d x0 y0 sx sy) (draw-solid-line (make-posn (+ (* from sx) x0) (+ (* (sessen from px py d) sy) y0)) (make-posn (+ (* to sx) x0) (+ (* (sessen to px py d) sy) y0)) GRAPH_COLOR)) (define (draw-sessen f px x0 y0 sx sy) (draw-a-sessen (/ (- x0) sx) (/ (- WX x0) sx) px (f px) (d/dx f px 0.00001) x0 y0 sx sy)) 45

46.

; draw-polyline (define (draw-polyline a-poly) (cond [(empty? (rest a-poly)) true] [else (and (draw-solid-line (first a-poly) (first (rest a-poly)) GRAPH_COLOR) (draw-polyline (rest a-poly)))])) ; draw-h-lines (define (draw-h-lines start skip stop width) (cond [(>= start stop) (draw-solid-line (make-posn 0 stop) (make-posn width stop) GRID_COLOR)] [else (and (draw-solid-line (make-posn 0 start) (make-posn width start) GRID_COLOR) (draw-h-lines (+ start skip) skip stop width))])) 46

47.

; draw-v-lines (define (draw-v-lines start skip stop height) (cond [(>= start stop) (draw-solid-line (make-posn stop 0) (make-posn stop height) GRID_COLOR)] [else (and (draw-solid-line (make-posn start 0) (make-posn start height) GRID_COLOR) (draw-v-lines (+ start skip) skip stop height))])) ; (x0, y0) is the origin. sx and sy represent one grid size (define (draw-grid x0 y0 sx sy) (and (draw-v-lines (- x0 (* (quotient x0 sx) sx)) sx (+ (* (quotient (- WX x0) sx) sx) x0) WY) (draw-h-lines (- y0 (* (quotient y0 sy) sy)) sy (+ (* (quotient (- WY y0) sy) sy) y0) WX) (draw-solid-line (make-posn 0 y0) (make-posn WX y0) AXIS_COLOR) (draw-solid-line (make-posn x0 0) (make-posn x0 WY) AXIS_COLOR))) 47

48.

; (x0, y0) is the origin. sx and sy represent one grid size (define (scaling a-list x0 y0 sx sy) (cond [(empty? a-list) empty] [else (cons (make-posn (+ (* (posn-x (first a-list)) sx) x0) (+ (* (posn-y (first a-list)) sy) y0)) (scaling (rest a-list) x0 y0 sx sy))])) ; ; f2: number -> number (define (f2 x) (- (* x x) 2)) ; (X0, Y0) reprerents the origin of the graph. ; SX and SY represent the size of one grid (define X0 (/ WX 2)) (define Y0 (/ WY 2)) (define SX 50) (define SY 50) (define PX 0.5) (start WX WY) (draw-grid X0 Y0 SX SY) (draw-polyline (scaling (samples f2 (/ (- X0) SX) (/ 1 SX) WX) X0 Y0 SX SY)) (draw-sessen f2 PX X0 Y0 SX SY) 48