-
Notifications
You must be signed in to change notification settings - Fork 5
/
5.html
1126 lines (1017 loc) · 41.8 KB
/
5.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
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
<!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>第5回</h1>
<span>
(於)ワークスアプリケーションズ<br>
中村晃一<br>
2013年10月10日
</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>
関数\(y = f(x)\)に対して,変数\(x,y\),導関数\(y',y'',\cdots\)からなる方程式を<strong>微分方程式</strong>と言います。多変数の場合は<strong>偏微分方程式</strong>となります。
</p>
<p>
例えば
\[ \frac{\mathrm{d}y}{\mathrm{d}x} = -2y \]
であるとか
\[ \frac{\partial u}{\partial t}=k\frac{\partial^2 u}{\partial t^2} \]
などが微分方程式・偏微分方程式の例となります。
</p>
</section>
<section style="font-size:80%">
<p>
一般に微分方程式の解は<strong>任意定数</strong>を含みます。
</p>
<p>
例えば
\[ \frac{\mathrm{d}y}{\mathrm{d}x}=-2y \qquad\cdots(1)\]
という微分方程式の解は\(y=\color{yellow}{C}\exp(-2x)\quad\text{($C$は定数)}\)と表される関数全てが解となります(後述)。
</p>
<div align="center"> <img width="400px" src="fig/ode-fig1.png"> </div>
</section>
<section>
<p>
ここでは,微分方程式の解のうち与えられた条件(<strong>初期条件・境界条件</strong>)を満たす解を近似的に求める手法を紹介します。
下図の様に独立変数を離散化して,与えられた初期値から出発して次々に関数の値を求めていく事になります。
</p>
<div align="center"> <img width="500px" src="fig/numerical-method.png"> </div>
<p style="font-size:80%">
以後,図のように変数\(x\)の刻み\(\Delta x\)を定数にした場合を考えます。精度や計算速度を改善する為に\(\Delta x\)を可変にする手法もありますが、この勉強会では扱いません。
</p>
</section>
<section>
<h2> 一階常微分方程式 </h2>
<p>
まず一変数関数の微分方程式(<strong>常微分方程式</strong>)のうち,一階導関数のみを含む
\[ \frac{\mathrm{d} f}{\mathrm{d} x} = g(x,f) \]
という形の微分方程式について考えます。
</p>
</section>
<section style="font-size:90%">
<h2> 陽的オイラー法 </h2>
<p>
\[ \frac{\mathrm{d} f}{\mathrm{d} x} = g(x,f) \]
の左辺を前進差分法で近似し,右辺に<strong>\(x=x_i\)時点の値</strong>を使うと
\[ \frac{f_{i+1}-f_i}{\Delta x}\approx g(x_i,f_i) \]
と近似出来ます。これを変形した
\[ \color{yellow}{f_{i+1} = f_i + g(x_i,f_i)\Delta x}\]
という漸化式で次の点の値を計算する方法を<strong>陽的オイラー法</strong>と言います。
</p>
</section>
<section style="font-size:90%">
<h2> 陰的オイラー法 </h2>
<p>
\[ \frac{\mathrm{d} f}{\mathrm{d} x} = g(x,f) \]
の左辺を前進差分法で近似し,右辺に<strong>\(x=x_{i+1}\)時点の値</strong>を使うと
\[ \frac{f_{i+1}-f_i}{\Delta x}\approx g(x_{i+1},f_{i+1}) \]
と近似出来ます。これを変形した
\[ \color{yellow}{f_{i+1} = f_i + g(x_{i+1},f_{i+1})\Delta x}\]
という漸化式で次の点の値を計算する方法を<strong>陰的オイラー法</strong>と言います。
</p>
<p>
陰的解法では,右辺に\(f_{i+1}\)が含まれる為そのままの形では使えません。\(f_{i+1}\)について解くか,方程式の数値解法を利用する事になります。
</p>
</section>
<section style="font-size:80%">
<div class="block" style="border-color:lightgreen">
<h4 style="color:lightgreen">例</h4>
微分方程式の初期値問題
\[ \frac{\mathrm{d} x}{\mathrm{d}t} = x,\ x(0) = 1\]
を2つのオイラー法で解きましょう。
</div>
</section>
<section style="font-size:80%">
<p>
【陽的オイラー法】<br>
漸化式は
\[ \frac{x_{i+1}-x_i}{\Delta t} = x_i\ \Leftrightarrow x_{i+1} = (1+\Delta t) x_i \]
となります。
</p>
<pre><code class="python" style="max-height:400px">>>> x = 1 # x(0) = 1
>>> dt = 0.1
>>> for i in range(10):
... t = i*dt
... print "{0} {1}".format(t, x)
... x = (1 + dt) * x
...
0.0 1
0.1 1.1
0.2 1.21
0.3 1.331
0.4 1.4641
0.5 1.61051
0.6 1.771561
0.7 1.9487171
0.8 2.14358881
0.9 2.357947691
</code></pre>
</section>
<section style="font-size:80%">
<p>
【陰的オイラー法】<br>
漸化式は
\[ \frac{x_{i+1}-x_i}{\Delta t} = x_{i+1}\ \Leftrightarrow x_{i+1} = \frac{x_i}{1-\Delta t} \]
となります。
<pre><code class="python" style="max-height:400px">>>> x = 1 # x(0) = 1
>>> dt = 0.1
>>> for i in range(10):
... t = i*dt
... print "{0} {1}".format(t, x)
... x = x/(1-dt)
...
0.0 1
0.1 1.11111111111
0.2 1.23456790123
0.3 1.37174211248
0.4 1.52415790276
0.5 1.69350878084
0.6 1.88167642316
0.7 2.09075158129
0.8 2.32305731254
0.9 2.58117479171
</code></pre>
</section>
<section style="font-size:80%">
<p>
ところで厳密解は\(x = \exp t\)となります(代入して確かめて下さい)。
</p>
<p>
それぞれの解をプロットしてみると下のようになります。このように陽的解法は厳密解のカーブの外側に,陰的解法は内側にずれていく傾向があります。
</p>
<div align="center"> <img width="500px" src="fig/ode-fig2.png"> </div>
</section>
<section>
<h2> オイラー法の精度 </h2>
<p> 陽的オイラー法
\[ f_{i+1} \approx f_i + g(x_i,f_i)\Delta x\]
の左辺を\(x_i\)の周りにテイラー展開すると
\[ \begin{aligned}
f_{i+1} &= f_i + \frac{\mathrm{d} f}{\mathrm{d} x}(x_i)\Delta x + \mathcal{O}(\Delta x^2) \\
&= f_i + g(x_i, f_i)\Delta x + \mathcal{O}(\Delta x^2)
\end{aligned}
\]
となるので,この近似は\(\mathcal{O}(\Delta x)\)の項まで一致しています。このことを陽的オイラー法の精度は1次であると言います。
</p>
<p>
陰的オイラー法の精度も同様に1次(練習問題)となります。
</p>
</section>
<section>
<div class="block" style="border-color:pink;font-size:80%">
<h4 style="color:pink"> 陽的・陰的オイラー法 </h4>
<p>微分方程式
\[ \frac{\mathrm{d} f}{\mathrm{d} x} = g(x,f) \]
の各点\(x_i\)における値を
\[ f_{i+1} = f_i + g(x_i,f_i)\Delta x\]
によって求める方法を<strong>陽的オイラー法</strong>
</p>
<p>
\[ f_{i+1} = f_i + g(x_{i+1},f_{i+1})\Delta x\]
によって求める方法を<strong>陰的オイラー法</strong>
という。
</p>
<p>
いずれも精度は1次である。
</p>
</div>
</section>
<section style="font-size:70%">
<h2> 改良オイラー法 </h2>
<p>
今度は陽的オイラー法と陰的オイラー法の更新式の平均をとって
\[ f_{i+1} = f_i + \frac{g(x_i, f_i) + g(x_{i+1},f_{i+1})}{2}\Delta x\]
としてみます。
</p>
<p>
しかし,右辺にある\(f_{i+1}\)はやはり面倒なので,陽的オイラー法
\[ f_{i+1} \approx f_i + g(x_i, f_i)\Delta x \]
で得られる近似値を代わりに使えば,
\[
f_{i+1} = f_i + \frac{g(x_i, f_i) + g(x_i + \Delta x, f_i + g(x_i, f_i)\Delta x)}{2}\Delta x \]
つまり
\[ \color{yellow}{\left\{\begin{array}{l}
k_1 = g(x_i, f_i) \\
k_2 = g(x_i + \Delta x, f_i + k_1\Delta x) \\
f_{i+1} = f_i + \frac{k_1 + k_2}{2}\Delta x
\end{array}\right.} \]
という漸化式を得る事ができます。
</p>
</section>
<section style="font-size:70%">
<h2> 改良オイラー法の精度 </h2>
<p>
合成微分則より
\[ \frac{\mathrm{d}^2 f}{\mathrm{d} x^2} = \frac{\mathrm{d} g}{\mathrm{d} x} = \frac{\partial g}{\partial x} + \frac{\partial g}{\partial f}\frac{\mathrm{d} f}{\mathrm{d} x} \]
であったので,\(g\)のテイラー展開を行うと
\[ \begin{aligned}
g(x_i + \Delta x, f_i + k_1\Delta x) &= g(x_i, f_i) + \frac{\partial g}{\partial x}(x_i, f_i)\Delta x + \frac{\partial g}{\partial f}(x_i, f_i)k_1\Delta x + \mathcal{O}(\Delta x^2)\\
&= \frac{\mathrm{d} f}{\mathrm{d} x}(x_i, f_i) + \frac{\mathrm{d}^2 f}{\mathrm{d} x^2}(x_i, f_i)\Delta x + \mathcal{O}(\Delta x^2)
\end{aligned} \]
となります。従って,改良オイラー法の右辺は
\[ f_i + \frac{\mathrm{d} f}{\mathrm{d} x}(x_i, f_i)\Delta x + \frac{1}{2}\frac{\mathrm{d}^2 f}{\mathrm{d} x^2}(x_i, f_i)\Delta x^2 + \mathcal{O}(\Delta x^3) \]
となります。
</p>
<p>
これは左辺のテイラー展開の\(\mathcal{O}(\Delta x^2)\)の項まで一致しているので二次の精度となります。
</p>
</section>
<section>
<div class="block" style="border-color:pink;font-size:80%">
<h4 style="color:pink"> 改良オイラー法 </h4>
<p>微分方程式
\[ \frac{\mathrm{d} f}{\mathrm{d} x} = g(x,f) \]
の各点\(x_i\)における値を
\[ \left\{\begin{array}{l}
k_1 = g(x_i, f_i) \\
k_2 = g(x_i + \Delta x, f_i + k_1\Delta x) \\
f_{i+1} = f_i + \frac{k_1 + k_2}{2}\Delta x
\end{array}\right. \]
によって求める方法を<strong>改良オイラー法</strong>という。
</p>
<p>
精度は2次である。
</p>
</div>
</section>
<section style="font-size:70%">
<div class="block" style="border-color:lightgreen">
<h4 style="color:lightgreen">例</h4>
微分方程式の初期値問題
\[ \frac{\mathrm{d} x}{\mathrm{d}t} = x,\ x(0) = 1\]
を改良オイラー法で解きます。
</div>
<pre><code class="python" style="max-height:400px">>>> def g(t, x):
... return x
...
>>> x = 1
>>> dt = 0.1
>>> for i in range(10):
... t = i*dt
... print "{0} {1}".format(t, x)
... k1 = g(t, x)
... k2 = g(t + dt, x + k1*dt)
... x += (k1 + k2)*dt/2
...
0.0 1
0.1 1.105
0.2 1.221025
0.3 1.349232625
0.4 1.49090205062
0.5 1.64744676594
0.6 1.82042867636
0.7 2.01157368738
0.8 2.22278892456
0.9 2.45618176164
</code></pre>
</section>
<section style="font-size:80%">
<p>
先ほどのグラフに改良オイラー法の結果を重ねると下のようになります。
</p>
<div align="center"> <img width="500px" src="fig/ode-fig3.png"> </div>
</section>
<section style="font-size:80%">
<h2> ルンゲ=クッタ法 </h2>
<p>
改良オイラー法を一般化して
\[ \begin{aligned}
k_1 &= g(x_i, f_i) \\
k_2 &= g(x_i + \alpha_2 \Delta x, f_i + \beta_{21}k_1\Delta x) \\
k_3 &= g(x_i + \alpha_3 \Delta x, f_i + \beta_{31}k_1\Delta x + \beta_{32}k_2\Delta x) \\
&\vdots \\
k_s &= g(x_i + \alpha_s \Delta x, f_i + \beta_{s1}k_1\Delta x + \beta_{s2}k_2\Delta x+\cdots+\beta_{s,s-1}k_{s-1}\Delta x) \\
f_{i+1} &= f_i + (\gamma_1k_1 + \gamma_2k_2 + \cdots + \gamma_sk_s)\Delta x
\end{aligned} \]
という形の更新式を考える事が出来ます。係数\(\alpha,\beta,\gamma\)はテイラー展開の項が一致するように決定していきます。
</p>
<p>
これは<strong>ルンゲ=クッタ法</strong>と呼ばれます。
</p>
</section>
<section style="font-size:80%">
<p>
ルンゲ=クッタ法には陰的な物も含め様々な手法が存在しますが,ここでは代表的な公式を一つ紹介します。
</p>
<div class="block" style="border-color:pink;font-size:80%">
<h4 style="color:pink"> 古典的ルンゲクッタ法 </h4>
<p>微分方程式
\[ \frac{\mathrm{d} f}{\mathrm{d} x} = g(x,f) \]
の各点\(x_i\)における値を
\[ \left\{\begin{array}{l}
k_1 = g(x_i, f_i) \\
k_2 = g(x_i + \frac{1}{2}\Delta x, f_i + \frac{1}{2}k_1\Delta x) \\
k_3 = g(x_i + \frac{1}{2}\Delta x, f_i + \frac{1}{2}k_2\Delta x) \\
k_4 = g(x_i + \Delta x, f_i + k_3\Delta x) \\
f_{i+1} = f_i + \frac{k_1 + 2k_2+2k_3+k_4}{6}\Delta x
\end{array}\right. \]
によって求める方法を<strong>古典的ルンゲクッタ法</strong>という。
</p>
<p>
精度は4次である。
</p>
</div>
</section>
<section style="font-size:80%">
<div class="block" style="border-color:lightgreen">
<h4 style="color:lightgreen">例</h4>
微分方程式の初期値問題
\[ \frac{\mathrm{d}x}{\mathrm{d}t} = x(1-x),\quad x(0) = 0.01 \]
を古典的ルンゲクッタ法で解きます。
</div>
<pre><code class="python" style="max-height:400px">def g(t, x):
return x*(1-x)
x = 0.01 # x(0) = 0.01
dt = 0.5
for i in xrange(20):
t = dt*i
print "{0} {1}".format(t, x)
k1 = g(t, x)
k2 = g(t + dt/2, x + k1*dt/2)
k3 = g(t + dt/2, x + k2*dt/2)
k4 = g(t + dt, x + k3*dt)
x += (k1 + 2*k2 + 2*k3 + k4)*dt/6
</code></pre>
</section>
<section>
<p>
厳密解は\( x = \frac{1}{1+99\exp(-t)} \)となります。
</p>
<div align="center"> <img width="500px" src="fig/logistic.png"> </div>
</section>
<section>
<h2> その他の方法 </h2>
<p>
これまで述べた以外にも様々な発想に基づく解法が存在しますが、時間の関係上割愛します。是非調べてみて下さい。
</p>
<ul>
<li> <strong>多段階法</strong> <br>
\(x=x_{i+1}\)における値の計算に,\(x=x_i\)時点の情報だけでなくさらに過去の時点の情報も利用する手法。
</li>
<li> <strong>予測子・修正子法</strong> <br>
一旦陽解法で解を一通り計算し,その結果を高精度な解法や陰解法の入力に用いて解を改善する手法。
</li>
<li> 加速手法との組み合わせ </li>
<li> など </li>
</ul>
</section>
<section>
<h2> 連立一階常微分方程式 </h2>
<p>
\[ \color{yellow}{\begin{aligned}
\frac{\mathrm{d} u}{\mathrm{d} x} &= f(x, u, v) \\
\frac{\mathrm{d} v}{\mathrm{d} x} &= g(x, u, v) \\
\end{aligned}} \]
の様な連立一階常微分方程式も同様に解くことが出来ます。例えば,陽的オイラー法であれば
\[ \begin{aligned}
u_{i+1} = u_i + f(x_i, u_i, v_i)\Delta x \\
v_{i+1} = v_i + g(x_i, u_i, v_i)\Delta x \\
\end{aligned} \]
の様にそれぞれ更新すれば良いです。
</p>
</section>
<section>
<h2> 高階常微分方程式 </h2>
<p>
二階以上の導関数を含む常微分方程式
\[ \color{yellow}{f(x, \frac{\mathrm{d}f}{\mathrm{d}x}, \frac{\mathrm{d}^2f}{\mathrm{d}x^2},\cdots, \frac{\mathrm{d}^mf}{\mathrm{d}x^m})=0} \]
は,
\[
u_1 = \frac{\mathrm{d} f}{\mathrm{d} x},
u_2 = \frac{\mathrm{d} u_1}{\mathrm{d} x},
\cdots
u_m = \frac{\mathrm{d} u_{m-1}}{\mathrm{d} x},
\]
と各導関数に変数を割り当てれば,連立一階微分方程式とする事が出来ます。
</p>
</section>
<section style="font-size:80%" graphics="spring">
<div class="block" style="border-color:lightgreen">
<h4 style="color:lightgreen">例</h4>
第1回にやったばねの運動方程式
\[ m\frac{\mathrm{d}^2 x}{\mathrm{d} t^2} = -kx \]
を解いてみましょう。
</div>
<svg>
</section>
<section style="font-size:80%">
<p>
\[ v = \frac{\mathrm{d} x}{\mathrm{d} t} \]
とおけばこの方程式は
\[ \left\{\begin{array}{l}
\frac{\mathrm{d} x}{\mathrm{d} t} = v \\
\frac{\mathrm{d} v}{\mathrm{d} t} = -\frac{k}{m}x
\end{array}\right. \]
という連立微分方程式となります。
</p>
<p>
ベクトルを利用すれば
\[ \frac{\mathrm{d}}{\mathrm{d} t}
\left(
\begin{array}{c}
x \\ v
\end{array}\right)
= \left(
\begin{array}{c}
v \\ -\frac{k}{m}x
\end{array}\right)
\]
と記述出来ます。
</p>
</section>
<section style="font-size:70%">
<p> ルンゲクッタ法でのコーディング例です。関数gで毎回ベクトルを生成するのは非効率ですが,読みやすさを優先します。</p>
<pre><code class="python" style="max-height:400px">
import numpy as np
k = 1.0
m = 1.0
def g(t, x):
return np.array([x[1], -k/m*x[0]])
x = np.array([1.0, 0.0]) # 初期条件 x(0) = 1, v(0) = 0
dt = 0.1
for i in range(100):
t = i*dt
print "{0} {1}".format(t, x[0])
k1 = g(t, x)
k2 = g(t + dt/2, x + k1*dt/2)
k3 = g(t + dt/2, x + k2*dt/2)
k4 = g(t + dt, x + k3*dt)
x += (k1 + 2*k2 + 2*k3 + k4)*dt/6
</code></pre>
</section>
<section>
<p>
計算結果です。\(\Delta t\)を変えたり,初期条件を変えたり,他の解法を試したりいろいろ遊んでみて下さい。
</p>
<div align="center"> <img width="500px" src="fig/spring-rk4.png"> </div>
</section>
<section>
<p>
【練習問題】<br>
先ほどの方程式に抵抗と外部からの強制振動を加えた以下の方程式を解いて下さい。
\[ m\frac{\mathrm{d}^2 x}{\mathrm{d} t^2} = -kx -\mu\frac{\mathrm{d}x}{\mathrm{d} t}+A\sin(Bt)\]
パラメータは適当に\(m = 1, k = 1, \mu = 0.2, A=0.5\)などとして,
\(B\)の値をいろいろ変えて実験してみて下さい。共振現象が観測出来るはずです。
</p>
</section>
<section>
<h2> 偏微分方程式 </h2>
<p>
\[ \frac{\partial u}{\partial t}=k\frac{\partial^2 u}{\partial x^2} \]
といった偏微分方程式も同様に差分化を行って解くのですが,解が\(\pm\infty\)に発散してしまわない為の条件の解析などが難しくなります。(フーリエ展開・超関数などの知識が必要)
</p>
<p>
今回は是非知っておいて欲しい偏微分方程式を1つだけ紹介します。
</p>
</section>
<section graphics="randomwalk">
<p>
確率\(\frac{1}{2}\)で数直線上を左右に移動する点を考えましょう。これを<strong>一次元ランダム・ウォーク</strong>と言います。
</p>
<svg></svg>
</section>
<section>
<p>
時刻\(t\)に数直線上の\(x\)の位置に点がある確率を
\[ p(x, t) \]
とします。点は時刻\(\Delta t\)毎に距離\(\Delta x\)移動するとすれば
\[ \color{yellow}{p(x, t + \Delta t) = \frac{1}{2}p(x-\Delta x, t) + \frac{1}{2}p(x+\Delta x, t)} \]
という漸化式が成り立ちます。
</p>
<div align="center"> <img width="500px" src="fig/randomwalk.png"> </div>
</section>
<section style="font-size:80%">
<p>
ここで
\[ \begin{aligned}
p(x,t+\Delta t) &= p(x,t) + \frac{\partial p}{\partial t}\Delta t + \mathcal{O}(\Delta t^2) \\
p(x+\Delta x, t) &= p(x,t) + \frac{\partial p}{\partial x}\Delta x + \frac{1}{2}\frac{\partial^2 p}{\partial x^2}\Delta x^2 + \mathcal{O}(\Delta x^3) \\
p(x-\Delta x, t) &= p(x,t) - \frac{\partial p}{\partial x}\Delta x + \frac{1}{2}\frac{\partial^2 p}{\partial x^2}\Delta x^2 + \mathcal{O}(\Delta x^3) \\
\end{aligned} \]
と展開出来るので,先ほどの漸化式は
\[ \frac{\partial p}{\partial t} = \frac{1}{2}\frac{\Delta x^2}{\Delta t}\frac{\partial^2 p}{\partial x^2}+\mathcal{O}(\Delta t) + \mathcal{O}(\frac{\Delta x^3}{\Delta t}) \]
となります。
</p>
<p>
ここで\(\frac{\Delta x^2}{\Delta t}\)が定数となるようにしながら\(\Delta x,\Delta t\rightarrow 0\)の極限をとると
\[ \color{yellow}{\frac{\partial p}{\partial t} = k\frac{\partial^2 p}{\partial x^2}} \]
という偏微分方程式方程式が得られます。
</p>
</section>
<section>
<p>
導出の流れから判るように,これは熱や粒子などが拡散していく様子を記述します。
</p>
<div class="block" style="border-color:pink;font-size:80%">
<h4 style="color:pink"> 拡散方程式 </h4>
\[ \frac{\partial u}{\partial t} = D\frac{\partial^2 u}{\partial x^2} \qquad(\text{$D$は定数})\]
</div>
</section>
<section style="font-size:80%">
<p>
これを解くには例えば前進差分と二階差分を使って
\[ \begin{aligned}
\frac{\partial u}{\partial t} &\approx \frac{u(x_i, t_{j+1})-u(x_i, t_j)}{\Delta t} \\
\frac{\partial^2 u}{\partial x^2} &\approx \frac{u(x_{i+1},t_j)-2u(x_i,t_j)+u(x_{i-1},t_j)}{\Delta x^2}
\end{aligned} \]
と近似して
\[ u(x_i, t_{j+1}) = u(x_i, t_j) + D\frac{\Delta t}{\Delta x^2}(u(x_{i+1},t_j)-2u(x_i,t_j)+u(x_{i-1},t_j)) \]
とします。
</p>
<p>
詳細は割愛しますが,解が発散しない為には\(D\frac{\Delta t}{\Delta x^2}\leq \frac{1}{2}\)という条件が必要です。
</p>
</section>
<section style="font-size:70%">
<div class="block" style="border-color:lightgreen">
<h4 style="color:lightgreen">例</h4>
拡散方程式を\(0 \leq x \leq 2\)の範囲で解いてみます。境界条件は
\[ \frac{\partial u}{\partial x}(0, t) = \frac{\partial u}{\partial x}(2, t) = 0 \]
,初期条件は\(u(0, 0) = 1,\quad u(x, 0) = 0 (x \neq 0) \)とします。
</div>
<pre><code class="python" style="max-height:400px">
import numpy as np
D = 1
dt = 0.002
nx = 21
dx = 2.0/(nx-1)
print D*dt/(dx*dx)
u = np.zeros(nx)
u[nx/2] = 1.0
new_u = np.zeros(nx)
for j in range(50):
if j%5 == 0:
for i in range(nx):
print "{0} {1}".format(i*dx, u[i])
print "\n\n"
new_u[0] = u[1] # 境界条件。端は水平
new_u[nx-1] = u[nx-2] # 境界条件
for i in range(1, nx-1):
new_u[i] = u[i] + D*dt/(dx*dx)*(u[i+1]-2*u[i]+u[i-1])
new_u, u = u, new_u
</code></pre>
</section>
<section>
<div align="center"> <img width="600px" src="fig/diffusion.png"> </div>
</section>
<section style="font-size:90%">
<h2> 正規分布 </h2>
<p>
ところで拡散方程式\( \frac{\partial u}{\partial t} = D\frac{\partial^2 u}{\partial x^2} \)に
\[ u(x,t) = \frac{1}{2\sqrt{\pi D t}}\exp\left(-\frac{x^2}{4Dt}\right) \]
を代入すると成り立つ事が確認出来ると思います。従って,これは拡散方程式の解の一つです。
</p>
<p>
これは<strong>正規分布</strong>という統計学において重要な分布と同じ式の形をしています。完全にランダムな変位が累積すると結果として正規分布が現れるという事が解ります。何故正規分布が統計学において重要となるのか、理解出来るのではないでしょうか。
</p>
</section>
<section>
<h2 class="chapter-title"> 積分法 </h2>
</section>
<section>
<h2> 積分とは </h2>
<p>
微分法は関数の瞬間的な変化の様子を記述し調べる事が出来る学問でした。
積分法は微分法と表裏一体の関係にあり<strong>変化の集積の様子</strong>が主な関心の対象となります。
</p>
</section>
<section style="font-size:90%">
<h2> 定積分 </h2>
<div align="center"> <img width="500px" src="fig/riemann-integral.png"> </div>
<p>
区間\([a,b]\)を\(n\)個の区間に分割し,\(t_i\)を区間\([x_{i-1},x_i]\)内の適当な点とします。この時
\[ \sum_{i=1}^nf(t_i)(x_i-x_{i-1}) \]
をこの分割に関する<strong>リーマン和</strong>と言います。
\(f(t_i)(x_i-x_{i-1})\)は上図の様に長方形の符号付き面積です。
</p>
</section>
<section style="font-size:90%">
<p>
ここで\(n\rightarrow\infty, \lim_{n\rightarrow\infty}\max\{x_i-x_{i-1}\} = 0\)とした時の極限が収束するならば,これを\(f(x)\)の区間\([a,b]\)での<strong>定積分</strong>と言い
\[ \color{yellow}{\int_a^bf(x)\mathrm{d} x = \lim_{n\rightarrow\infty}\sum_{i=1}^nf(t_i)(x_i-x_{i-1})} \]
と表します。
</p>
<p>
但し,\(a > b\)の時は
\[\int_a^bf(x)\mathrm{d} x = -\int_b^af(x)\mathrm{d} x\]
と定義します。
</p>
<img align="right" width="400px" src="fig/riemann-integral2.png">
<p>
直感的には右図のような曲線\(y=f(x)\)と\(x\)軸によって囲まれた部分の符号付き面積の事です。
</p>
</section>
<section>
<h2> 定積分の公式 </h2>
<p>
定義より以下の公式を示す事が出来ます。
</p>
<div class="block" style="border-color:pink;font-size:90%">
\[ \begin{aligned}
&\int_a^b\{f(x)+g(x)\}\mathrm{d} x = \int_a^bf(x)\mathrm{d}x + \int_a^bg(x)\mathrm{d} x \\
&\int_a^bkf(x)\mathrm{d} x = k\int_a^bf(x)\mathrm{d} x\quad\text{($k$は定数)}\\
&\int_a^bf(x)\mathrm{d} x = \int_a^cf(x)\mathrm{d} x + \int_c^bf(x)\mathrm{d} x
\end{aligned} \]
</div>
</section>
<section style="font-size:90%">
<h2> 定積分についての平均値の定理 </h2>
<div class="block" style="border-color:pink;font-size:90%">
<p>
\(f(x)\)が区間\([a,b]\)で連続であるとき,ある\(c \in (a,b)\)が存在して
\[ \int_a^b f(x)\mathrm{d} x = f(c)(b-a) \]
が成り立つ。
</p>
</div>
<p style="font-size:70%">
【証明】<br>
区間\([a,b]\)での\(f(x)\)の最大値・最小値を\(M,m\)とすると,定積分の定義より
\[ m(b-a) \leq \int_a^b f(x)\mathrm{d} x \leq M(b-a) \]
となる。ところで\(f(x)(b-a)\)も連続で最大値が\(M(b-a)\),最小値が\(m(b-a)\)であるから,
\[ \int_a^b f(x)\mathrm{d} x = f(c)(b-a) \]
を満たす\(c\)が\(c \in (a,b)\)に存在する(中間値の定理)。
</p>
</section>
<section style="font-size:90%">
<h2> 不定積分 </h2>
<p>
\[ \frac{\mathrm{d}}{\mathrm{d} x}F(x) = f(x) \]
を満たす\(F(x)\)を\(f(x)\)の<strong>原始関数</strong>と言います。
</p>
<p>
例えば
\[ \frac{\mathrm{d}}{\mathrm{d} x}x^2 = 2x \]
なので\(x^2\)は\(2x\)の原始関数です。
</p>
<p class="fragment">
ここで\(f(x)\)の原始関数全体を\(f(x)\)の<strong>不定積分</strong>と呼び
\[ \color{yellow}{\int f(x)\mathrm{d} x} \]
と表します。
</p>
</section>
<section style="font-size:90%">
<p>
\(f(x)\)の原始関数\(F_1,F_2\)について
\[ \frac{\mathrm{d}}{\mathrm{d}x}\{F_2(x)-F_1(x)\} = f(x)-f(x) = 0 \]
となるので\(F_2(x)-F_1(x)=C\)と表せます。つまり原始関数は定数項が異なるもののみとなります。
</p>
<div class="block" style="border-color:pink;font-size:90%">
<h4 style="color:pink"> 不定積分 </h4>
\(f(x)\)の原始関数の一つを\(F(x)\)とすると
\[ \int f(x)\mathrm{d} x = F(x) + C \quad\text{($C$は任意定数)} \]
と表される。
</div>
</section>
<section style="font-size:80%">
<div class="block" style="border-color:pink;font-size:90%">
<h4 style="color:pink"> 基本関数の不定積分 </h4>
\[ \begin{aligned}
&\int x^\alpha\mathrm{d} x = \frac{1}{\alpha + 1} x^{\alpha + 1}+C\qquad (\alpha \neq -1) \\
&\int \frac{1}{x}\mathrm{d} x = \ln|x|+C \\
&\int \sin x \mathrm{d} x = -\cos x+C,\quad \int \cos x \mathrm{d} x = \sin x + C \\
&\int \tan x\mathrm{d} x = -\ln|\cos x| + C \\
&\int \exp x \mathrm{d} x = \exp x + C,\quad \int \ln x \mathrm{d} x = x(\ln x - 1) + C
\end{aligned} \]
</div>
<p>
これらが成立する事は微分法を利用して簡単に確認する事が出来ます。
</p>
</section>
<section>
<div class="block" style="border-color:pink;font-size:90%">
<h4 style="color:pink"> 置換積分の公式 </h4>
<p> \(x = g(t)\)のとき
\[ \int f(x)\mathrm{d} x = \int f(g(t))g'(t)\mathrm{d} t \]
</p>
</div>
<div class="block" style="border-color:pink;font-size:90%">
<h4 style="color:pink"> 部分積分の公式 </h4>
\[ \int f(x)g'(x)\mathrm{d} x = f(x)g(x) - \int f'(x)g(x)\mathrm{d} x \]
</div>
<p>
これらも非常によく使われる公式です。両辺を微分する事によってその成立を容易に確認する事が出来ます。
</p>
</section>
<section>
<h2> 微積分学の基本定理 </h2>
<p>
以上で見たように元々<strong>定積分</strong>と<strong>微分・不定積分</strong>は全く別に定義されています。
これらを結びつけるのが<strong>微積分学の基本定理</strong>です。
</p>
<div class="block" style="border-color:pink;font-size:90%">
<h4 style="color:pink"> 微積分学の基本定理 </h4>
<p>
\[ \int_a^xf(t)\mathrm{d} t \]
は\(f(x)\)の原始関数である。すなわち
\[ \frac{\mathrm{d}}{\mathrm{d} x}\int_a^xf(t)\mathrm{d} t = f(x) \]
また,\(F(x)\)が\(f(x)\)の原始関数であるとき
\[ F(b)-F(a) = \int_a^b f(x)\mathrm{d} x \]
</p>
</div>
</section>
<section style="font-size:80%">
<p>
【証明】<br>
\[ G(x) = \int_a^x f(t)\mathrm{d} t \]
とおくと
\[ \begin{aligned}
\frac{G(x+\Delta x)-G(x)}{\Delta x} &= \frac{1}{\Delta x}\left\{ \int_a^{x+\Delta x}f(t)\mathrm{d} t - \int_a^x f(t)\mathrm{d} t\right\} \\
&= \frac{1}{\Delta x}\int_x^{x+\Delta x}f(t)\mathrm{d} t \\
&= f(c) \qquad \text{($c$は$x$と$x+\Delta x$の間に存在)}
\end{aligned} \]
となるので\(\Delta x \rightarrow 0\)の極限を両辺でとって
\[ \frac{\mathrm{d}}{\mathrm{d}x}G(x) = \frac{\mathrm{d}}{\mathrm{d} x}\int_a^xf(t)\mathrm{d} t = f(x) \]
となる。
</p>
</section>
<section style="font-size:80%">
<p>
【証明続き】<br>
\[ G(x) = \int_\alpha^x f(t)\mathrm{d} t \]
は\(f(x)\)の原始関数であったので,任意の原始関数\(F(x)\)は
\[ F(x) = \int_\alpha^x f(t)\mathrm{d} t + C \]
と表せる。従って
\[ F(b)-F(a) = \int_\alpha^bf(t)\mathrm{d} t - \int_\alpha^a(t)\mathrm{d} t = \int_a^b f(x)\mathrm{d} x \]<span style="float:right">□</span> </p>
</p>
</section>
<section>
<div class="block" style="border-color:lightgreen">
<h4 style="color:lightgreen">例</h4>
\[ \int_0^\frac{\pi}{4}\tan^2 x\mathrm{d} x \]
</div>
<p>
\(\tan^2 x = \frac{1}{\cos^2 x}-1 \)で,\((x)' = 1,\ (\tan x)'=\frac{1}{\cos^2 x}\)だから,\(C\)を任意定数として
\[ \int \tan^2 x\mathrm{d} x = \tan x - x + C \]
従って
\[ \int_0^\frac{\pi}{4}\tan^2 x\mathrm{d} x = (\tan\frac{\pi}{4}-\frac{\pi}{4})-(\tan 0-0) = \color{yellow}{1 - \frac{\pi}{4}} \]
</p>
</section>
<section style="font-size:90%">
<h2> 広義積分 </h2>
<p>
統計の計算では
\[ \int_{-\infty}^\infty f(x)\mathrm{d} x \]
の様な積分が頻繁に登場します。このような積分は<strong>広義積分</strong>と呼ばれるものの一つです。
</p>
<p>
この場合の定義は単純で
\[ \lim_{t\rightarrow\infty}\int_a^t f(x)\mathrm{d} x \]
が存在するときにこれを
\[ \int_a^\infty f(x)\mathrm{d} x \]
と表します。下端についても同様です。
</p>
</section>