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

sp-13. 数値微分と数値積分 （Scheme プログラミング） URL: https://www.kkaneko.jp/pro/scheme/index.html 金子邦彦 1

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

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

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

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

接線の傾き 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

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

区間 [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

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

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

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

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

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

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

13-2 パソコン演習 14

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

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

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

「例題１．接線の傾き」の手順 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) ☆ 次は，例題２に進んでください 18

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

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

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

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

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

(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

(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

(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

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

「例題２．台形則による数値積分」の手順 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

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

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

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

(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

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

例題３．「台形則の計算の精度」の手順 (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)))])) これは，例題２と同じ ;; 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

例題３．「台形則の計算の精度」の手順（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

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

13-3 課題 38

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

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

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

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

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

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

; 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

; 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

; 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

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