数値解析の基礎

目次

初回である今回は,C言語の復習と数値解析の基礎知識習得を兼ねて,簡単なプログラムを作成してもらいます.この授業で扱うプログラムは,次回以降,難度を増していくので,今の内にしっかりとC言語の復習をしておいてください.

計算機による数値解析

様々な装置の挙動や自然の振る舞いなどを予測する際,我々は理論的なモデル(方程式)を仮定し,そのモデルを数学的に解くことで振る舞いを予測します.簡単な現象であれば,それを表すモデルも簡単なものとなり,それらは解析的に(=式変形などをしていくことで理論的に正確な解を得ること)解くことができます.

ところがモデルが複雑になってくると,解析的な解を得ることは極めて難しくなります.そのような場合に,計算機の助けを借りて数値的に(=具体的な数値を当てはめながら解を予測すること)問題を解くことがあります.そのような問題解決手法を,数値解析とか数値シミュレーションと呼んでいます. 一般に,数値解析プログラムには,次のような条件が必要とされます.

いずれも重要な条件ですが,数値解析の目的を考えた場合,「誤差が少ないこと」は特に重要です.計算機の扱える数値には限界があるため,数値解析では必ず誤差が生じてしまいます.そのため,数値解析では常に誤差の存在を意識する必要があります.

今回の課題では,C言語の構造を簡単に復習しながら,数値解析で生ずる誤差を実際に体験してもらいたいと思います.

復習: C言語のプログラム

まずは,以下のサンプルプログラムを使って,C言語の基本的な構造を簡単に復習していきましょう(なお,次のプログラムリストでは,各行の左側に “1” というように行番号がついていますが,これは説明のためにつけたものであり,実際にプログラムを入力する際には不要です

1
2
3
4
5
6
7
8
9
10
11
12
/* ヘッダファイルのインクルード */
#include <stdio.h>

int main(int argc, char **argv) {
    // コマンドラインから変数を読みとる
    int a, b;
    scanf("%d %d", &a, &b);

    // 割り算
    const double c = (double)a / (double)b;
    printf("%d / %d = %f\n", a, b, c);
}

1行目はプログラムを読みやすくするためのコメントです.C言語には現在2種類のコメントがあります.C言語の古い形式であるC90では /**/ で囲まれた部分をコメントとしていましたが,98年に行われた言語仕様の改定を取り入れたC99という形式では // と書かれると,その行内でその場所以降がコメントとして扱われます.これは上記のコードの5行目9行目にも使われています.

2行目は,stdio.hというヘッダファイル(拡張子が.hのファイル)を取り込むための命令です.このように,#で始まる命令は,プリプロセッサ(=前処理指令)と呼ばれ,コンパイル作業に先だって実施されます.上記の例にある #includeというプリプロセッサは,その後に指定されたファイル(この場合は stdio.h)の内容で上記の #include から始まる行を置き換える指令です.よって上記の例では,コンパイル直前にstdio.h というファイルの中身が2行目と置き換わります.

ちなみにstdio.hというのはStandard I/Oの略であり,日本語では「標準入出力」と訳します.I/OというのはInput, Outputの略です.このヘッダファイルはキーボードやファイルからデータを読み込んだり,ファイルや画面にデータを出力するために必要なマクロコマンドや関数プロトタイプ宣言が書かれています.データの入出力を必要とする際には,必ずこのヘッダファイルをインクルードしておく必要があります.

stdio.h 以外にも,C言語には用途ごとに様々なヘッダファイルが用意されています.また,自分だけのヘッダファイルを作ることも可能です.それらのヘッダファイルが必要な場合には,必要に応じてファイル冒頭部分に書き込みます.例えばsinexpなどの数学関数を使いたければ,数学関数関連のヘッダファイルであるmath.hを読み込みます.具体的には,

#include <math.h>

と書きます.

5行目では,関数内で用いるローカル変数を定義しています.ローカル変数はその変数が宣言されている直近の {}の間でだけ用いることができる変数です.例えばこのプログラムでは,ローカル変数 abmain関数の範囲を表す{}の間でだけ有効となります.このローカル変数と対をなす概念にグローバル変数があり,グローバル変数はソースコードの全ての{}の外側で,多くの場合,ソースコードの冒頭で宣言される変数です.ですが,現在のプログラミング・パラダイムにおいてはグローバル変数を乱用することはあまり良いとされていません.それは,変数の出どころがわかりづらく,特に大規模なプログラムにおいてバグを招きやすく,また,バグが起こった際のデバッグも難しくなるためです.グローバル変数の仕様は極力最小限に止めるように心がけましょう.

なお,旧来のC言語の形式であるC90では関数の最初でしか変数が宣言できませんでしたが,C99以降,この制限は緩和されており,現在は変数を使用する直前に宣言も行うのが良いとされています.ただ,ごくまれにC90以前のバージョンしか使えない環境なども存在するので,上記のローカル変数の宣言位置については頭の片隅に置いておくと良いでしょう.

10行目では割り算の計算結果を新しく宣言した変数 c に代入しています.ここで変数 cは単なるdouble型ではなく, const doubleと定数形式で宣言されていることに注目してください.const doubleはそれ以降,値が変更できないという点を覗き double 型と同じように動作します.ここで const をつけているのは変数 cがこれ以降変化する必要がないからです.このように変化しない変数に const をつけることは特にC言語やcで大規模なプログラムを書く時には非常に重要で,バグを減らしたりデバッグを容易にする効果があります.自身で大きなプログラムを作成する際には,できる限り定数を使うように心がけましょう.

ちなみに,C言語においてmainという名前の関数は特別な存在です.C言語のプログラムが実行されると,一番初めにmain関数が呼び出されて実行されるという約束になっています.そのため,通常のC言語プログラムには,必ずmain関数を用意してあげる必要があります.この main 関数はC言語においては戻り値の型を明示せずに

main() {
    ...
}

のように書くこともできますが,ここでは,より明確に戻り値の型である int を明記し,また main 関数の引数にはコマンドライン引数を受け取っています.引数の argc はプログラム名を含むコマンド全体をスペースで区切った時に,何個の連続した文字列が現れるかを示しており argv はスペースで区切った時に得られる文字列を argc 個分格納しています.例えば,

./a.out 1 2 hello

のようにしてプログラムを実行した時には argc = 4 であり argv./a.out, 1, 2, helloの4つの文字列を格納した配列になっています.

練習問題

3つのdouble型小数a, b, cをキーボードから入力し,最初の2数の差(a-b) と,3番目の値 c の値が等しいか否かを判断するプログラムを作成せよ.(キーボードから値を読み込むにはscanfを使う)

ヒント (プログラムの一部分の例):

double a, b, c;

if (a - b == c) {
   printf("a-b == c\n");
} else {
   printf("a-b != c\n");
}

プログラムが完成したら,例えば,次の3つの小数を入力し,動作結果を確認してください.

丸め誤差

上記の演習プログラムに,0.6,0.4,0.2の3つの数を入力すると,0.6-0.4 = 0.2 であるにも関わらず,a-b !=c という表示が出たと思います.これは,計算機の丸め誤差が原因です.

一般に計算機が扱うことのできる数字の桁数は限られています.そのため,(2進数において)小数点以下の桁数が多い数(端的な例としては無限小数)を,そのまま計算機内で扱うことはできません.そのような場合,計算機内では適当な桁で丸められて近似表現されます.これにより生ずる誤差を丸め誤差と呼んでいます.

例えば,一見きりの良い数に思える0.2は,2進数では 0.001100110011….という無限循環小数になります.そのため,計算機内に記憶される際に,適当な桁で切り捨てられて(もしくは,コンパイラ等の仕様によっては切り上げられて)記憶されます.そのために,上記の0.6-0.4 != 0.2などといった不思議な結果が生じます.

この例からもわかるように,丸め誤差に関して注意したいのは,10進数できれいな形で表される数が2進数では必ずしもきれいな形にならない,という点です.そのため,プログラマが任意に数を選べるときには,丸め誤差を少なくするために,2進数できれいに表される数,すなわち,$1 / 2^n$(あるいは,$2^n$)となる値を使った方が良い結果を生みます.たとえば,本授業の後半で出てくる数値積分では,積分区間を複数の微小区間に分割し,各区間の近似計算結果の和をとる,という計算が行いますが,このとき,区間分割数を$2^n$にとると,誤差を低減することができます.

もう一つ,丸め誤差を低減するポイントとしては,少しでも長い桁数を表現できる変数型を使うことが挙げられます.C言語には,小数を扱う変数型として floatdoubleが用意されています.これらを比較すると,double型の方がはるかに多い桁数を扱うことができます.

より厳密にはfloat型は32ビット(2の32乗種類の数が表現できる),double型は64ビット(2の64乗種類の数が表現できる)で表されていて,それぞれが表すことのできるの値の絶対値の範囲は次の表のようになっています.

  最小値 最大値
float型 $1.17 \times 10^{-38}$ $3.40 \times 10^{38}$
double型 $2.23 \times 10^{-308}$ $2.23 \times 10^{308}$

これらの最大値,最小値などは,標準ヘッダファイルfloat.hに記載されていますので,興味のある人はfloat.hを覗いてみてください.よって,プログラミングに慣れないうちには,

としておくと良いでしょう.ただ,依然として並列計算や低レベルなアセンブリ言語などを使う場合はfloat型を使うほうがdouble型を使うときの数倍高速に動作する場合もあり,float型を使う場面がまったくないではありません (昨今の深層学習のプログラムなどはfloat型で書かれていることが多い).

ちなみに,古いプログラムの教科書では「float型を使うと処理速度が上がる」と書かれている場合があり,処理速度を向上させたいときには float型を使うことが推奨されていたりします.これに関しては,古い計算機ではfloatを使った方が良いと言えるものの,数値演算プロセッサ(=小数演算を専門に行うチップ)が標準装備されている現在のパソコンでは,逆にdouble型を使った方が高速な場合もあります.良くあるやり方としては,

typedef Float double;

のように新しい型を定義しており,必要に応じて

typedef Float float;

と書き直すなどの方法が考えられます.本講義ではそれほど速度を重要視するプログラムは作成しないため,精度を重視して double で統一して説明をすることとします.

浮動小数の値比較

丸め誤差以外にも,数値演算には様々な誤差が伴います.そのように様々な誤差を伴う数値演算において,

if (a == 0.123) {
    ...
}

といった厳密な条件判定をすることは意味がありません (ただし0.0については,そのように書くことも多い).

演習1で見たように,このような条件判定は意味の無い間違った結果を返してきます.意味のある正しい条件判定をするためには,演算で発生する誤差の評価を行い,その誤差幅を考慮して条件判定をする必要があります.

たとえば,2数を比較する際には,単に

if (x == y) {
    ...
}

とするのではなく,誤差幅を適切に見積もった上で,

#include <math.h>
#define EPSILON	1.0e-12

/*** 途中省略 ***/

if (fabs(x - y) < EPSILON) {
    ...
}

とする必要があります.上記のコードでfabsは小数の絶対値を与える関数である (単なるabsも存在するが,こちらはint型用なので注意する).fabsを使用するにはmath.hをインクルードするのを忘れないこと.この例では,許容誤差幅としてEPSILONという値を定義し\(1.0 \times 10^{-12}\)を許容誤差幅として設定しています.

課題1 (配点: 2点)

以後,課題は完成したらチェックを受けること

練習問題で作成したプログラムが,上記説明に基づき正しい動作をするように改良し,その動作を確認せよ.ただし,各数値 a, b, cscanf を用いてターミナルから入力可能となるようにし,いくつかの値について正しく動作することを確かめること.

オーバーフロー,桁落ち誤差

復習: C言語における関数の定義

C言語のプログラムは,多数の小さな関数の集まりによって記述されます.一般に,プログラミング初心者のうちは,一つの関数に大量の処理を詰め込む傾向が見られ,百行を超えるような処理をmain()関数一つですませようとしている例にも良く出会います.

そのように一つの関数に多くの処理をつめこんだプログラムは,大変読みづらく,開発効率の低下にもつながります.プログラムを書くときにはmain()関数も含め全ての関数をなるべくコンパクトに記述するようにしてください.

一つの関数の長さの目安としては,20〜30行程度,どんなに長くとも100行以内に収めるのが良いでしょう.100行を越えるような長い関数は,処理単位ごとに複数の関数にわけてください.なお,逆に短い分にはどれだけ短くても構いません.中身が数行しか無いような関数も,実際には多く使われます.

さて(抽象的な表現ですが)一般的な関数の記述は次のようになります.

return-type function-name (parameter declarations, if any) {
	declarations and statements
}

return-type には,その関数が戻り値(その関数が持つ値)の型を記述します.戻り値が無い場合には void と書きます.

function-name の部分には関数の名称を parameter declarations の部分には,関数に与える引数(パラメータ)を定義します.引数(パラメータ)はいくつあっても良いので,必要な数だけ,引数宣言を並べることができます.引数が無い場合には単に()と書くか,引数がないことを明示して(void)と書きます.引数のスコープ(=参照可能範囲)はその関数の内部になります.

関数の中身にはその関数内部で用いる変数 (これをローカル変数と呼ぶ)の定義と,それらを用いた実際の演算処理であるstatementsを書いていきます.古いC言語の仕様ではローカル変数の宣言は必ず関数の先頭で行う必要がありましたが,C99というC言語の仕様以降は,ローカル変数が関数内のどこであっても宣言できるようになりました.

ローカル変数と対になる概念として,どの関数の内部でもないスコープに定義されたグローバル変数があります.ローカル変数の名前は,グローバル変数と重複してつけることができます.例えば,グローバル変数にnumという名前の変数があったとして,同じnumという名前の変数を引数やローカル変数として宣言することができます.

このように同じ名前の変数が複数ある場合,よりスコープが近い(狭い)ものが優先されます.グローバル変数とローカル変数の両者でnumという変数を宣言した場合には,その関数の中ではローカル変数が優先して参照されます.その場合,グローバル変数は参照できません.そのため,一般的にはグローバル変数と同じ名前のローカル変数を用意することは避けるほうがよいでしょう.なお,同一のスコープで,名前の同じ複数の変数を宣言することはできません.

では,実際に関数宣言の一例としてCombination (組み合わせの場合の数) \(_n C_r\) を求める関数 combination を考えてみましょう.この関数は,nrの2つの引数(パラメータ)を必要とします.これらの引数はいずれも整数であり,この関数が返す計算結果もまた整数です.よって,この関数の記述は次のようになります.

long combination(int n, int r) {
    ...
}

一般に,コンビネーションの値は大きな値をとることが多いので,ここでは関数の return-type として long 型を使っています.

プロトタイプ宣言

関数を使う場合にプロトタイプ宣言と呼ばれる関数の事前宣言を行うことがあります.次のサンプルプログラムを見てみましょう.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
#include <stdio.h>

long combination(int n, int r);  // 関数のプロトタイプ宣言

int main(int argc, char **argv) {
    int a, b;
    long c;

    scanf("%d, %d", &a, &b);
    c = combination(a, b);
    printf("%d", c);
}

long combination(int n, int r) {
    /* 省略 */
}

この例で3行目に記述されているのが,プロトタイプ宣言です.プロトタイプ宣言は,コンパイラに事前に関数の仕様を伝えることで,エラーの検出を容易にするためのものです.一般にプロトタイプ宣言は,プログラムの冒頭部分(グローバル変数の宣言を書く部分と一緒)に書かれるのが通例で,C99以降の仕様ではプロトタイプ宣言が行われている箇所以降で,この関数を定義・仕様することができるようになります.

プロトタイプ宣言の書き方は,上記のサンプルにあるとおり,関数記述の1行目(15行目)をそのまま書き写すだけです.ただし,関数宣言部とプロトタイプ宣言とでは,セミコロンの有無が違います.15行目の行末には,セミコロン(;)をつけてはいけませんが,プロトタイプ宣言の行末には,セミコロン(;)が必ず必要です

関数宣言の次は,関数の呼び出し方(=使い方)を見てみましょう.関数の使い方は,とても簡単です.上記プログラムでは,11行目が,関数 combination() を呼び出している場所です.このように,関数を呼び出すには関数名を書いて,それに続く () の中に,プロトタイプ宣言に宣言された通りの型・個数の引数を与えればOKです.引数が無い(=voidである)関数の場合は,()の中に何も書く必要はありません.ただし,()を省略することはできません.

関数が返値を持つ場合には(return-typeが void 以外の場合には),代入式の右辺に持ってきたり,他の関数の引数として用いることもできます.

/* 関数値を変数に代入する */
a = combination(3, 4);

/* 関数値を他の関数の引数として使う.この場合は,func1()の返値を,combination()の引数としている */
combination(func1(2), 3);

整数型とオーバーフロー

C言語には,整数用の型としてintが用意されています. intには,修飾子として

などを適用することができます.修飾子をつけた場合の変数型名は,short intなどのように修飾子と元の型をつなげたものになります.ただし,それだと長いので,通常は short intint の部分を省いて単にshortlong int は単に long と書き表します.これらの整数型のうち,数値計算では intlong を使うのが一般的です.

intlongが,それぞれ何ビットの記憶領域を使うかは,使う計算機やコンパイラによって異なります.一般的には,intは16ビットか32ビット,longは32ビットであることが多いです(各変数型の大きさや扱える値の最大値は,標準ヘッダファイル limits.h の中に記載されています).16ビットの場合は,-32768〜32767の範囲,32ビットの場合は,-2147483648〜2147483647の範囲の整数を扱うことができます.

演算結果が,これら各変数型の扱える範囲を超えてしまうと,意味の無い値が変数に格納されてしまいます.そのような現象をオーバーフローと呼び,それにより生じる誤差を桁あふれ誤差,あるいは,オーバーフロー誤差と呼びます.前記の丸め誤差の場合は,誤差の大きさは小さいので被害も少ないのですが,オーバーフローが起こると,全く異なる計算結果が出てきてしまうので,オーバーフローには注意が必要です.

通常,intはその計算機でもっとも扱いやすい大きさ,となっていることが多いため,オーバーフローの心配が無い場合には,intを使います.ただし,オーバーフローの恐れが有る場合には,long(あるいは,最近のコンパイラでは,より大きなlong longという型もあります)を使っておくのが安全でしょう.

なお,設計演習室のCygwin環境では,intlongもどちらも32ビットであるため,扱える数値の範囲は一緒です.よって,今回の環境ではintlongを使い分ける必要はありません.ただし,別なコンパイラや,別の計算機でコンパイルすると,ビット数が異なる可能性があるので注意してください.

誤差の評価方法

前述の丸め誤差以外に,数値計算における主要な誤差には,次のものがあります.

それぞれの誤差について説明する前に,まずは誤差の評価方法について考えてみましょう.誤差の評価には,大きく分けて絶対誤差相対誤差という二つの尺度があります.例えば,実数 \(x\) の近似値を \(\tilde{x}\) としたとき,

と表されます.

例えば1.23456を小数第4位以下を切り捨てて,1.234とした場合の絶対誤差は,0.00056であり,相対誤差は,0.00056/1.23456 です.

桁落ち誤差

値がほぼ等しい2つの数の差を作ると,計算結果において絶対値が小さくなった分だけ相対誤差が大きくなってしまいます.この要因により増加する誤差を桁落ち誤差と呼びます.

たとえば,十進数で 1.23456と1.23432という二つの数を考えます.この二つの数が丸め誤差を含んでいるとする(二つの数がそれぞれ真の値から小数第6位で四捨五入されている)とすると,丸め誤差による絶対誤差の限界は$5.0\times 10^6$と評価できます.

両者を引き算すると,0.00024 が得られます.元の誤差範囲を考えると,ここでの絶対誤差の限界は$1.0 \times 10^{-5}$と評価できます.元の数と比較すると,絶対誤差は2倍になっただけです.ところが,相対誤差を考えてみると,引き算により値そのものが小さくなったため,元の数に比べ,約4桁も増加してしまいます.言い換えると,有効数字の桁数が4桁も減っています.有効数字桁数の減少は,その後の演算において更なる誤差の増加につながります.

このような桁落ち誤差を防ぐには,桁落ちを起こす引き算が生じないように式や計算方法を変形することが重要です.たとえば,十分小さいxに対して,sqrt(1+x)-1を計算する場合,式を変形してx/(sqrt(1+x)+1)とすることで,桁落ち誤差を防ぐことができます.

情報落ち誤差

一方,絶対値の大きく異なる2数の加減算では,小さい方の値が事実上無視されてしまうことがあります.これを情報落ち誤差と呼びます.このようなことが起こるのは,数値計算において,値が丸められることにより計算後の値が実際の値とずれてしまうことに起因します.

数値計算の分野では1.0に対して小さな小数$\varepsilon$を足し算した時に,その値が正しく計算できる最小の$\varepsilon$を計算機イプシロン (machine epsilon)と呼びます.例えば,C言語のfloatの場合には計算機イプシロンはおよそ$1.19 \times 10^{-7}$となっており,1.0にこれより小さな値を足し算したときには計算結果が変わらず1.0になってしまいます.

これは相対的な大小関係が維持されている場合もおよそ同様で,加算される二つの値の比が機械イプシロンを下回る場合には加算後に小さい方の値が無視されるけっかとなってしまいます.

このようなことが起こらないように,加減算を繰り返す場合には,値の大きさに注意を払う必要があります.例えば,単調減少する級数の和を求めるときに値の大きいものから順に足していくと,後半の小さい数は,それまでの和に比べて遙かに小さい値となるため,情報落ち誤差が顕著に発生します(後半の演算は,実質的に無視されてしまいます).

その反対に値の小さなものから順に足し算をしていけば,新しく足される大きな値とこれまでの値の和の比が計算機イプシロンを超えていれば(数値計算の範囲で)正しく値を求めることができます.したがって,そのような複数の数の和を求める場合には,値の小さいものから加算していくようにすることで情報落ち誤差を防ぐことが重要です

打ち切り誤差

無限級数の計算を有限項で打ち切った場合に生ずる誤差や,反復計算を有限回で打ち切ったときの誤差は,打ち切り誤差と呼ばれます.一例として$e^x$ を $x=0$ 近傍でテイラー展開(マクローリン展開)すると

\[e^x = 1 + \frac{x}{1!} + \frac{x^2}{2!} + \frac{x^3}{3!} + \cdots + \frac{x^{k-1}}{(k-1)!} + \frac{x^k}{k!} + \cdots\]

となります.この式を用いて \(e^x\) を計算することを考えてみましょう.この級数は無限級数ですが,計算機では無限の繰り返しを扱うことはできませんので,計算するにはどこか適当な項で計算を打ち切る必要があります.当然ながら,途中で打ち切った計算結果は,真の値と比べて誤差を持ちます.これが打ち切り誤差です.

このような反復計算では,誤差の範囲を明確にするため,許容する誤差範囲(下の例では $\varepsilon$ )をあらかじめ定めておいて,演算結果が誤差範囲に収まった時点で計算を打ち切るようにプログラムを書きます(ただし,収束しない可能性もある際には,無限ループに陥らないよう,反復回数に対しても制約を与えておくことが必要).

\[\vert S_k - S_{k-1} \vert < \varepsilon\]

ただし $S_k$は,第$k$項までの和を表します.

課題2 (配点: 4点)

上記のテイラー展開を用いて \(e^x\) を計算する関数 myexp() を作成せよ (myexp = my exponential function).この際,自作の myexpmath.h に定義されている exp() 関数で計算した真値との誤差が \(1.0 \times 10^{-12}\) 以下になるように関数を作成せよ. なお,自作関数の計算誤差を見積もるために myexp 内で math.hexp() を使ってはならない.

ヒント

$x$として負数を与えた際,そのまま計算してしまうと大きな桁落ち誤差がでる.桁落ち誤差が小さくなるように,$x$が負の場合の計算式を変形せよ.計算式をほんの少し変えるだけで誤差が圧倒的に小さくなる.

関数が完成したら \(e^{40}\) ならびに \(e^{-40}\) を計算し,math.hで定義されたexp()の計算結果とそれを比較せよ.その際,以下のテンプレートプログラムを参考にすると良い.

myexp.c

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define EPS 1.0e-12

double myexp(double x);

int main(int argc, char **argv) {    
    printf("exp(x)の値を計算します。\n");
    printf("xの値を入力してください: ");

    double x;
    if (scanf("%lf", &x) != 1) {
        printf("[ ERROR ] 入力された値は数値では内容です。\n");
        return 1;
    }
    printf("\n");

    double a = myexp(x);
    double b = exp(x);
    double rel = fabs(a - b) / fabs(b);
    printf("あなたの値: %.8e\n", a);
    printf("正しい値: %.8e\n", b);
    printf("相対誤差: %.8e\n", rel);
    printf("あなたの計算は「%s」です\n", rel < EPS ? "正しい" : "間違い");
}

double myexp(double x) {
    // この関数を実装してください
    return 0.0;
}

なお,実際には以下のような値となる.

ただし,プログラムの書き方で微妙に異なった値になる場合があるので注意すること.

課題3 (配点: 4点)

テイラー展開を用いて \(\log x\) を求めるプログラムを作成せよ.

なお少し考えると$\log x$はマクローリン展開できないため,そのまま式を評価するのでは,級数により値を求めることは難しいと考えられる.どのような関数を変わりに評価するのが良いか各自考えてほしい.

なお,正解の例は

\[\begin{aligned} \log(0.5) &= -6.93147181 \times 10^{-1} \\ \log(50) &= 3.91202301 \end{aligned}\]

のようになる.なお,\(x = 100000\)のような大きい値を入力すると,誤差が \(10^{-12}\) を少しだけ上回るので,比較的小さな \(x\) にだけ正しい値が計算できていれば正解とする.

この問題には2段階でヒントをつけたので,自分のレベルに合わせて適宜ヒントを参照せよ.

ヒント1

$\log x$は0周りでのテイラー展開(=マクローリン展開) ができないため,その代わりに$\log (1 + x)$をマクローリン展開しておき,実際に $\log x$を求めたいときには$x - 1$をマクローリン展開の計算に用いれば良いように思える.

ところが$\log (1 + x)$のマクローリン展開は,

\[\log(1 + x) = \frac{x}{1} - \frac{x^2}{2} + \cdots + (-1)^{n+1}\frac{x^n}{n} + \cdots\]

となるため,後から足される項に大きな値が現れて,打ち切りをどこですれば良いのかの判断が非常に難しくなる.実際に良い制度で求まる項数が存在したとしても,もし1つ項を少なくすれば,誤差はとても大きいと予想されるし,保守的に項を多めにとると,項の大きさが発散してオーバーフローを引き起こすためである.

この問題を解決するための方法は,やはり別の関数をマクローリン展開することであり,この場合には$\log \left( \frac{1 + x}{1 - x} \right) = \log(1+x) - \log(1-x)$の各項をマクローリン展開して,さらに変形すると良い.

ヒント2

前述の通り$\log \left( \frac{1 + x}{1 - x} \right) = \log(1 + x) - \log(1 - x)$の各項をマクローリン展開すると,それぞれ

\[\begin{aligned} \log(1 + x) &= +\frac{x}{1} - \frac{x^2}{2} + \frac{x^3}{3} - \frac{x^4}{4} + \frac{x^5}{5} - \cdots \\ \log(1 - x) &= -\frac{x}{1} - \frac{x^2}{2} - \frac{x^3}{3} - \frac{x^4}{4} - \frac{x^5}{5} - \cdots \end{aligned}\]

のようになる.したがって,

\[\log \left( \frac{1 + x}{1 - x} \right) = 2 \left( \frac{x}{1} + \frac{x^3}{3} + \frac{x^5}{5} + \cdots \right)\]

今,求めたいのは$\log x$の値であったので,仮に$y = \frac{1 + x}{1 - x}$と置き直せば,$x = \frac{y - 1}{y + 1}$となるので,

\[\log y = 2 \left( \frac{1}{1} \left( \frac{y - 1}{y + 1} \right)^{1} + \frac{1}{3} \left( \frac{y - 1}{y + 1} \right)^{3} + \cdots + \frac{1}{2n + 1} \left( \frac{y - 1}{y + 1} \right)^{2n + 1} + \cdots \right)\]

と書き直せる.この級数の各項は,$y > 0$より$n$が大きくなるに従って0に収束していくことが分かるので,この公式を用いることで,数値安定的に$\log x$の値を求めることができる.

次回までの予習

次回の授業までにC言語の配列とポインタについても各自復習しておいてください.