-
Notifications
You must be signed in to change notification settings - Fork 5
/
6.html
930 lines (842 loc) · 36.4 KB
/
6.html
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
<!doctype html>
<html lang="ja">
<head>
<meta charset="utf-8">
<title>プログラマの為の数学勉強会</title>
<!-- For reveal.js -->
<link rel="stylesheet" href="lib/reveal/css/reveal.css">
<link rel="stylesheet" href="lib/reveal/css/theme/night.css">
<link rel="stylesheet" href="lib/reveal/lib/css/ir_black.css">
<!-- For Graphics -->
<link rel="stylesheet" href="css/graphics.css">
<style>
.reveal .chapter-title {
margin-top: 3em;
}
.reveal {
font-size: 32px;
line-height: 1.4em;
}
.reveal .slides {
text-align: left;
}
.reveal section img {
border: none;
background: 0;
margin-left: 1em;
margin-right: 1em;
box-shadow: none;
}
.reveal strong {
color: yellow;
}
.reveal sup {
font-size: 40%;
}
.reveal table {
margin-top: 0.5em;
margin-bottom: 0.5em;
border: 2px solid lightblue;
}
.reveal pre {
font-size: 0.7em;
}
.reveal pre code {
max-height: 600px;
}
.reveal .note {
font-size: 50%;
}
.reveal .controls div.navigate-up,
.reveal .controls div.navigate-down {
display: none;
}
.reveal .block {
border: solid 2px;
position: relative;
border-radius: 8px;
margin-top: 0.5em;
margin-bottom: 0.5em;
padding: 1em 0.8em 1em 0.8em;
}
.reveal .block:after {
content: "";
display: block;
clear: both;
height: 1px;
overflow: hidden;
}
.reveal .answer {
color: #111111;
}
.reveal .block h4 {
position: absolute;
top: -0.5em;
margin: 0 auto;
background: #111111;
font-weight: bold;
}
</style>
<!-- Setup libraries for RequireJS-->
<script src="lib/require.js"></script>
<script>
requirejs.config({
baseUrl: "js",
paths: {
d3: "../lib/d3/d3.v3.min",
numeric: "../lib/numeric-1.2.6",
MathJax: "https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.1/MathJax.js?config=TeX-AMS_HTML"
},
shim: {
d3: { exports: "d3" },
numeric: { exports: "numeric" },
MathJax: { exports: "MathJax" }
}
});
</script>
<!-- Initialize MathJax -->
<script type="text/x-mathjax-config">
require(["MathJax"], function (MathJax){
MathJax.Hub.Register.StartupHook("AsciiMath Jax Config", function () {
var AM = MathJax.InputJax.AsciiMath.AM;
AM.symbols.push(
{input:"mathbi",tag:"mstyle",atname:"mathvariant",atval:"bold-italic",
output:"mathbi",tex:null,ttype:AM.TOKEN.UNARY}
);
});
MathJax.Hub.Config({
showProcessingMessages: false,
skipStartupTypeset: true,
tex2jax: {
inlineMath: [ ["\\(","\\)"] ],
displayMath: [ ["\\[","\\]"] ]
}
});
});
</script>
</head>
<body>
<div class="reveal">
<div class="slides">
<section style="text-align: center">
<h1> プログラマの為の<br>数学勉強会<br>第6回</h1>
<span>
(於)ワークスアプリケーションズ<br>
中村晃一<br>
2013年10月17日
</span>
</section>
<section>
<h2>謝辞</h2>
<p>
この会の企画・会場設備の提供をして頂きました<br>
㈱ ワークスアプリケーションズ様<br>
にこの場をお借りして御礼申し上げます。
</p>
</section>
<section>
<h2> この資料について </h2>
<p>
<ul>
<li> <a href="http://nineties.github.com/math-seminar">
http://nineties.github.com/math-seminar
</a>に置いてあります。 </li>
<li> SVGに対応したブラウザで見て下さい。主要なブラウザで古いバージョンでなければ大丈夫だと思います。</li>
<li> 内容の誤り、プログラムのバグは<a href="http://twitter.com/9_ties">@9_ties</a>かkoichi.nakamur AT gmail.comまでご連絡下さい。</li>
<li> サンプルプログラムはPythonで記述しています。 </li>
</ul>
</p>
</section>
<section>
<h2 class="chapter-title"> 数値積分法 </h2>
</section>
<section>
<h2> 数値積分法 </h2>
<p>
解析的に定積分を求める事が出来ない・難しい関数では数値計算によって定積分の値を求める事になります。
</p>
<div align="center"> <img width="600px" src="fig/riemann-integral2.png"> </div>
</section>
<section graphics="interpolation" style="font-size:90%">
<h2> 準備:ラグランジュ補間 </h2>
<p>
データ列\((x_0, y_0),\ \cdots\ ,(x_n, y_n)\)が与えられた時に,\(x \neq x_i\)における\(y\)の値の近似値を与える関数\(y=f(x)\)を得る事を<strong>補間</strong>と言います。
</p>
<p>
下のグラフのように\(n+1\)個の点全てを通る\(n\)次多項式(これは唯一つに定まる)によって間を埋める補間法を<strong>ラグランジュ補間</strong>と言います。
</p>
<svg></svg>
</section>
<section style="font-size:90%">
<p>
例えば\(x\)座標の異なる3つのデータ\((x_0, y_0),(x_1, y_1),(x_2, y_2)\)に対して
\[ \begin{aligned}
p_2(x) = &\frac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)}y_0 + \frac{(x-x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)}y_1\\
& + \frac{(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)}y_2
\end{aligned} \]
という多項式関数を考えると,これは
\[ p_2(x_0) = y_0,\ p_2(x_1) = y_1,\ p_2(x_2) = y_2 \cdots(1)\]
を満たします。そして3点を通る放物線は一意に定まる事から,これが(1)を満たす唯一の二次多項式関数という事になります。
</p>
</section>
<section style="font-size:90%">
<p>
同様にして,一般にデータ列\((x_0,y_0),\cdots,(x_n,y_n)\)に対するラグランジュ補間の公式は
\[ \color{yellow}{p_n(x) = \sum_{i=0}^n \frac{(x-x_0)\cdots(x-x_{i-1})(x-x_{i+1})\cdots(x-x_n)}{(x_i-x_0)\cdots(x_i-x_{i-1})(x_i-x_{i+1})\cdots(x_i-x_n)}y_i} \]
となります。
</p>
<p>
\(g(x) = (x-x_0)(x-x_2)\cdots(x-x_n)\)とおくと
\[ \color{yellow}{p_n(x) = \sum_{i=0}^n \frac{g(x)}{g'(x_i)}\frac{y_i}{x-x_i} } \]
と書くことも出来ます。後で使います。
</p>
</section>
<section graphics="interpolation2">
<h2> ラグランジュ補間法の誤差 </h2>
<p>
データ列\((x_0, f(x_0)),\cdots,(x_n, f(x_n))\)に対するラグランジュ補間多項式を\(p_n(x)\)とし,\(f(x)\)と\(p_n(x)\)の誤差を考えます。
</p>
<svg></svg>
</section>
<section style="font-size:80%">
<p>
\(f\)は必要なだけ微分可能であるとします。また\(x_i\)は全て異なるとします。
</p>
<p>
まず,ある\(x\)に対する\(f(x)\)と\(p_n(x)\)の誤差を「\((x-x_0)\cdots(x-x_n)\)の\(r\)倍」という形で表して、
\[ f(x) = p_n(x) + r(x-x_0)\cdots(x-x_n)\]
とおきます。ただし\(x \neq x_i\)とします。
</p>
<p>
ここで関数
\[ g(y) = f(y) - p_n(y) - r(y-x_0)\cdots(y-x_n) \]
を考えると定義より
\[g(x_0)=g(x_1)=\cdots=g(x_n) = g(x) = 0\]
となります。すると,ロルの定理(第３回資料)より
\[ g'(y) = 0\]
となる相異なる点が少なくとも\(n+1\)点存在します。再びロルの定理を用いれば
\[ g''(y) = 0\]
となる相異なる点が少なくとも\(n\)点存在します。
</p>
<p>
これを繰り返していくと,最終的にある\(c\)が\(x_0,\cdots,x_n,x\)を含む閉区間内に存在して
\[ \color{yellow}{g^{(n+1)}(c) = 0}\]
となります。
</p>
</section>
<section style="font-size:90%">
<p>
\[ g(y) = f(y) - p_n(y) - r(y-x_0)\cdots(y-x_n) \]
において\(p_n(y)\)は\(n\)次式なので\(n+1\)回微分すると0。
また
\[ r(y-x_0)\cdots(y-x_n) = ry^{n+1} + \cdots \]
なのでこれを\(n+1\)回微分すると
\[ r\times (n+1)! \]
従って
\[ \begin{aligned}
g^{(n+1)}(c) = 0\ &\Leftrightarrow\ f^{(n+1)}(c)-r\times (n+1)! = 0\\
&\Leftrightarrow\ \color{yellow}{r = \frac{f^{(n+1)}(c)}{(n+1)!}}
\end{aligned}
\]
と表されます。
</p>
</section>
<section>
<div class="block" style="border-color:pink;font-size:90%">
<h4 style="color:pink"> ラグランジュ補間 </h4>
<p>
\( x_i \)が全て異なるとき
\[ \small{p_n(x) = \sum_{i=0}^n \frac{(x-x_0)\cdots(x-x_{i-1})(x-x_{i+1})\cdots(x-x_n)}{(x_i-x_0)\cdots(x_i-x_{i-1})(x_i-x_{i+1})\cdots(x_i-x_n)}f(x_i)} \]
が\(n\)次のラグランジュ補間多項式となる。
</p>
<p>
この時\(f(x)\)と\(p_n(x)\)の誤差はある\(c\)が\(x_0,\cdots,x_n, x\)を含む閉区間内に存在して
\[ \small{\frac{f^{(n+1)}(c)}{(n+1)!}(x-x_0)(x-x_1)\cdots(x-x_n)} \]
と表される。
</p>
</div>
</section>
<section style="font-size:80%">
<div class="block" style="border-color:lightgreen">
<h4 style="color:lightgreen">例</h4>
<p>
\(\sin 0 = 0,\ \sin\frac{\pi}{6} = \frac{1}{2},\ \sin\frac{\pi}{2} = 1\)を利用して\(\sin x\)を補間します。
</p>
</div>
<p>
\[ p_2(x) = \frac{(x-0)(x-\frac{\pi}{2})}{(\frac{\pi}{6}-0)(\frac{\pi}{6}-\frac{\pi}{2})}\frac{1}{2} + \frac{(x-0)(x-\frac{\pi}{6})}{(\frac{\pi}{2}-0)(\frac{\pi}{2}-\frac{\pi}{6})} = -\frac{3}{\pi^2}x^2 + \frac{7}{2\pi} x \]
となります。
</p>
<div align="center"> <img width="450px" src="fig/interpolation.png"> </div>
</section>
<section>
<h2> ニュートン＝コーツ系の公式 </h2>
<p>
では,本題の数値積分法の手法を紹介します。
</p>
<p>
積分区間\([a,b]\)に等間隔に取った点で\(f(x)\)をラグランジュ補間し,それを積分して得られる公式を<strong>ニュートン=コーツの積分公式</strong>と言います。
</p>
</section>
<section style="font-size:80%">
<h2> 台形公式 </h2>
<p>
二点\((x_i, f_i),(x_{i+1}, f_{i+1}))\)をラグランジュ補間すると
\[ \begin{aligned}
f(x) &\approx \frac{x-x_{i+1}}{x_i-x_{i+1}}f_i + \frac{x-x_i}{x_{i+1}-x_i}f_{i+1}\\
&= -\frac{f_i}{\Delta x}(x-x_{i+1}) + \frac{f_{i+1}}{\Delta x}(x-x_i)
\end{aligned} \]
となります。
</p>
<div align="center"> <img width="450px" src="fig/trapezoidal.png"> </div>
</section>
<section style="font-size:80%">
<p>
これを積分すれば
\[ \begin{aligned}
\int_{x_i}^{x_{i+1}} f(x)\mathrm{d} x &\approx \left[ -\frac{f_i}{2\Delta}(x-x_{i+1})^2 + \frac{f_{i+1}}{2\Delta x}(x-x_i)^2\right]_{x_i}^{x_{i+1}} \\
&= \color{yellow}{\frac{f_i + f_{i+1}}{2}\Delta x}
\end{aligned} \]
と近似できます。
</p>
</section>
<section style="font-size:80%">
<p>
次に剰余項
\[ R = \int_{x_i}^{x_{i+1}}\frac{f''(c)}{2}(x-x_i)(x-x_{i+1})\mathrm{d}x \]
についてですが,\(f''(x)\)の\([x_i,x_{i+1}]\)での最大値・最小値を\(M,m\)とすれば,この区間で\((x-x_i)(x-x_{i+1})\leq 0\)だから
\[ \small{\frac{M}{2}(x-x_i)(x-x_{i+1})\leq \frac{f''(c)}{2}(x-x_i)(x-x_{i+1}) \leq \frac{m}{2}(x-x_i)(x-x_{i+1})} \]
よって,各辺を積分して
\[ -\frac{M}{12}\Delta x^3 \leq R \leq -\frac{m}{12}\Delta x^3\ \Leftrightarrow\ m\leq -\frac{12}{\Delta x^3}R \leq M \]
となります。ところで\(f''(x)\)が連続ならば,\(m\)から\(M\)までの全ての値を取るのだからある\(c' \in [x_i, x_{i+1}]\)が存在して
\[ f''(c') = -\frac{12}{\Delta x^3}R\ \Leftrightarrow\ R = -\frac{f''(c')}{12}\Delta x^3 \]
と表せます。
</p>
</section>
<section style="font-size:80%">
<p>
区間\([a,b]\)を\(n\)等分した各区間に先ほどの公式を利用して加えれば
\[ \begin{aligned}
\int_a^b f(x)\mathrm{d}x &\approx \frac{f_0+f_1}{2}\Delta x + \frac{f_1+f_2}{2}\Delta x + \cdots + \frac{f_{n-1}+f_n}{2}\Delta x \\
&= \color{yellow}{\frac{\Delta x}{2}(f_0+2f_1+2f_2+\cdots + 2f_{n-1} + f_n)}
\end{aligned} \]
という公式を得ます。
</p>
<p>
また,区間内での\(|f''(x)|\)の上限を\(M\)とすれば各区間での誤差の絶対値が\( \frac{M}{12}\Delta x^3 \)以下となるので,総誤差は
\[ \frac{M}{12}\Delta x^3\times n = \frac{M(b-a)}{12}\Delta x^2 \]
以下となります。
</p>
<div align="center"> <img width="450px" src="fig/trapezoidal2.png"> </div>
</section>
<section>
<div class="block" style="border-color:pink;font-size:90%">
<h4 style="color:pink"> 台形公式</h4>
<p>
\( \Delta x = \frac{b-a}{n},\ x_i = a+i\Delta x,\ f_i=f(x_i) \)とおき
\[ \int_a^b f(x)\mathrm{d} x \approx \frac{\Delta x}{2}(f_0+2f_1+2f_2+\cdots+2f_{n-1}+f_n) \]
と近似する公式を台形公式という。
</p>
<p>
この時の絶対誤差は\(|f''(x)|\)の\([a,b]\)での上限を\(M\)とすれば
\[ \frac{M(b-a)}{12}\Delta x^2 \]
以下となる。
</p>
</div>
</section>
<section style="font-size:80%">
<div class="block" style="border-color:lightgreen">
<h4 style="color:lightgreen">例</h4>
<p> \(\int_0^1\frac{1}{1+x^2}\mathrm{d} x\)を計算してみます。厳密値は\(\frac{\pi}{4}=0.785398\cdots\)です。</p>
</div>
<pre><code class="python" style="max-height:400px">>>> a = 0.0
>>> b = 1.0
>>> n = 10
>>> dx = (b-a)/n
>>> def f(x):
... return 1.0/(1+x*x)
...
>>> s = f(a) + f(b)
>>> for i in range(1,n):
... s += 2*f(a+i*dx)
...
>>> s *= dx/2
>>> s
0.7849814972267897
</code></pre>
<p>
</p>
</section>
<section style="font-size:80%">
<div class="block" style="border-color:lightgreen">
<h4 style="color:lightgreen">例続き</h4>
<p> \(\Delta x\)を変えながら誤差の変化を見てみましょう。 </p>
</div>
<p>
理論的には\(\Delta x\)を\(\frac{1}{10}\)倍する毎に,誤差は\(\frac{1}{100}\)倍程度に減るはずです。
</p>
<img align="right" width="450px" src="fig/trapezoidal-error.png">
<p>
実際には右の実験結果のように\(\Delta x\)を小さくしすぎるとこの通りになりません。これは演算回数が増加する事によって累積誤差が増加する為です。精度をさらに向上させる為には、別の手法を利用する必要があります。
</p>
</section>
<section>
<h2> シンプソン公式 </h2>
<p>
台形公式と同様にして,２次のラグランジュ補間を行って積分する事でより高精度の公式を得る事が出来ます。
</p>
<div align="center"> <img width="600px" src="fig/simpson.png"> </div>
</section>
<section style="font-size:90%">
<div class="block" style="border-color:pink;font-size:90%">
<h4 style="color:pink"> シンプソン公式</h4>
<p>
\(n\)は偶数とする。\( \Delta x = \frac{b-a}{n},\ x_i = a+i\Delta x,\ f_i=f(x_i) \)とおき
\[ \begin{aligned}
\int_a^b f(x)\mathrm{d} x \approx &\frac{\Delta x}{3}(f_0+4f_1+2f_2+4f_3+2f_4+\\
&\cdots+2f_{n-2} +4f_{n-1}+f_n)
\end{aligned} \]
と近似する公式をシンプソン公式という。
</p>
<p>
この時の絶対誤差は\(|f^{(4)}(x)|\)の\([a,b]\)での上限を\(M\)とすれば
\[ \frac{M(b-a)}{180}\Delta x^4 \]
以下となる。
</p>
</div>
<p>
詳しい導出は省略しますが,台形公式と同様です。
</p>
</section>
<section style="font-size:80%">
<div class="block" style="border-color:lightgreen">
<h4 style="color:lightgreen">例</h4>
<p> \(\int_0^1\frac{1}{1+x^2}\mathrm{d} x\)を計算してみます。厳密値は\(\frac{\pi}{4}=0.785398\cdots\)です。</p>
</div>
<pre><code class="python" style="max-height:400px">>>> a = 0.0
>>> b = 1.0
>>> n = 10
>>> dx = (b-a)/n
>>> def f(x):
... return 1.0/(1+x*x)
...
>>> s = f(a) + f(b)
>>> for i in range(1,n,2):
... s += 4*f(a+i*dx)
...
>>> for i in range(2,n,2):
... s += 2*f(a+i*dx)
...
>>> s *= dx/3
>>> s
0.7853981534848038
</code></pre>
<p>
</p>
</section>
<section>
<div class="block" style="border-color:lightgreen">
<h4 style="color:lightgreen">例続き</h4>
<p> \(\Delta x\)を変えながら誤差の変化を見てみましょう。 </p>
</div>
<div align="center"> <img width="450px" src="fig/simpson-error.png"> </div>
</section>
<section>
<h2> ガウス求積法 </h2>
<p>
<strong>ガウス求積法</strong>と呼ばれる,ニュートン=コーツ系の公式とは異なる手法を紹介します。
</p>
</section>
<section style="font-size:80%">
<h2> ルジャンドル多項式 </h2>
<p>
\[ P_2(x) = \frac{3}{2}x^2-\frac{1}{2} \]
という多項式は特別な性質を持っていて,任意の\(a,b\)に対して
\[ \int_{-1}^1 (ax+b)P_2(x)\mathrm{d} x = 0 \]
となります(計算してみましょう)。
</p>
<p>
同様に
\[ P_3(x) = \frac{5}{2}x^3 - \frac{3}{2}x \]
という多項式は任意の\(a,b,c\)に対して
\[ \int_{-1}^1 (ax^2+bx+c)P_3(x)\mathrm{d} x = 0\]
を満たします。
</p>
</section>
<section style="font0-size:80%">
<p>
一般に任意の\(n-1\)次以下の多項式\(f(x)\)に対して
\[ \int_{-1}^1 f(x)P_n(x)\mathrm{d} x = 0\]
を満たす\(n\)次多項式\(P_n(x)\)は(定数倍の不定性がまだ残っていますが)<strong>ルジャンドル多項式</strong>と呼ばれます。
</p>
<p>
詳しい導出は省略しますが,
\[
P_1(x) = x,\ P_2(x) = \frac{3}{2}x^2-\frac{1}{2},\ P_3(x) = \frac{5}{2}x^3-\frac{3}{2}x \]
などとなり,一般には漸化式
\[ (n+1)P_{n+1}(x) = (2n+1)xP_n(x)-nP_{n-1}(x) \]
で得る事が出来ます。
</p>
</section>
<section style="font-size:90%">
<p>
\(P_n(x)\)をルジャンドル多項式,\(P_n(x)=0\)の解を\(\alpha_1,\cdots,\alpha_n\)とします。
</p>
<p>
ここで\(f(x)\)を任意の\(2n-1\)次多項式とすると,割り算を行って
\[ f(x) = q(x)P_n(x) + r(x) \qquad \text{($q(x),r(x)$は$n-1$次以下)} \]
と表す事が出来ますが,ルジャンドル多項式の性質より
\[ \int_{-1}^1 f(x)\mathrm{d} x = \int_{-1}^1 r(x)\mathrm{d} x \]
と被積分関数の次数を下げる事が出来ます。
</p>
</section>
<section style="font-size:85%">
<p>
ここで\(r(x)\)を\(x = \alpha_1,\cdots,\alpha_n\)での値を使って補間すると,\(r(\alpha_i)=f(\alpha_i)\)である事に注意すれば
\[ r(x) = \sum_i\frac{P_n(x)}{P'_n(\alpha_i)}\frac{f(\alpha_i)}{x-\alpha_i} \]
となります。\(r(x)\)は\(n-1\)次以下なので,これは厳密に成立します。
</p>
<p>
すると
\[ \int_{-1}^1 f(x)\mathrm{d} x = \int_{-1}^1 r(x)\mathrm{d} x = \sum_i f(\alpha_i)\int_{-1}^{1}\frac{P_n(x)}{P'_n(\alpha_i)}\frac{\mathrm{d}x}{x-\alpha_i} \]
となります。右辺の積分は定数なので\(w_i\)とおけば
\[ \int_{-1}^1 f(x)\mathrm{d} x = \sum_i f(\alpha_i)w_i \]
となります。
</p>
<p>
つまり,<strong>\(2n-1\)次多項式の積分値をたった\(n\)点での値を使って厳密に求める事が出来る</strong>という事になります。
</p>
</section>
<section style="font-size:80%">
<div class="block" style="border-color:pink;font-size:90%">
<h4 style="color:pink"> ガウス求積法</h4>
<p>
\[ \int_{-1}^1 f(x)\mathrm{d} x \approx \sum_{i=1}^n f(\alpha_i)w_i \]
但し\(\alpha_i\)はルジャンドル多項式\(P_n(x)\)の根であり,
\[ w_i = \int_{-1}^{1}\frac{P_n(x)}{P'_n(\alpha_i)}\frac{\mathrm{d}x}{x-\alpha_i} \]
となる。各\(n\)に対する\((\alpha_i, w_i)\)の値は
\[ \small{\begin{aligned}
n = 1 &: \left(0, 2\right)\quad n = 2 : \left(\pm\sqrt{\frac13}, 1\right)\quad n = 3 : \left(0, \frac89\right),\ \left(\pm\sqrt{\frac35},\frac59\right)\\
n = 4 &: \left(\pm\sqrt{\frac{3-2\sqrt{6/5}}7}, \frac{18+\sqrt{30}}{36}\right),\ \left(\pm\sqrt{\frac{3+2\sqrt{6/5}}7}, \frac{18-\sqrt{30}}{36}\right)
\end{aligned}} \]
などとなる。
</p>
<p>
\(f(x)\)を\(2n-1\)次多項式で補間した場合と同じ精度となる。
</p>
</div>
</section>
<section style="font-size:80%">
<div class="block" style="border-color:lightgreen">
<h4 style="color:lightgreen">例</h4>
<p> \(\int_0^1\frac{1}{1+x^2}\mathrm{d} x\)を計算してみます。厳密値は\(\frac{\pi}{4}=0.785398\cdots\)です。</p>
</div>
<p>
\(x = (t + 1)/2\)と置換して,
\[ \int_0^1\frac{1}{1+x^2}\mathrm{d} x = \int_{-1}^1 \frac{1}{1+(t+1)^2/4}\frac{1}{2}\mathrm{d} t = \int_{-1}^1 \frac{2}{4+(t+1)^2}\mathrm{d} t\]
としてから公式を利用すれば良いです。
</p>
<pre><code class="python" style="max-height:400px">>>> def f(x):
... return 2.0/(4+(x+1)**2)
...
>>> a = 0.577350269189626 # 1/sqrt{3}
>>> f(-a) + f(a) # n = 2の公式
0.7868852459016393
>>> b = 0.774596669241483 # sqrt{3/5}
>>> f(-b)*5.0/9 + f(0)*8.0/9 + f(b)*5.0/9 # n = 3の公式
0.785267034990792
</code></pre>
</section>
<section style="font-size:80%">
<h2> ロンバーグ積分法 </h2>
<p>
これまでとは全く別の考え方を用いる手法を紹介します。
</p>
<p>
定積分\(\int_a^b f(x)\mathrm{d}x \)を分割数\(n=2^i\)の台形公式で求めた値を\(T^0_i\)と表す事にします。
\[ T^0_i = \frac{h_i}{2}(f_0+2f_1+2f_2+\cdots+f_{2^i})\qquad \left(h_i = \frac{b-a}{2^i}\right)\]
すると
\[ T^0_0,\ T^0_1,\ T^0_2, T^0_3, \cdots, \]
という数列が得られます。理論的にはこの数列で\(i\rightarrow\infty\)とすれば厳密な値に収束しますが,実際には累積誤差が生じるのでそうはいきません。
</p>
</section>
<section style="font-size:80%">
<p>
しかし,\(T^0_n\)の値を直に計算しなくても,\(T^0_0,\ T^0_1,\ T^0_2,\ \cdots\)の値からそれを予測する事が可能です。この考え方を<strong>外挿(補外)</strong>と言います。
</p>
<div align="center"> <img width="450px" src="fig/romberg.png"> </div>
</section>
<section style="font-size:70%">
<p>
難しくなるので詳しい導出はしませんが,台形近似の剰余項は\(i\)に依存しない定数\(\alpha_1,\alpha_2,\cdots\)を用いて
\[ \int_a^bf(x)\mathrm{d} x = T^0_i + \alpha_1h_i^2 + \alpha_2h_i^4 + \alpha_3h_i^6 + \cdots \]
と表す事が出来ます(<strong>オイラー=マクローリンの公式</strong>と言います)。
</p>
<p>
ここで\(h_i = \frac{h_{i-1}}{2}\)であるので,
\[ \int_a^bf(x)\mathrm{d} x = T^0_i + \alpha_1\frac{h_{i-1}^2}{4} + \alpha_2\frac{h_{i-1}^4}{4^2} + \alpha_3\frac{h_{i-1}^6}{4^3} + \cdots (1)\]
となりますが,一方
\[ \int_a^bf(x)\mathrm{d} x = T^0_{i-1} + \alpha_1h_{i-1}^2 + \alpha_2h_{i-1}^4 + \alpha_3h_{i-1}^6 + \cdots (2)\]
でもあるので,(1)を4倍して(2)を引けば
\[ 3\int_a^b f(x)\mathrm{d} x = 4T^0_i - T^0_{i-1} + \mathcal{O}(h_{i-1}^4)\ \Leftrightarrow\ \int_a^b f(x)\mathrm{d} x = \frac{4T^0_i - T^0_{i-1}}{3} + \mathrm{O}(h_{i-1}^4) \]
となって剰余項から\(\mathcal{O}(h_{i-1}^2)\)の項を消してしまう事が出来ます。
</p>
<p>
同様にして,剰余項が\(\mathcal{O}(h^4)\)の公式を2つ組み合わせれば\(\mathcal{O}(h^4)\)の項も消してしまう事が出来,以上を繰り返す事で精度を高める事が出来ます。
</p>
</section>
<section>
<div class="block" style="border-color:pink;font-size:90%">
<h4 style="color:pink"> ロンバーグ積分法</h4>
<p>
\(T^0_i\)を分割数\(2^i\)の台形公式で求めた値とする。
\[ T^j_i = \frac{4^jT_i^{j-1}-T_{i-1}^{j-1}}{4^j-1}\ (i = 1, 2, \cdots, n\ j = 1, 2, \cdots, i) \]
によって外挿を行い,\(T^n_n\)を求める。
</p>
</div>
<p style="font-size:80%"> \(T^0_i\)は実際には漸化式でもっと効率よく計算出来ますが,分り易さの為に説明を省きます。
</p>
<p style="font-size:80%">
また,今回利用した補外法は一般に<strong>リチャードソン補外</strong>と呼ばれます。
</p>
</section>
<section style="font-size:80%">
<div class="block" style="border-color:lightgreen">
<h4 style="color:lightgreen">例</h4>
<p> \(\int_0^1\frac{1}{1+x^2}\mathrm{d} x\)を計算してみます。厳密値は\(\frac{\pi}{4}=0.785398\cdots\)です。</p>
</div>
<pre><code class="python" style="max-height:400px">from numpy import *
def f(x):
return 1.0/(1+x*x)
a = 0.0
b = 1.0
N = 5 # ここを変えてみて下さい
T = zeros( (N+1, N+1) ) # N+1 x N+1の２次元配列。
T[0,0] = (f(a) + f(b))*(b-a)/2
for i in range(1, N+1):
# T[i,0]は台形公式で計算
n = 2**i
dx = (b-a)/n
s = f(a) + f(b)
for j in range(1,n):
s += 2*f(a + j*dx)
T[i,0] = s*dx/2
# 外挿
for j in range(1, i+1):
T[i,j] = T[i, j-1] + (T[i, j-1] - T[i-1,j-1])/(4**j-1)
print T[N,N]
</code></pre>
</section>
<section>
<p>
ロンバーグ積分の精度は以下のようになり,分割数\(128\)では浮動小数点数で表せる範囲での誤差\(0\)を達成しています。
</p>
<div align="center"> <img width="450px" src="fig/romberg-error.png"> </div>
</section>
<section>
<h2 class="chapter-title"> 積分法の応用 </h2>
</section>
<section style="font-size:80%">
<h2> 面積 </h2>
<p>
以前に説明したように,以下の領域の符号付き面積が
\[ \int_a^b f(x)\mathrm{d} x \]
によって得られます。
</p>
<div align="center"> <img width="400px" src="fig/riemann-integral2.png"> </div>
</section>
<section style="font-size:80%">
<p>
一般に,下図のように２曲線で囲まれた領域の面積は
\[ \int_a^b|f(x)-g(x)|\mathrm{d} x \]
となります。
</p>
<div align="center"> <img width="400px" src="fig/riemann-integral3.png"> </div>
<p class="fragment">
普段図形の面積を計算する場面なんてないと考える人もいるかもしれませんが、少くとも<strong>確率計算が面積計算に帰着する</strong>為、統計学の計算において必要となります。<span style="font-size:80%">(変数の数によっては面積ではなくて体積などとなりますが,区別せずに面積と呼ぶことにします。一般的には測度と呼ばれます。)</span>
</p>
</section>
<section style="font-size:80%">
<p>
例えば,\([0, 1]\)から実数\(x,y\)を無作為に選んだ場合
\[ y \geq x^2 \]
となる確率は
\[ \frac{\text{(色を塗った部分の面積)}}{\text{(正方形の面積)}} = \frac{\int_0^1(1-x^2)\mathrm{d}x}{1} = \frac{2}{3} = 0.6666\cdots\]
と計算されます。
</p>
<img align="right" width="300px" src="fig/parabola.png">
<pre><code class="python" style="max-height:400px">
>>> from random import random
>>> N = 1000
>>> count = 0
>>> for i in xrange(N):
... x = random()
... y = random()
... if y >= x*x:
... count += 1
...
>>> count
665
</code></pre>
<p>
上のコードはこの結果を乱数を用いた実験で確認するものです。また、上の様なコードで確率を近似的に求める事で,逆に右図の面積の近似値を求める事が出来ます。これは<strong>モンテカルロ法</strong>と呼ばれる手法です。統計の回に説明する予定です。
</p>
</section>
<section style="font-size:80%">
<div class="block" style="border-color:lightgreen">
<h4 style="color:lightgreen">例:ビュフォンの針</h4>
間隔\(h\)で平行線がいくつも引いてある所に,長さ\(\ell\)の針を全く無作為に落としたとき,平行線と交差する確率を求めて下さい。\(\ell < h\)とします。
</div>
<p style="font-size:70%">
針の中心と一番近い平行線の距離を\(d\)とし,図の角度を\(\theta\)とすれば\( 0\leqq d \leqq \frac{h}{2},\quad 0\leqq \theta \leqq \frac{\pi}{2} \)
の範囲の実数が無作為に選ばれる事になります。
</p>
<p style="font-size:70%">
ところで,平行線と針が交わる条件は
\[ d \leqq \frac{\ell}{2}\cos\theta \]
であるので,右下図の面積の比率を求めれば
\[ \frac{1}{\frac{h}{2}\frac{\pi}{2}}\int_0^\frac{\pi}{2}\frac{\ell}{2}\cos x\mathrm{d} x = \frac{2\ell}{h\pi} \]
が求める確率となります。
</p>
<div align="center"><img width="200px" src="fig/buffons-needle.png">
<img width="300px" src="fig/buffons-needle2.png">
</div>
</section>
<section style="font-size:80%">
<div class="block" style="border-color:lightgreen">
<h4 style="color:lightgreen">例:確率分布の求積</h4>
<p>
得点分布が正規分布で近似できるとき,偏差値が\(a\)以上\(b\)以下の人の割合は
\[ \int_a^b \frac{1}{\sqrt{2\pi \times 100}}\exp\left\{-\frac{(x-50)^2}{2\times 100}\right\}\mathrm{d} x \]
となります。では偏差値が40以上60以下の人の割合を求めて下さい。
</p>
</div>
<p>
以下の面積を求めよという問題です。このようにして,確率分布に関する確率計算は面積計算の問題に帰着します。
詳しくは統計の回にやります。
</p>
<div align="center"> <img width="400px" src="fig/normal-distr.png"> </div>
</section>
<section>
<p>
この積分は解析的に求める事が出来ないので数値積分を行います。シンプソン法を利用すれば以下のようになり,約\(68.27\)%という事が判ります。
</p>
<pre><code class="python" style="max-height:400px">
>>> from math import pi, sqrt, exp
>>> def f(x):
... return 1/(sqrt(2*pi*100))*exp(-(x-50)**2/200)
...
>>> a = 40
>>> b = 60
>>> dx = 1.0
>>> s = 0.0
>>> for i in range(0, 20, 2):
... s += f(a+i*dx) + 4*f(a+(i+1)*dx) + f(a+(i+2)*dx)
...
>>> s *= dx/3
>>> s
0.6826900317769402
</code></pre>
</section>
<section>
<h2> 今回はここで終わります。 </h2>
<p>
重積分はヤコビアンの計算など線型代数の知識が必要となるので後の回にやります。また,微分方程式の厳密解法に関しても線形空間の知識や固有値に関する知識が必要となりますので後の回に回します。ということで微積分法はここで一旦終了します。
</p>
<p>
次回から<strong>線型代数</strong>をやっていきます。初回は<strong>ベクトル・行列・行列式</strong>の基礎知識,<strong>連立一次方程式</strong>の解法,<strong>ガウス消去法</strong>などを紹介します。
</p>
</section>
</div>
</div>
<script src="lib/reveal/lib/js/head.min.js"></script>
<script src="lib/reveal/js/reveal.js"></script>
<script>
Reveal.initialize({
width: 960,
height: 640,
controls: true,
progress: false,
history: true,
overview: false,
touch: true,
center: false,
rollingLinks: false,
transition: "page",
transitionSpeed: "default",
// When scale != 1, positions of mouse events will be incorrect.
minScale: 1.0,
maxScale: 1.0,
dependencies: [
{ src: "lib/reveal/lib/js/classList.js", condition: function() { return !document.body.classList; } },
{ src: "lib/reveal/plugin/highlight/highlight.js", async: true, callback: function() { hljs.initHighlightingOnLoad(); } },
{ src: "lib/reveal/plugin/zoom-js/zoom.js", async: true, condition: function() { return !!document.body.classList; } },
{ src: "lib/reveal/plugin/notes/notes.js", async: true, condition: function() { return !!document.body.classList; } }
]
});
// register event listeners
require(["MathJax"], function (MathJax) {
// Delay typesetting of slides
function typeset (idx) {
for (var i = idx - 2, n = idx + 2; i <= n; i++) {
var slide = Reveal.getSlide(i);
if (!slide) continue;
if (!slide.typeset) {
MathJax.Hub.Typeset(slide);
slide.typeset = true;
}
}
}
function initializeGraphics (idx) {
for (var i = idx - 2, n = idx + 2; i <= n; i++) {
var slide = Reveal.getSlide(i);
if (!slide) continue;
var graphics = slide.getAttribute("graphics");
if (graphics && !slide.initialized) {
slide.initialized = true;
(function () {
var p = slide;
require([graphics], function(g) {
if (g.initialize) g.initialize(p);
});
})();
}
}
}
function start (slide) {
var graphics = slide.getAttribute("graphics");
if (graphics) {
require([graphics], function(g) { if (g.start) g.start(slide); });
}
}
function stop (slide) {
var graphics = slide.getAttribute("graphics");
if (graphics) {
require([graphics], function(g) { if (g.stop) g.stop(slide); });
}
}
function simpleEvent (type) {
var event = document.createEvent("HTMLEvents");
event.initEvent(type, true, true);
return event;
}
Reveal.addEventListener("slidechanged", function (event) {
typeset(event.indexh);
initializeGraphics(event.indexh);
start(event.currentSlide);
stop(event.previousSlide);
});
Reveal.addEventListener("fragmentshown", function (event) {
var slide = Reveal.getCurrentSlide();
var graphics = slide.getAttribute("graphics");
if (graphics) {
require([graphics], function(g) { if (g.proceed) g.proceed(slide); });
}
});
typeset(Reveal.getIndices().h);
initializeGraphics(Reveal.getIndices().h);
start(Reveal.getCurrentSlide());
});
</script>
</body>
</html>