数値計算・数値解析 - ライブラリ使用のすすめ

数値計算とは？

数値計算はコンピュータを利用して物理現象等をより詳しく理解するための一つの手法です。コンピュータの誕生から現在にいたるまでコンピュータ利用の原点とも言えるもので、その重要性はコンピュータが高速化された現在、より一層高まってきていると言えます。

数値計算の大まかな流れ

数値計算を行いたい問題がある。
問題を数学的(方程式など)に表現する。
問題を解くための数値計算プログラムを作成する。
- 利用するアルゴリズムを考える、もしくは既に知られているアルゴリズムを選択する。
- プログラムを作成し、アルゴリズムを実装する。その際に数値計算ライブラリ(LAPACK, nAGライブラリ等)の利用も可能。
- コンパイル・リンクを行う。
実行して結果を得る。
結果を吟味する。（グラフ化、可視化なども利用）
問題がなければ完了。

主な数値計算のカテゴリ

数値計算には主に以下のようなカテゴリがあります。

偏微分方程式
常微分方程式
積分
連立一次方程式
固有値問題
方程式の根
フーリエ変換
曲線、局面フィッティング
最適化問題
乱数発生
回帰分析
分散分析
多変量解析
時系列解析
etc.

数値計算で利用される言語

「数値計算と言えばFortran」と言われるほど、この分野ではFortranが広く利用されています。これはFortranが数値計算プログラムを記述するために設計された言語であるという点や Fortran言語自体の歴史が古く過去の資産（研究者が過去に作成した計算コード）が豊富である等の理由によります。

しかしながら現代においてはアプリケーションに組み込まれて使われる事も多く、Fortran以外の言語（C言語、C++、Java、Ruby、Javascript等）での数値計算も盛んに行われるようになってきています。

数値計算と誤差

数値計算においてその結果には誤差が伴います。これはコンピュータ内部での数の表現が有限桁であることが原因で丸め誤差や桁落ちと呼ばれる誤差等が発生するためです。このため数値計算プログラムの作成を行う場合これらの誤差に注意してアルゴリズム選択や実装を行う必要があります。

丸め誤差の例

コンピュータ内部では、単精度を32ビット表し、倍精度を64ビットで表します。このようにコンピュータ内部での実数は、有限の桁数で表現されるため、真の値を正確に表現できない場合が多くあります。

また、内部では１０進数ではなく２進数で表現がなされているため以下のFortranコードの例のように十進数では仮数部１桁で表すことが可能な数字であっても内部表現は近似値となります。

real a
a = 0.1

と記述した場合に、aには0.1ではなく、単精度の浮動小数点表現で表すことができる0.1に最も近い値（0.100000001490116119384765625）が代入されることになります。これにより生じる誤差を丸め誤差と言います。

例えば以下のコード例では単精度の0.1を100回足して、その結果を出力しています。このコードを、ある特定の環境で実行した場合に10.0000019が結果として出力されました。

丸め誤差を示すサンプルコード（Fortran)

program rounding_error
  implicit none
  real :: sum = 0
  integer i
  do i = 1, 100
    sum = sum + 0.1
  end do
  print *, sum
end program rounding_error

出力例：
  10.0000019

情報落ちの例

絶対値の大きな数と絶対値の小さな数を加減算した際に、浮動小数点表現が有限の桁数であるため、絶対値の小さな数が結果に反映されないことで発生します。

例えば以下の例では100000と0.001を足してその結果を表示していますが、絶対値の小さい数字0.001が結果に反映されていません。

情報落ちのサンプル

program trailing_loss
  implicit none
  real :: a = 100000
  real :: b = 0.001
  print *, "a   = ", a, "       b   = ", b
  print *, "a+b = ", a+b
end program trailing_loss

出力例：
 a   =    1.0000000E+05        b   =    1.0000000E-03
 a+b =    1.0000000E+05

桁落ちの例

絶対値のほぼ等しい二つの数の差を求めたときに発生し、結果の有効桁数が小さくなります。

計算過程において上位の桁がゼロになった際に左詰が行われますが（その際下位の桁は0に設定される）、これにより相対的に大きな誤差が現れます。

例えば以下の例では 0.1234 - 0.1231 を単精度で計算させてその結果を表示しています。この例では 3.0000508E-04 となっています。

単精度の浮動小数表現で表すことができる、0.0003に最も近い値は 3.000000142492353916168212890625E-04 ですので、かなり大きな誤差であると言えます。

桁落ちのサンプル

program cancellation
  implicit none
  real a, b
  a = 0.1234
  b = 0.1231
  print *, a, b, a-b
end program cancellation

出力例：
   0.1234000   0.1231000   3.0000508E-04

打ち切り誤差の例

打ち切り誤差は計算をある時点で打ち切ってしまったことにより生じる誤差です。これは数値計算において、与えられた時間内に計算を終える事が求められるためで、計算を続ければより精度の高い結果が得られるような場合であっても計算を打ち切り、そこまでに求まった最善の値を結果とすることにより生じる誤差です。

例えばcos(x)をテイラー展開を用いて求める場合、次の無限級数を用いることが可能です。

これを下記のように有限個の項で計算することで打ち切り誤差が生じます。

数値計算プログラムに求められること

数値計算プログラムはまず第一に正しい計算結果を返す必要があります。そのためには上記で説明された「誤差」も考慮されている必要があります。

次に数値計算プログラムに求められることは計算スピードです。何度もパラメータを変えて計算させたい場合や、より厳密な解を求めたい場合など、スピードが早ければ早いほど数値計算の可能性が広がります。計算スピードは計算手順（アルゴリズム）や実装方法の工夫で高めることが可能ですが、ある程度以上のスピードアップには比較的高度な知識と経験が必要です。

更に現在ではマルチプロセッサ、マルチコアマシンが比較的安価に手に入りますので、このような環境で複数のスレッドを立ち上げて並列処理を行って計算スピードを高める事も可能です。コンパイラが自動的にある程度の並列化を行ってくれる部分もありますが、本当にパフォーマンスを得たい場合には手作業での並列化が必須となります。

数値計算ライブラリ利用のすすめ

数値計算プログラムを作成する際に、既に数値計算の専門家により作成された「数値計算ライブラリ」を利用する事も可能です。私たちnAGは数値計算ライブラリの利用を推奨いたします。その理由は以下のとおりです。

数値計算ライブラリを利用することでその部分のコードを自作しなくて良い。これにより結果の検証とバグとの戦いに費やされる時間と労力から解放される。
良質な数値計算ライブラリは計算結果の検証が十分に行われている。これにより計算結果の検証に費やす時間を削減することができる。
良質な数値計算ライブラリはアルゴリズム選択、実装共に優れている。そのためパフォーマンスにも優れている。
マルチコア、マルチプロセッサ向けに並列化済みの数値計算ライブラリを利用すれば、更にパフォーマンスが得られる。

このように数値計算ライブラリを利用するメリットは大いにあります。 nAGの数値計算ライブラリは、結果の正しさ、スピード、並列化の面で優れた科学技術計算・統計計算ライブラリです。

nAGが提供する数値計算ライブラリ