2024年度「数値解析」のページ

お知らせ

資料

• 浮動小数点数
  • 「浮動小数点数」
• 補間法
  • 「補間法」
  • 「Vandermonde多項式」
• 非線形一変数方程式の数値解法
  • 「二分法」
  • 「ニュートン法」
  • 「ニュートン法 補足：収束の速さ」
  • 三角形の辺の長さからマルファッティの円の半径を計算するプログラム（短い解説つき）
• 連立線型方程式の数値解法
  • 「LU分解法」
  • LU分解法をCの関数として実装したプログラム
• 連立非線型方程式の数値解法
  • 「多変数のニュートン法」
• 収束の加速
  • 「エイトケン加速」
• 数値積分
  • 「台形公式とシンプソン公式」
• 常微分方程式の数値解法
  • 「常微分方程式の初期値問題」
• 「Cの入門書のようなもの（仮）」（題名募集中）
  • 2023年4月10日暫定版

練習問題

問題の解説などは受講生がCプログラムを書けることを想定しているが、課題のプログラムは他言語で書いても良い。ただし、浮動小数点数の計算が可能な言語に限る。

第1回課題

1. \(\sum_{n=1}^{2^{28}}(1/n^2)\) を計算するプログラムを、 以下のそれぞれのアルゴリズムを実装する形で書き、実行し、誤差の考察をせよ。 それぞれ、floatとdoubleで実装して、違いについても実験・考察せよ。

   1. 初項から正順に足していくアルゴリズム
   2. 絶対値の小さなものから順に足していくアルゴリズム
   3. 再帰的に部分和を約半分ずつ分割して足していくアルゴリズム
   4. カハンのアルゴリズム

   ヒント1

   \[ \sum_{n=1}^{\infty}\frac{1}{n^2} = \frac{\pi^2}{6} \]

   ヒント2

   \(2^{28}=16^7\)であるので、\(2^{28}\)はCでは 0x10000000 と書ける。

   \(2^{28}\)がintの範囲におさまるかどうかは処理系依存である。移植性を重視するなら、ループカウンタにint型の変数は使わずにlong型を使うが安全である。

   int32_t型を使うと、さらに安全確実になる。int32_t型をサポートする処理系ならば、そのビット幅は確実に32ビットである。なお、これを使うには、
   #include <stdint.h>
   が必要である。

   ヒント3

   一般に、\(n\)を\(1\)から\(N\)まで動かして\(1/n^2\)を使って何かするには、例えば、次のようなコードを書くとうまくいく。

   for (long n = 1; n <= N; n ++) {
           double  t = n;
           double  x = 1 / (t * t);
           /* ここで x を使って何かする */
   }

   • ループカウンタは整数型の変数にしたい。浮動小数点数をループカウンタに使うと誤差の蓄積が怖い。
   • \(1/n^2\)を計算するには浮動小数点数型が必要。整数型の割り算は小数点以下切り捨てなので、整数型のまま割り算してもうまくいかない。

   そのため、整数型のループカウンタの値をいったん浮動小数点数に変換して使用する形にすれば良い。Cなどのプログラミング言語では、整数型の変数から浮動小数点数型の変数に代入を行うと、型変換を行うコードを自動的に追加してくれる。

第2,3回課題

1.   \[ \begin{aligned} f(-2) &= 1/4, \\ f(-1) &= 1/2, \\ f(0) &= 1, \\ f(1) &= 2, \\ f(2) &= 4 \end{aligned} \] から以下のそれぞれの補間法で\(f(-3)\), \(f(-1/2)\), \(f(3/2)\), \(f(10)\)の近似値を計算し、\(f(x)=2^x\) との誤差を評価せよ。

   1. ラグランジュ補間
   2. ニュートン補間

第4,5回課題

1. 浮動小数点数 \(a\) を入力し、\(x\) についての方程式 \(xe^x = a\) の実数解を出力するプログラムを書け。ただし、\(a\gt0\) と仮定してよい。（仮定しないプログラムを書いたらかっこいい）
   1. 二分法を用いよ。
   2. ニュートン法を用いよ。
   3. 簡易ニュートン法を用いよ。

   ヒント1

   \(a\gt0\) と仮定し、\(f(x)=xe^x-a\)　とおく。 このとき、以下が成り立つ。

   • \(f(0)\lt0\) かつ \(f(a)\gt0\) が成り立つ。
   • \(x\gt0\) のとき \(f'(x)\gt0\) かつ \(f''(x)\gt0\) が成り立つ。

   したがって、解は区間\((0,a)\)にちょうど一つあることと、\(f\)は\((0,+\infty)\)で単調増加かつ下に凸であることがわかる。

   ヒント2

   ヒント1の事実を利用すると、初期値を以下のように選ぶとうまくいくことがわかる。

   二分法の場合

   初期区間を\((0,a)\)にとる。

   ニュートン法の場合

   初期値を\(a\)にとる。

2. 浮動小数点数 \(a\) を入力し、\(x\) についての方程式 \(x^3-3ax^2-3x+a = 0\) の実数解をすべて出力するプログラムを書け。ただし、\(a\gt0\) と仮定してよい。（仮定しないプログラムを書いたらかっこいい）
   1. 二分法を用いよ。
   2. ニュートン法を用いよ。
   3. 簡易ニュートン法を用いよ。

   ヒント1

   \(f(x)=x^3-3ax^2-3x+a\)　とおくと \(a\gt0\) のとき以下が成り立つ。 \[ \begin{gathered} f(-2)\lt0, \\ f(0)\gt0, \\ f(1)\lt0, \\ f(a)\lt0, \\ f(3a+2)\gt0. \end{gathered} \]

第6,7回課題

1. LU分解法をCの関数として実装したプログラムを他の言語に移植せよ。

第8,9回課題

1. 三角形\(ABC\)の外心を\(O\)とし、外接円半径を\(r_O\)とする。外心は頂点から等距離にあることから、\[\tag{☂︎}\overline{OA}^2=\overline{OB}^2=\overline{OC}^2={r_O}^2 \] が成り立つ。 頂点と外心の座標を \(A=(x_A,y_A)\), \(B=(x_B,y_B)\), \(C=(x_C,y_C)\), \(O=(x_O,y_O)\) とおくと、(☂︎)は \[ \tag{☂︎☂︎} \left\{ \begin{aligned} (x_O-x_A)^2+(y_O-y_A)^2-{r_O}^2 &= 0 \\ (x_O-x_B)^2+(y_O-y_B)^2-{r_O}^2 &= 0 \\ (x_O-x_C)^2+(y_O-y_C)^2-{r_O}^2 &= 0 \end{aligned} \right. \] と書ける。 したがって、(☂︎☂︎)を\(x_O\)と\(y_O\)と\(r_O\)の方程式とみなして解くと、三角形の頂点の座標から外心の座標と外接円半径を計算できる。

   1. (☂︎☂︎)にニュートン法を適用して数値的に解くプログラムを書け。また、簡易ニュートン法を適用して数値的に解くプログラムも書け。
   2. (☂︎☂︎)を代数的に解いてその結果を実装したプログラムを書け。

   ヒント

   aは、 \[ \begin{gathered} x^{(n)}= \begin{bmatrix} x_O^{(n)} \\ y_O^{(n)} \\ r_O^{(n)} \end{bmatrix} \\ f(x_O, y_O, r_O)= \begin{bmatrix} (x_O-x_A)^2+(y_O-y_A)^2-{r_O}^2 \\ (x_O-x_B)^2+(y_O-y_B)^2-{r_O}^2 \\ (x_O-x_C)^2+(y_O-y_C)^2-{r_O}^2 \end{bmatrix} \end{gathered} \] とおいて、多変数のニュートン法を適用すればよい。

第10,11回課題

1. 以下のそれぞれの数列について、その極限値を素朴に計算するプログラムとエイトケン加速を適用するプログラムとエイトケン加速を2回適用するプログラムを書き、収束の速さを比較せよ。

   1. \[ \tag{✈︎} \begin{aligned} x_0 &= 1 \\ x_{n+1} &= 1 + \frac{1}{x_n} \end{aligned} \]
   2. \[ \tag{⚓︎} \sum_{k=1}^{2^n} \frac{1}{k^2} \]

第12,13回課題

1. \[ \int_0^1 \sqrt{1-x^2}\,dx \] を数値積分で計算するプログラムを書け。 区間分割数をいろいろと変えて、精度の変化を観測せよ。
   1. 台形公式を使って計算せよ。
   2. シンプソン公式を使って計算せよ。
2. \[ \int_0^1 \frac{dx}{1+x} \] を数値積分で計算するプログラムを書け。 区間分割数をいろいろと変えて、精度の変化を観測せよ。
   1. 台形公式を使って計算せよ。
   2. シンプソン公式を使って計算せよ。
3. \[ \int_0^1 \frac{dx}{1+x^2} \] を数値積分で計算するプログラムを書け。 区間分割数をいろいろと変えて、精度の変化を観測せよ。
   1. 台形公式を使って計算せよ。
   2. シンプソン公式を使って計算せよ。

第14,15回課題

1. 以下の常微分方程式を以下のそれぞれの方法で\(0\le t\le 2\)の範囲で数値的に解き、精度を比較せよ。 \[ \frac{dx}{dt}=2xt,\qquad x(0)=1. \]

   1. オイラー法
      1. 刻み幅\(0.1\)
      2. 刻み幅\(0.01\)
      3. 刻み幅\(10^{-3}\)
      4. 刻み幅\(10^{-4}\)
      5. 刻み幅\(10^{-5}\)
      6. 刻み幅\(10^{-6}\)
   2. 修正オイラー法
      1. 刻み幅\(0.1\)
      2. 刻み幅\(0.01\)
      3. 刻み幅\(10^{-3}\)
      4. 刻み幅\(10^{-4}\)
      5. 刻み幅\(10^{-5}\)
      6. 刻み幅\(10^{-6}\)
   3. ルンゲ・クッタ法
      1. 刻み幅\(0.1\)
      2. 刻み幅\(0.01\)
      3. 刻み幅\(10^{-3}\)
      4. 刻み幅\(10^{-4}\)
      5. 刻み幅\(10^{-5}\)
      6. 刻み幅\(10^{-6}\)

この常微分方程式は解析的に解ける。

講義録画

• （準備中）