2023 年度「計算数学ａ」 2023-04-28

今回のプリントのうち、乱数以外は、重要ではあるが今後この授業で使うことはないので、時間の都合などのある人は深入りしなくてもよい。途中までをこなし、あとは後日、でもよい。

§4.1 数値解析・数値計算とは

方程式の近似解を求めたり、定積分の近似値を求めたり、微分方程式の近似解を求めたりすることをいう。

解の公式のない方程式の場合、定積分が初等関数の組み合わせで書けない場合、解が初等関数の組み合わせで書けない微分方程式の場合を想定するが、ここでは実際は答えがわかるものを中心にプログラムを書く。

また、非常に素朴なアルゴリズムのみを説明する。

§4.2 ニュートン法

方程式 f(x) = 0 の近似解を求める方法の一つに、ニュートン法がある。

a を解に近い値としよう。点 (a, f(a)) における y = f(x) のグラフの接線が、 x 軸と交わる点の x 座標を a' とすると、 a' は a よりも精度のよい近似値になっている場合がある。絵を書いてみて確かめよ。これをくり返して解の近似値を求めるのがニュートン法である。

a' を a と f(a) と f'(a) の式で書く方法は、諸君には簡単であろうから省略する。各自考えよ。

x² - 2 = 0 を解くことで、√2 の近似値を求めるプログラム例。

#include <stdio.h>

int main() {
  int n;
  double x = 2;  /* 最初にとった、解に近い値 */

  for (n = 0; n < 20; n++) {
    printf("%.20f\n", x);
    x = x - (x*x-2)/(2*x);
  }
}

§4.3 練習問題――初期値の選び方と収束先

x³ - x = 0 の解をニュートン法で求めるプログラムを書け。

解は三つある。そのうちのどれに収束するかは、最初にとった値に依存する。いろいろ試せ。

±1/√5 から始めると、計算の上では ±1/√5 を交互にくり返す。試してみよ。

§4.4 定積分の近似値

次の公式は高等学校の数学 III で習った区分求積法である。

n-1

(

)

(

)

∫

lim

---

∑

---

lim

---

∑

---

f(x)dx

n→∞

k=0

n→∞

k=1

このことから、n を十分大きくとれば積分の近似値が ∑f(k/n)/n として求められそうだとわかる。

f(x) = x² の 0 から 1 までの定積分の近似値を求めるプログラム。

#include <stdio.h>

#define MESH (1.0/65536)        /* 1/n のこと */

int main() {
  double x, sum;

  sum = 0;
  for (x = 0; x < 1; x = x + MESH) {
    sum = sum + x*x * MESH;
  }
  printf("%.10f\n", sum);

  sum = 0;
  for (x = MESH; x <= 1; x = x + MESH) {
    sum = sum + x*x * MESH;
  }
  printf("%.10f\n", sum);
}

#define は、プログラムの中の MESH を (1.0/65536) に置き換えてからコンパイルせよ、の意味である。 MESH の値を変えて実験する際、この (1.0/65536) だけを書き換えればよいので楽である。この値は「（2 のベキ乗）分の 1」がよさそうである。

1.0/65536 をかっこでくくるのは、一種の定跡・定石である。かっこのある場合とない場合とで、（ここには現れないが） 1/MESH がどう置き換わるか、興味のある者は考えてみよ。

答えは二つ出力される。一つめは k = 0 から n-1 までの公式を、二つめは k = 1 から n までの公式を用いた。 x² は積分範囲 [0, 1] で単調増加関数だから、一つめの値は実際よりも小さく、二つめの値は実際よりも大きい。

MESH の値を変えて実験してみよ。「（2 のベキ乗）分の 1」以外にしたらどうなるか。

§4.5 練習問題――標準正規分布

高等学校の数学 B の教科書の、確率・統計の単元に載っているように、標準正規分布は f(x) = exp(-x²/2)/√(2π) で与えられる。これを、0 から 1, 2, 3 まで積分するプログラムを書け。（この関数の不定積分は初等関数の組み合わせで表せないことが知られている。）

#include <math.h> を忘れずに。円周率は #define PI 3.141592653589793238462643383279 をコピーして利用せよ。（double 型で計算するにはもっと少ない桁数でじゅうぶんであるが。）

一つのプログラムで 1, 2, 3 まで積分した値を出力させるものも書けるが、むずかしければ、三つのプログラムにすればよい。（あるいは、一つだけ出力させてみればよい。）

正規分布表はネット上でみつかる。それに近い値になるだろうか？

§4.6 微分方程式

「微分しても元のままの関数を求めよ」、すなわち「y' = y を満たす x の関数 y = y(x) を求めよ」のような問題を解くことを「微分方程式を解く」という。y' = y は微分方程式である。

上の問題の答えは y = C exp(x)（C は定数）である。これが答えであることはすぐわかる。逆に、このような y = y(x) があったとして、f(x) = exp(-x) y を考える。積の微分法により f'(x) = -exp(-x) y + exp(-x) y' となるが、この右辺は y' = y により 0 となる。よって f(x) は定数関数、ゆえに y = C exp(x) の形に書ける。

次が、y' = y を仮定「x = 0 のとき y = 1」のもとに数値的に解くプログラムである。このことを、「初期条件 y(0) = 1 のもとで」ということもある。 x の値、数値的に求めた y の値、理論的に求めた y の値（すなわち exp(x)) を並べて出力する。

#include <stdio.h>
#include <math.h>

#define MESH (1.0/4096)

int main() {
  double y, x;

  y = 1;
  for (x = 0; x < 1; x = x + MESH) {
    printf("%f %f %f\n", x, y, exp(x));
    y = y + y * MESH;
  }
}

x が 0 の瞬間の y の増加率は 1 である。次の瞬間――厳密な言い方ではないが――には、y は 1 より大きい。よって、y の増加率も 1 より大きい。しかし、このプログラムでは、x が 0 から MESH まで増えるあいだ、 y の増加率はずっと 1 であると仮定している。次は x が MESH から 2*MESH までのあいだであるが、 y の増加率は x が MESH のときの y の値、とする。以下同様。

§4.7 練習問題――単振動と振り子

変位 x を時刻 t の関数とする。 x の t による 2 階の微分が -x に等しい関数、すなわち x" = -x を満たす t の関数 x = x(t) は単振動 x = C₁sin(t) + C₂cos(t) （C₁, C₂ は定数）であることが知られている。ばねを引きのばして手を離した場合の運動はこれである。

微分方程式 x" = -x を、初期条件 x(0) = 1, x'(0) = 0 のもとで数値的に解くプログラムを書け。 x' の値を入れておく変数 v が必要である。（v は velocity（速度）の v である。）

MESH だけ時間がたつと、 v は v +（v の微分）* MESH に、 x は x +（x の微分）* MESH に変わる。理論的な解は x = cos(t), v = -sin(t) である。それと並べて出力されるようにせよ。

次は時間に余裕のある人だけが取り組めばよいだろう。振り子の運動は、重力加速度を 1 とすると x" = -sin(x) に従う。これの解は初等関数の合成で書けないことが知られている。これを解くプログラムを書け。振幅が大きすぎないよう、x(0) = 0.2 ぐらいで試せ。（x = 0.2 * cos(t), v = -0.2 * sin(t) に近いが、一致はしないはずである。）

§4.8 乱数

(4.8.1) 「乱数」が数学でどう定義されるかは、ここでは考えない。計算機でいう乱数とは、種（たね）と呼ばれる整数を元に数式を用いて発生させる、乱数のように見える数列のことである。

(4.8.2) 関数 rand() は、小カッコの中には何も書かずに呼び出す。呼び出されるたびに、 0 以上 RAND_MAX 以下の整数（int 型）に値をとる乱数の項を、一つずつ返す。次のプログラムでは、一行に一つずつ、全部で 10 個の整数が出力される。その並びが乱数、というわけである。

#include <stdio.h>
#include <stdlib.h>

int main() {
  int i;

  for (i = 0; i < 10; i++) {      /* 10 項からなる乱数列を発生 */
    printf("%d\n", rand());
  }
}

(4.8.3) 2 行目。関数 rand() を使うには stdlib.h を #include することが必要である。どの関数を使う際にどの .h ファイルが必要になるのかは、この授業ではそのつど教える。また、K&R2 の付録に書いてある。

(4.8.4) RAND_MAX はコンパイラごとに違う値を取り得る定数である。次のプログラムで調べられる。

#include <stdio.h>
#include <stdlib.h>

int main() {
  printf("%d\n", RAND_MAX);
}

この演習室の gcc では、その値は 32767 (= 2¹⁵ - 1) である。意外と小さいことに注意。

(4.8.5) (4.8.2) のプログラムは、毎回同じ乱数列を出力する。それでは使いづらいので、現在時刻を乱数の種にすることを考える。（#include <time.h> は time() 関数を使うために必要である。）

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

int main() {
  int i;

  {   /* この中カッコの中（乱数の種の設定）はわからなくてもよい */
    unsigned seed = (unsigned)time(NULL);   /* 現在時刻を取得して */
    printf("乱数の種は %u です.\n", seed);  /* その値を出力 */
    srand(seed);                            /* それを乱数の種に */
  }

  for (i = 0; i < 10; i++) {
    printf("%d\n", rand());
  }
}

ここのセンターのパソコンでは、ほかの多くのコンピュータと同じく、 time() は世界協定時での 1970 年 01 月 01 日 00 時 00 分 00 秒からの経過秒数を返す。ただし、秒未満は無視される。

また、1 秒に一回ぐらいのテンポで動かすと、乱数の種が 1 だけ変わっても、乱数の第 1 項はあまり変化しないことに気づくだろう。

(4.8.6) rand() % 100 とすると、 0 以上 100 未満の整数に値をとる乱数が得られる。厳密に言えば、RAND_MAX の値によっては、 0 の出現率と 99 の出現率はわずかながら異なるが、この授業ではそこまでは気にしないことにしよう。

(4.8.7) 練習：(4.8.5) のプログラムを、(4.8.6) を参考にして改造し、さいころのように、1 から 6 までの整数に値をとる乱数が発生されるようにせよ。また、実行するごとにその乱数列が異なることを確認せよ。

§4.9 モンテカルロ法

(4.9.1) 例で示す。

(4.9.2) 関数 rand() は 0 以上 RAND_MAX 以下の整数に値をとる乱数を返すのだった。その値を RAND_MAX + 1 で割ると、 0 以上 1 未満の小数に値をとる乱数が得られる。それを二度発生させ、その値を x, y とすると、 x² + y² が 1 以下となる確率は、「半径 1 の円の四分の一の面積」割る「一辺の長さが 1 の正方形の面積」で、 π/4 となる。このことを利用して、実験結果から π の概算値を求めることができる。

(4.9.3) 「x = rand() / (RAND_MAX + 1);」でなく「x = 1.0 * rand() / (RAND_MAX + 1);」となっている理由は (3.6.2) を参照。

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

#define N 10000000

int main() {
  int n, in;
  double x, y;

  srand((unsigned)time(NULL));  /* 乱数の種を現在時刻で初期化 */

  in = 0;                       /* 何個が半径 1 の四分の一円にはいったかを数えるカウンター */
  for (n = 0; n < N; n++) {
    x = 1.0 * rand() / (RAND_MAX + 1);
    y = 1.0 * rand() / (RAND_MAX + 1);
    if (x*x + y*y < 1) {
      in++;
    }
  }
  printf("%f\n", 4.0 * in / N); 
}

(4.9.4) 次のようにしても、ここのコンピュータでは同じこと。（(RAND_MAX + 1) * (RAND_MAX + 1) が int 型に収まるため。）

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

#define N 10000000

int main() {
  int n, in, x, y;

  srand((unsigned)time(NULL));  /* 乱数の種を現在時刻で初期化 */

  in = 0;                       /* 何個が半径 1 の四分の一円にはいったかを数えるカウンター */
  for (n = 0; n < N; n++) {
    x = rand(); y = rand();
    if (x*x + y*y < (RAND_MAX + 1) * (RAND_MAX + 1)) {
      in++;
    }
  }
  printf("%f\n", 4.0 * in / N); 
}

付：C 言語を選んだ理由

私が全貌を学んだことのある唯一のプログラム言語だから、が一番の理由です。そうでない言語を選ぶと、中学二年生が一年生に英語を教えていて「God Save the King は God が三人称単数なのに save に s がついていないから間違い」と言ってしまうのに似た現象が起きかねない、と考えています。

C 言語のうち、この授業で扱う範囲は、 1989 年に ANSI C という規格で仕様が明確に定められているので、動作がはっきり決まっています。（ただし、ANSI C は int 型で扱える整数の範囲を -32768 から 32767 まででもよいとしているので、すべてのプログラム例がすべてのコンピューター、すべてのコンパイラで動くわけではありません。）

いまはほかの言語が多く使われていますが、C 言語は決して現役を引退した言語ではありません。

金沢大学大学院自然科学研究科博士前期課程入学試験の、数物科学専攻（計算科学コース）の問題のうち、プログラムを書けという問題は、 C 言語または Fortran で、と指定されています。

岩瀬順一