2021 年度「計算数学ａ」 2021-10-15

今回とりあげることがらは、重要ではあるが、今後、この授業で使うことはないので、時間の都合などのある人は深入りしなくてもよい。途中までをこなし、あとは後日、でもよい。

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

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

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

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

§3.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);
  }
}

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

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

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

±1/√5 から始めると、計算の上では ±1/√5 を交互にくり返すはずなのだが、プログラムでは誤差のため収束してしまうようである。

§3.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」以外にしたらどうなるか。

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

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

#include <math.h> を忘れずに。円周率は #define PI 3.141592653589793238462643383279 をコピーして利用せよ。（double 型で計算するにはもっと少ない桁数でじゅうぶんであるが。また、字母 π の呼び名は pi であり、これを「パイ」と読むのは英語式である。）

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

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

§3.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 の値、とする。以下同様。

§3.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) に近いが、一致はしないはずである。）

付：コマンドライン引数の使い方

プロンプトの出ているところに「a 3 4.56」または「a.exe 3 4.56」のように打ち込んで動かし、プログラムの側でそれらの値を利用する方法がある。興味のある者は、以下の例で学べ。 x に 3 が、y に 4.56 がはいるだけのプログラムである。この「3 4.56」をコマンドライン引数と呼ぶ。

#include <stdio.h>
#include <stdlib.h> /* atoi(), atof() を使うのに必要 */

int main(int argc, char *argv[ ]) {     /* 小かっこの中は、常にこう書く。きょうは説明しない */
  int x;
  double y;

  if (argc != 3) {                      /* argc は、コマンドラインに打ち込んだものの個数 */
    printf("使い方が違います.\n");
  } else {
    x = atoi(argv[1]);                  /* int 型に収めるときは atoi() */
    y = atof(argv[2]);                  /* double 型に収めるときは atof() */
    printf("x = %d, y = %f.\n", x, y);  /* きちんと収まったか、確かめる（だけ）*/
  }
}

岩瀬順一