環境システム科学演習−　MATLABを用いた海洋情報処理・作図演習

目次    授業ホーム    ホーム

２）　行列・連立1次方程式、最小自乗法(Least Square)

MATLABはその名の通り、Matrix Laboratoryであるので行列の演算が簡単に行える。ここでは、行列演算を用いて連立１次方程式を解く演習と、それに関連して最小自乗法を用いたデータの直線あるいは曲線近似を行う。

２-１　行列演算と連立１次方程式

目次    授業ホーム    ホーム

連立一次方程式は例えば以下のようなものである。

a₁₁ x₁ + a₁₂ x₂ + a₁₃ x₃ = b₁

a₂₁ x₁ + a₂₂ x₂ + a₂₃ x₃ = b₂

a₃₁ x₁ + a₃₂ x₂ + a₃₃ x₃ = b₃

これらは未知数3、方程式3なので条件を満たせば解くことが出来る。未知数ｘを消去しつつ求めていく方法で解くことが出来る。またこの問題は行列で表記すれば以下のようになる。

AX=b

ここで、

である。行列演算を用いてこの問題を解くにはAの逆行列を両辺の左からかければよい。即ち

A^-1AX=A^-1b

X=A^-1b

これによってXを求めることが出来る。Matlabでこれを行うにはMatlabの演算子の一つ\あるいは \で行う(同じキーボードキーであるがフォントが日本語なら　\、　英数字なら　\　になるので注意)。即ち、

>> X=A\b

あるいは

>> X=inv(A)*b

としても同様な結果が得られる。ここでinv()によってAの逆行列が求められた。もし逆行列が存在しない場合、問題を解くことが出来ない。

例えば以下の問をMATLABを使って解くことを考える。壁付近の流れを次元解析する。今、下図のように壁近傍の流れを考える。流体には粘性があるため流体粒子の一部は壁についたまま離れない。このため壁に近づけば近づくほど流れが弱まることが予想される。壁付近の運動を記述する関数を想像してみると、以下のようになるだろう。即ち、

F(U,ρ,τ,ν,z)=0

ここで、Uは壁に平行な流れ、ρは流体密度、τは壁に流体が及ぼす応力、νは分子粘性係数、zは壁からの距離である。関数Ｆの形はまだわからないが、運動はこれらの変数で記述できると予想した。今流れUと壁からの距離zを、他の変数ρ,τ,νを用いて無次元化することを考える。無次元化とは即ち単位を持たない変数にするということである。言い換えれば、ρ,τ,νを適当に用いてUとzの次元を表そうということである。ρ,τ,νは単位で書き表すと、ρ[kg/m³]、τ[N/m²]、ν[m²/s]である。ニュートン、Nは[kg m/s²]であるので、τ[kg/ms²]ともいえる。さて、kgは質量Ｍ、mは長さL、sは時間Tとして、以下のような次元行列を作ろう。

となる。この行列は行列式が0でないので逆行列を持つ。（Backinghamの定理はこの問題を変数3の作る時限行列が逆行列を持つのでランク3と呼び全変数の数5からランク数を引いて、5-3=2つの無次元変数が存在できることを述べている。）

3つの変数ρ,τ,νを使ってUの次元[LT^-1]を表すことを考える。それぞれを何乗して掛け合わせたり割ったりすればUと同じ単位、次元になるだろうか。この問題は以下のように記述できる。

U[M⁰LT^-1] = ρ[M¹L^-3T⁰]^x1τ[M¹L^-1T^-2]^x2ν[M⁰L²T^-1]^x3

次元 Mに関する式：　   x₁ + x₂          = 0

次元 Lに関する式：　-3x₁ - x₂ + 2x₃ = 1

次元 Tに関する式：　           -2x₂ -1x₃ = -1

今上に示した次元行列をAとし、未知数x₁,x₂,x₃をX、右辺をbとして、この連立方程式を解けばよい。つまり、

この連立方程式を解くプログラムを以下のように書き実行してみよう。

計算の結果x1は-0.5、x₂は0.5、x₃は0となる。つまりρ-1/2とτ1/2乗でUと同様な次元がつくれる事がわかる。Uを無次元化してえられる無次元変数はしたがって以下のとおりである。(8)の分母を物理学では摩擦速度と呼ぶ、u*で表すことが多い。

課題　次元解析

さて、上記の問題ではUをρ,τ,νを使って無次元化した。次は、zを無次元化してみよう。zを無次元化するプログラムを書き結果をまとめて提出しよう。

２-２　最小自乗法

M個の観測値が観測時刻tで存在するとしその観測値を格納する変数をyとする。yは(M×1)の行列で、MATLABでこれを列ベクトルと呼ぶため、以下ベクトルと呼ぶ。また以下断りがなければ大文字は行列を、小文字は列ベクトルを表すこととする。観測値をあるモデルで近似することを考えよう。即ち、

ここでθは観測値を表すモデルである。n(t)はモデルと実際の観測値との差、誤差である。まず(10)のようにモデルを直線(線形関数)と設定する。(10)式は、以下のように行列で書くことが出来る。

ここで誤差の自乗が最小となるとき最もそのモデルが忠実にデータを再現できると考えてよいだろう。ここで誤差nの自乗は以下のように書ける。

n^Tは転置を表し、行列やベクトルを転置する。また後述する行列の法則から以下が導ける。

この誤差Jを最小にするにはどうすればよいだろうか。（14）から誤差の2乗はxに関して2次式となっている。また2次の項である右辺第1項はE同士の掛けあわせであるから常に正の係数を持つことがわかる。2次の項が正の係数を持つ場合、その微分値が0となるｘで最小値が存在する。そこで(14)の微分を考える。

        ここまでの導出に用いた法則は以下のものである。ここでr、q、x、 yはベクトル。E、Aは行列でAは正方行列である。

誤差の微分(15)を0とすると以下の関係が得られる。

観測値をモデルとの誤差を最小にして再現するモデルの係数で構成されるベクトルｘは(18)を解くことで得られることがわかる。

では実際に最小自乗法を用いて直線近似をしてみよう。以下のサンプルプログラムを書き実行しよう。

ここでones(n,m)はn行m列の1を出力する。またhold onは一度描いた図をそのままにして次の作図を上に重ねるという命令だ。tは観測時刻、yは観測値である。Eはモデル行列である。t2は得たモデル係数ｘを使ってモデルで観測値を再現する為に準備したモデル用の時刻、y2はt2に対応するモデルの結果である。

課題　以下のデータは海洋観測によって得られた各水深における水温データである。水深をt、水温をyとし、水温の鉛直分布を3次式、y=a+bt+ct²+dt³　の曲線で近似せよ。