0. はじめに
ここでは Excel 数値計算ソフトウェア XLPack の基本的な使い方を説明します. 例題ワークシートをダウンロードすることができますので, 実際に使ってみながら説明を読むことをお勧めします.
本章では全体の概要を説明します. 1 章では例題ワークシートの入手方法と使い方を説明します. 2 章以降では分野ごとに例題を用いて使い方を説明します. 例題はすべて XLPack Basic (6.0) を使って実行可能です.
なお, 解法の詳細についてはここでは触れませんが, 別途「数値計算解法編」で解説していく予定ですのでそちらもご覧ください.
XLPackの構成
XLPackは, ワークシート関数ライブラリ, ソルバーおよびVBAサブルーチン／関数ライブラリの３つで構成されています. これらは, プログラミングをせずに数値計算するためのツールとプログラミングをサポートするツールの２つに分けられます.
プログラミングせずに数値計算
ワークシート関数またはソルバーを使うとプログラミングを必要とせずに計算を行うことができます.
ワークシート関数ライブラリ
ワークシート関数はExcelワークシートに直接記述できる関数です. 大きく分けると, 値を１つ返すものと, ベクトルあるいは行列のようにひとまとまりの値を返すものの２種類があります.
前者の主な例は特殊関数で, 通常のExcel関数のように数式の一部として使います. 他にはデータ入力の数値積分関数があります.
後者は通常のExcel関数とは違い配列型の結果を返す関数で, 主なものは線形計算(連立一次方程式, 固有値・固有ベクトル, 最小二乗法, 特異値など)の関数です. 他に, 補間, 代数方程式およびFFTがこのような関数となります. これらをワークシートに入力するには, 結果を出力するセル領域を選択し, 数式バーのfxをクリックするか, 「Shift」+「F3」のショートカットを押すと「関数の挿入」が表示されるのでこれを使って関数を選択します. パラメータを設定し終わったら, 「OK」をクリックしたり, 「Enter」を押したりはせずに, 「Ctrl」+「Shift」+「Enter」を押します.
XLPackソルバー
Excelに標準添付されているソルバーのように, メニュー操作で指定されたセルを参照することにより計算をすすめるアドインです. 数式とそれに必要なデータをワークシートに入力するだけで計算を行うことができます.
対象となるのは非線形方程式, 非線形最適化, 非線形最小二乗法, 数値積分および常微分方程式です.
VBA数値計算プログラミングをサポート
VBAサブルーチン／関数ライブラリ
VBAサブルーチン／関数ライブラリはVBAプログラムから呼び出すことができるサブルーチン／関数群です. VBAによる数値計算プログラミングが必要なときにこれらを使うことにより, 数値計算の専門的なプログラミングはライブラリに任せて, 主に入出力部分だけを作ればいいようになっています.
サンプルワークシートをひな形として活用すれば容易に使うことができます. また, マニュアルには簡単なものを除きほとんどのサブルーチン／関数の使用例が掲載されているので, これを参考にすることもできます.
対象とする問題
XLPackは一般的な数値計算のほぼすべての分野をカバーしています. 以下それぞれについて簡単に説明します.
特殊関数
連立一次方程式
固有値・固有ベクトル
近似
非線形方程式
非線形最適化
数値積分
常微分方程式
高速フーリエ変換 (FFT)
注 – 偏微分方程式およびその関連分野(大規模疎行列が係数の連立一次方程式など)については開発中です. また, 統計計算は今のところ対象外としています.
特殊関数
ベッセル関数, エアリー関数, 楕円積分, 指数および対数積分, ガンマ関数, ベータ関数, ディガンマ関数, 誤差関数などの特殊関数はワークシート関数およびVBA関数が提供されます. これらは式の中で組み込み関数のように使うことができます.
連立一次方程式
例えば未知数が2つの場合(2元連立一次方程式)は次のようなものでした. x1とx2が未知数です.
a11x1 + a12x2 = b1 a21x1 + a22x2 = b2
これを行列を使って表すと次のようになります.
Ax = b
ただし, A = (a11 a12)
(a21 a22)
b = (b1)
(b2)
x = (x1)
(x2)
A を係数行列, b を右辺ベクトル, x を解ベクトルといいます.
A の逆行列A-1がわかれば, 左から方程式の両辺に掛けてやると,
A-1Ax = A-1b
となり, A-1A が消えて次のようにxを求めることができます.
x = A-1b
ExcelではMINVERSEという逆行列を求める関数があるので, これを使って逆行列を求め右辺ベクトルbに左から掛けてやることにより解ベクトルxを求めることができます. しかし, 計算量が多いことと精度が悪くなりやすいことから, 逆行列を求めないで直接計算した方がよいとされています.
そのため, 逆行列を求めないで計算するガウスの消去法にもとづいた解法が広く使われています. 係数行列が一般行列の場合にはLU分解, 対称行列の場合にはコレスキー分解が標準的な解法となっています. まず係数行列Aを次のように分解します.
A = LU, Lは下三角行列, Uは上三角行列 -- LU分解 A = LDLT, Lは下三角行列, Dは対角行列 -- (修正)コレスキー分解
次にLU分解の場合ならば, 下のように y = Ux とおいて, Ly = b から y を求め, 最後に Ux = y から x を求めます.
LUx = b Ly = b Ux = y
ここで, L と U は三角行列なので, 代入操作だけで容易に解を計算することができます. 対称行列のコレスキー分解の場合にも同様にして解を計算することができます.
前者のサブルーチンはDgesv, 後者はDposvです. ワークシート関数はそれぞれ, WDgesv, WDposvです.
連立一次方程式を解く際には係数行列の条件数を求めておくと有益です.
条件数は次のように定義されます.
cond(A) = ||A||・||A-1||
条件数の値は, 最も条件が良い単位行列 I の場合 1, 最も条件が悪い特異行列の場合 ∞ で, 普通はその間の値となります. 条件数は定義どおりに計算すると逆行列を求めなければならず計算量が多くなるため, 通常は簡便に計算できる推定値が使わます.
これを使って, 求められた解xの相対誤差を次のようにして推定することができます.
||x - x||/||x|| ≦ cond(A)・ε
ここで, 太字のxは真の解, εは計算機イプシロン(*)です. 条件数が非常に大きい場合を悪条件であるといい, 条件数が1015程度以上になるとほとんど有効桁がない可能性があることがわかります.
ワークシート関数 WDgesvおよびWDposvでは条件数の推定値を返すようになっています. VBAプログラムの場合には, DlangeとDgeconおよびDlansyとDpoconを使って推定することができます.
(*) 計算機イプシロンは実数を計算機の浮動小数で近似したための誤差で, Excel (IEEE形式の倍精度) では1.11×10-16となります
固有値・固有ベクトル
行列Aに関する固有値問題とは,
Ax = λx, x ≠ 0
を満たすスカラーλ とベクトルx を求めるもので, λ を固有値, x を固有ベクトルといいます. Aがn×n行列の場合にはn個の固有値とそれぞれに対応する固有ベクトルを持ちます. 多くの応用において現れるのは対称行列の固有値を求める問題であるため, ここでは対称行列について考えます. Aが対称行列ならば固有値は実数になります.
固有値は次式を満たします.
det(A - λI) = 0
これは特性方程式と呼ばれ, λに関するn次代数方程式になります. そこで, この代数方程式を解いてλを求めればよいように思われますが, 数値計算ではもっと効率のよい方法が使われます.
まず行列の固有値を変えないようにして行列Aを三重対角行列 (対角要素とその隣の上下の要素だけが0でない行列) に変換します. この変換は有限回の演算で完了します. そして次に, 三重対角行列の固有値をQR法などの反復法により求めます. 三重対角行列だとこの反復法の過程の計算効率がよく, 最初から反復法を適用するよりも全体としても優ります.
この計算法によるVBAサブルーチンDsyevあるいはワークシート関数WDsyevを使うことにより, 対称行列のすべての固有値・固有ベクトルを求めることができます.
近似
近似式を求める問題では最小二乗近似がよく使われています. これは実験データへのモデル関数のあてはめと捉えることができます. モデル関数が線形 (ある基底関数の一次結合で表される) の場合は線形最小二乗法, そうでない場合は非線形最小二乗法で, それぞれ全く違う解法です. 最小二乗法と混同しやすいのが補間法ですが, こちらは飛び飛びに与えられた関数値からその間の値を求める (内挿する) ものです.
線形最小二乗法
m個の測定データ
(xi, yi) = (xi + exi, yi + eyi), i = 1, ..., m
が与えられているものとします. ここでそれぞれのデータ(xi, yi)は, 真値(xi, yi)に対して測定誤差(exiおよびeyi)を含んでいます.
各点の真値はある関数
yi = F(xi)
を満たすものとします.
この関数 F(x)を n個の基底関数φjの一次結合 f(x)により近似することにします.
F(x) ≒ f(x) = c1φ1(x) + c2φ2(x) + ... + cnφn(x)
この f(x) をモデル関数といいます. データ点が十分に多いとき(m ≫ n), モデル関数がデータをできるだけうまくあてはめるようなパラメータcj (j=1～n) を求めることにします.
モデル関数の基底関数φjとして最もよく使われるのは,
φj(x) = xj-1
すなわち, 多項式近似です.
f(x) = c1 + c2x + ... + cnxn-1
多項式近似以外でも一次結合で表すことができれば線形最小二乗法で扱うことができます. 例えば次のような場合です.
三角関数近似：
f(x) = c1sin(x) + c2sin(2x) + ... + cnsin(nx)
指数関数近似 (λjが定数の場合)：
f(x) = c1exp(λ1x) + c2exp(λ2x) + ... + cnexp(λnx)
さて, m個の測定データをモデル関数に代入して縦に並べると次の連立一次方程式になります (観測方程式といいます).
Ac = y
ただし, A は m×n係数行列で, Aij = φj(xi) (i = 1～m, j = 1～n)
c は n 解ベクトルで, 要素はパラメータcj (j=1～n)
y は m 右辺ベクトルで, 要素はデータyi (i=1～m)
この方程式は未知数の数よりも方程式の数の方が多い (m > n) ので厳密に解くことはできない (不能な方程式) ので, 近似解を求めることにします. ひとつの方法として, ノルム ||Ac – y|| の2乗が最小になる解を求める(*)方法が最小二乗法です. その場合, 解cを最小二乗解といいます.
最小二乗解は ||Ac – y||2の微分が0になる解として求めることができて, 次の方程式の解となります.
ATAc = ATy
これを正規方程式と呼び, 未知数の数 = 方程式の数 = n となるので解くことができます. ただし, 正規方程式の係数行列の条件数は Cond(ATA) = (Cond(A))2 と大きくなりやすく精度が悪くなることがあり, この解法はあまり使われていません.
代わってQR分解 (A = QR) が広く使われています. Qは m×m直交行列 (QTQ = Iとなる行列), Rは m×n上三角行列 (第n+1～m行は0) です. この分解ができれば元の連立一次方程式は次のように変形できます.
Ac = y QRc = y QTQRc = QTy Rc = QTy
最後の式において, 左辺の R は上三角行列なので代入操作だけで容易にc を求めることができます. こうして求めたcは最小二乗解になっています.
次に, パラメータ数n のとり方ですが, むやみに大きくすればよいというものではありません. 例えば, 明らかに直線にのっているデータを10次多項式で近似しても意味がありません. nが大きすぎると係数行列A の列が一次従属になってしまい, QR分解が失敗します. これをランク落ちといいます (ランク数は一次独立な列の数です). このような場合にはモデル関数の見直しが必要です.
この計算法によるVBAサブルーチンDgelsあるいはワークシート関数WDgelsを使うことにより, 精度よく最小二乗解を求めることができます.
(*) すなわち, 残差 ri = f(xi) – yi の二乗和 Σri2 (i = 1～m) を最小にするということです
非線形最小二乗法
線形最小二乗法と異なり, モデル関数がn個のパラメータを含む一般関数 f(x) の場合を考えます.
F(x) ≒ f(x; c1, c2, ... , cn)
モデル関数がデータをうまくあてはめるようにパラメータcjを決める方法のひとつが非線形最小二乗法で, 残差
ri = f(xi) - yi
の二乗和 Σri2 (i = 1～m) を最小にするようにパラメータcjを決めます. 線形最小二乗法の場合には連立一次方程式を解くことによりパラメータを求めましたが, 非線形最小二乗法では反復計算により解を求める必要があります. 初期値として与えられた近似解を出発点に解を改良する反復を行い, 最終的に最小二乗解に収束させます. 適切な初期値を与えないとうまく解が求まらない場合があります.
非線形最小二乗法には, マルカート法のVBAサブルーチン Lmdif1 を使うことができます. Lmdif1は微分 (ヤコビ行列) を前進差分近似で計算し, ユーザーによる微分計算を不要にした使いやすいサブルーチンです.
Lmdif1 はXLPackソルバーから使うこともできます.
補間
補間法はしばしば最小二乗法と混同されます. 最小二乗法ではデータ点に誤差が含まれていることを前提とし, 近似式はぴったりデータ点を通るわけではありません. しかし, 補間法は誤差が含まれていないデータ点を滑らかにつなぐ方法で, 補間式はぴったり全データ点を通ります. データ点に誤差が含まれていると不自然に曲がった近似式が得られます.
いろいろな関数値を表にした数表というのがあります. これを使う際に, 表にのっていない点の値を求めるために補間 (内挿) が使われていました. 現在では電卓やコンピュータがあるためその需要はなくなりましたが, 例えば, 計算するのに多大な労力が必要な関数の値を素早く求めたい場合などに使われます.
主な補間法としては以下があげられます.
- 多項式補間 (ラグランジュ補間)
n個の点を通る (n-1)次多項式を補間式として使う方法です
- 区分的多項式補間
全体を小さい区間に分割しそれぞれの区間で多項式補間を行う方法です. その中で最もよく使われるのは3次スプライン補間です. 連続するデータ点の組を3次多項式で近似し, つなぎ目が滑らかになるように2次までの導関数が連続になるように構成されます
3次スプライン補間係数と補間値を求めるためにはVBAサブルーチン Pchse と Pchfe あるいはワークシート関数 WPchse と WPchfe が使用できます.
非線形方程式
非線形方程式は一般に有限回の計算で解ける方法はなく, 初期値を与えてそれを改良して解に近づけていく反復法で解きます.
代数方程式
非線形関数が多項式の場合を代数方程式といいます.
p(x) = a0xn + a1xn-1 + ... + an-1x + an = 0
2, 3, 4次方程式については解の公式により解くことができます. 5次以上は反復法で解く必要があります. 代数方程式には多くの解法が提案されており使い分けるとよいとされています. なお, 実数解を1つだけ求めればよい場合には下の1変数の非線形方程式の解法を用いる方法もあります.
VBAサブルーチン Rpzero2 およびワークシート関数 WRpzero2 は, 実数係数の代数方程式をニュートン法で解くプログラムです. 複素解が求められます.
1変数の非線形方程式
一般の1変数の非線形関数の実数のゼロ点を求める問題です.
f(x) = 0
代表的な解法はニュートン法です. 収束が速く, これをベースにしている解法が数値計算には多くあります. 他には二分法などゼロ点を含む区間を縮小していく解法があり, 収束は速くありませんが確実な方法です.
VBAサブルーチン Dfzero は二分法と補間法を組み合わせたルーチンで, ゼロ点が含まれる範囲を与えると二分法よりも高速に解を見つけてくれます.
Dfzero はXLPackソルバーから使うこともできます.
多変数の連立非線形方程式
多変数の非線形連立方程式は数値計算としては解くのが難しい問題の一つです.
f1(x1, x2, ..., xn) = 0
f2(x1, x2, ..., xn) = 0
...
fn(x1, x2, ..., xn) = 0
線形連立一次方程式においては, 方程式が特異でない限り必ず唯一の解が存在しました. しかし, 非線形連立方程式では実数解が存在しないことがあり, また複数存在することもあります.
使われる解法はニュートン法 (の多変数への拡張) およびその修正版です. これらは与えられた初期値の近傍の解を見つけるもので, 大域的な解を求めるものではありません. そのため, 所望の解を得るためにはできるだけ近くの初期値を与える必要があります.
非線形連立方程式を解くためには, 修正版のひとつであるパウエルのハイブリッド法によるVBAサブルーチン Hybrd1 を使うことができます. Hybrd1 は微分(ヤコビ行列)を前進差分近似により計算する使いやすいルーチンです.
Hybrd1 はXLPackソルバーから使うこともできます.
非線形最適化
非線形関数の最小点 (符号を逆にすれば最大点) を求める問題です. これも反復法で解く必要があります.
1変数関数の非線形最適化
1変数の非線形関数の最小点を求める問題です.
min f(x)
よく使われる解法は黄金分割探索法のように最小点を含む区間を縮小していく方法です. 非線形方程式の二分法に似て収束は速くありませんが確実な方法です.
VBAサブルーチン Dfmin は黄金分割探索と放物線補間を組み合わせたルーチンで, 与えられた区間に含まれる最小点を見つけます. 非線形方程式の場合の Dfzero と似たサブルーチンです.
Dfmin はXLPackソルバーから使うこともできます.
多変数関数の非線形最適化
多変数の非線形関数の最小点を求める問題です. これも数値計算としては解くのが難しい問題の一つです.
min f(x1, x2, …, xn)
非線形最適化には種々の解法があります. シンプレックス法などの直接探索法は関数値だけを用います. 共役勾配法, 準ニュートン法などは関数値だけでなく微分値も用います. なお, これらの方法により求められるのは局所的最小点であり, 大域的最小点を求めることはできません. 局所的最小点はいくつも存在する可能性があり, 求めたい最小点に近い初期値を与えないと別の最小点に収束することがあります.
非線形最小二乗法や非線形連立方程式の解を求めるために残差二乗和を目的関数値として非線形最適化を適用することが考えられますが, それぞれの問題の特徴を使わないことになるので, やはりそれぞれの専用の解法を用いるべきです.
非線形最適化問題には, 準ニュートン法のVBAサブルーチン Optif0 を使うことができます. Optif0 は微分(ヤコビ行列)を差分近似により計算する使いやすいルーチンです.
Optif0 はXLPackソルバーから使うこともできます.
数値積分
f(x) は区間 [a, b]上で定義された関数とし, その定積分
∫ f(x) dx [a, b]
を数値的に求めます.
2種類の数値積分ルーチンがあります. 1つは関数入力タイプで, 被積分関数値が要求に応じて区間内の任意の点で計算できる場合に使われます. もう1つはデータ入力タイプで, 区間を包含するいくつかの点での被積分関数値が知られている場合に使われます.
関数入力タイプ
古くから使われている多くの求積公式があります. 例えば, 台形公式, シンプソン公式, ニュートン・コーツ公式, ガウス公式などです. しかし, これらの公式を使って得られた結果の誤差を推定するのは簡単ではありません. 最近では, 積分値と誤差の推定値の両方を計算するのに効果的な新しい公式が使われています. ガウス・クロンロッド公式はそのひとつです. ユーザーは推定誤差から結果の精度について有益な情報を得ることができます. また, 推定誤差は適応求積アルゴリズムを使う際に重要な役割を持っています. 適応求積ルーチンは, 推定誤差にもとづいて区間分割を自動的に調節し与えられた目的精度を満たすように積分値を計算します. 信頼性が高く効率的もよく, 汎用的に使うことができます.
VBAサブルーチン Qk15 は分点固定 (15点) のガウス・クロンロッド公式による求積ルーチンです. 15回の関数呼び出しで積分値とその推定誤差を計算します. VBAサブルーチン Qag はガウス・クロンロッド公式による適応求積ルーチンです. 被積分関数が滑らかな場合にはQk15で十分な精度が得られます. 被積分関数が滑らかでない場合や, 関数呼び出しが増えてでも所定の目的精度を得たい場合には, 適応求積分ルーチンQagを使うとよいかもしれません.
積分範囲の片方あるいは両方が有限でない場合にはVBAサブルーチン Qagi を使うことができます. これもガウス・クロンロッド公式による適応求積分ルーチンです.
Qag および Qagi はXLPackソルバーから使うこともできます.
データ入力タイプ
データ点の補間を行いその補間式の積分を計算することにより積分値を求めることができます. 3次スプライン補間を使ったVBAサブルーチン Pchse と Pchia, あるいは, ワークシート関数 WPchse と WPchia がこのために使用できます.
常微分方程式
常微分方程式は科学技術計算においてはしばしば現れます. 解析的に解けることは少なく数値解法に頼ることになります.
初期値問題
1階の常微分方程式は次の形に書くことができます.
y' = f(y, t)
常微分方程式の初期値問題とは, 初期点 t0 における値 y0 が与えられたとき, 最終点 tf における値 y(tf)を求めるものです. 次式は初期条件とよばれます.
y(t0) = y0
1階の連立方程式
1階の常微分方程式の初期値問題は, スカラーyをベクトルyに代えてこのまま連立方程式に拡張することができます.
y' = f(y, t) y(t0) = y0
例えば, 2変数の場合について書き下すと次のようになります.
y1' = f1(y1, y2, t) y2' = f2(y1, y2, t)
この場合, 初期条件は2つ必要になります.
y1(t0) = y1,0 y2(t0) = y2,0
高階の常微分方程式
高階の常微分方程式, すなわち, 2階以上の微分を含んだ常微分方程式は, 1階の連立常微分方程式の形にして扱うことができます. 例えば, 次の方程式を考えます.
y''' = f(y, y', y'', t)
これは, y1 = y, y2 = y’, y3 = y” とおけば, 次のように1階の連立常微分方程式になります.
y1' = y2 y2' = y3 y3' = f(y1, y2, y3, t)
この場合には次のように3つの初期条件が必要になります.
y(t0) = y0 y'(t0) = y'0 y''(t0) = y''0
同様にして高階の連立方程式を解くこともできます.
常微分方程式の初期値問題の解法
常微分方程式の初期値問題の解法には, オイラー法やルンゲ・クッタ法などの1段解法, アダムス法などの多段解法があります. 数値積分と同様に, 解を求めるだけではなく誤差の推定が効率よくできる種々の解法も開発されています. 例えば, 埋め込み型ルンゲ・クッタ法では1つの公式で2つの異なった次数の近似解を求めることができ, うち1つは誤差推定のために使われます. これにより, 推定誤差をもとに刻み幅を自動的に調節しながら解を求める適応ルーチンも作られています.
VBAサブルーチン Derkf は埋め込み型ルンゲ・クッタ法のひとつである4(5)次ルンゲ・クッタ・フェールベルグ法による適応ルーチンで, RKF45という名前で広く使われているものとほぼ同じです. 1階の連立常微分方程式に対応しています.
高速フーリエ変換 (FFT)
離散フーリエ変換 (Discrete Fourier Transform (DFT)) は, 離散データに対するフーリエ変換を数値的に求めるもので, ディジタル信号処理, 画像処理など多くの分野で使われています. DFTを少ない演算回数で求める数値計算法が高速フーリエ変換 (Fast Fourier Transform (FFT))です.
f(x)は周期2πの関数とします. 区間[0,2π]をN等分して, その分点xj=2πj/Nにおける関数値 f(xj)を fj と表します(j=0～N-1).
f(x)が実数関数の場合にはDFTは次式で表されます.
Fk = (1/2)(ak - i bk)
ただし, 係数 ak および bk は次のとおりです.
ak = (2/N)Σ fj cos(2πkj/N) (Σは j=0～N-1) (k=1～N2) bk = (2/N)Σ fj sin(2πkj/N) (Σは j=0～N-1) (k=1～N2) ただし, N2=N/2-1(Ｎが偶数の場合), N2=(N-1)/2(Nが奇数の場合)
ここで b0 = 0 となります. また, Nが偶数の場合 bN/2 = 0 となります. a0とaN/2は次のように1/2して表します.
a0 = (1/N)Σ fj (Σは j=0～N-1) aN/2 = (1/N)Σ fj cos(jπ) = (1/N)Σ(-1)j fj (Σは j=0～N-1)
f(x)が実数関数の場合には fj=conj(fj) が成り立ち, DFTの結果は Fk=conj(FN-k) となります. 従って, Fkはk=0～N/2についてだけ求めておけばよいことになります (conj()は共役複素数を表します).
DFTの逆変換を離散フーリエ逆変換(Inverse Discrete Fourier Transform (IDFT))とよび, 次式で表されます.
fj = Σ(ak cos(2πkj/N) + bk sin(2πkj/N)) (Σはk=0～N/2) (j=0～N-1)
注 – DFT/IDFTの定義において, 正規化係数(1/N)や符号は文献により異なることがあるので注意してください
VBAサブルーチン Rfft1f, Rfft1b, Rfft1i あるいはワークシート関数 WRfft1f, WRfft1b を使い1次元の実フーリエ変換および実フーリエ逆変換を求めることができます.
