日本語プログラミング言語なでしこで数値計算（4）LU分解で連立方程式(NxN行列aとNベクトルbでax=bのxを求める）を解く。奥村晴彦さんのC言語による最新（標準）アルゴリズム事典を参考に。そのうち行列演算まとめてプラグインにします。

さて、(1)～(3)で以下のようなものをやった。

日本語プログラミング言語なでしこで数値計算（１）　まずは何はともあれ4段4次のルンゲクッタ法でローレンツ方程式を計算してみる。

日本語プログラミング言語なでしこで数値計算（2）　カオス（ローレンツ方程式）の次はフラクタルだということでマンデルブロ集合を描いてみる。

日本語プログラミング言語なでしこで数値計算（3）　コロナウイルスのような感染症流行を表す微分方程式、SIRモデルをルンゲクッタ法で計算する。西浦博さんの「感染症疫学のためのデータ分析入門」を参考に。　

今回は、本来は数値計算というならば最初にやった方がよかった連立方程式を。。。

まあLU分解でピボット選択があるやつがいいだろう。一から書くより何かを移植したい。

手元にあるのは奥村晴彦さんのC言語による最新アルゴリズム事典か（改訂版では標準アルゴリズム事典になっていた）

Numerical Recipes in Cか

Rosetta codeか。

https://rosettacode.org/wiki/LU_decomposition

まあ奥村先生のにしよう。

リンクはこちら。

https://n3s.nadesi.com/id.php?1890

こんな感じで計算できます。

ソースコード：

テキストでもコード書いておきます。そのうち行列の演算（四則演算、逆行列、行列式）とかまとめてプラグインにします。

#LU分解で連立方程式を計算する。
#参考文献：C言語による最新アルゴリズム事典/標準アルゴリズム事典（奥村晴彦さん）
# ax=bを計算する。a: NxN行列, b: N成分ベクトル
a = [[2,  5,  7],
　　 [4, 13, 20],
     [8, 29, 50]]
b = [23, 
     58, 
     132]

xにaとbの連立方程式答えを代入する。
xをベクトル表示する。

●(aとbの)連立方程式答えとは
　　n=aの表行数
　　weightにnのベクトル作成を代入する。
　　ipにnのベクトル作成を代入する。
　　xにnのベクトル作成を代入する。
　　det = 0
　　kを0からn-1まで繰り返す。
    　　ip[k]=k
    　　u=0
    　　jを0からn-1まで繰り返す。
        　　t = ABS(a[k][j])
       　　もし 、(t>u)ならば
　　      　　u=t
        　 ここまで
        ここまで  
        weight[k] = 1/u
    ここまで

    det=1
    kを0からn-1まで繰り返す。
        u=-1
        iをkからn-1まで繰り返す。
　         ii = ip[i]
           t = ABS(a[ii][k])*weight[ii]
           もし、 (t>u) ならば
               u = t
               j = i
           ここまで
        ここまで。
        ik = ip[j]
        もし、(j<>k) ならば
　　　    ip[j] = ip[k]
          ip[k] = ik
          det = -det
        ここまで
　      u=a[ik][k]
        det = det * u
        iをk+1からn-1まで繰り返す。
　　　     もし、i<=n-1 かつ i<>kならば
　             ii = ip[i]
               a[ii][k] = a[ii][k]/u
               t = a[ii][k] 
　　　      ここまで
            jをk+1からn-1まで繰り返す。
　　　　       もし、j<=n-1 かつ j<>kならば
　　　　          a[ii][j] = a[ii][j] - t*a[ik][j]
　　　　       ここまで
　　　      ここまで。
　　    ここまで
    ここまで　

    iを0からn-1まで繰り返す。
　     ii = ip[i]
       t=b[ii]
      jを0からi-1まで繰り返す。
　　     もし、i<>j  かつ j >=0ならば
　　　      t = t-a[ii][j]*x[j]
         ここまで
　    ここまで。
      x[i] = t
   ここまで。

   iをn-1から0まで1ずつ減らし繰り返す
　　   t=x[i]
       ii = ip[i]
       jをi+1からn-1まで繰り返す。
　　　　  もし、j<=n-1 かつ j<>iならば
　　　　　　　t = t-a[ii][j]*x[j]　
          ここまで
　　   ここまで
　　   x[i] = t/a[ii][i]
   ここまで
   xを戻す
ここまで


●(NとMの）行列作成とは
   A=[]
   iを０からN-１まで繰り返す
     A[i] = []
     jを０からM-1まで繰り返す
　　　A[i][j]=0
　　ここまで。
　ここまで。
   Aを戻す。
ここまで

●(Nの）ベクトル作成とは
   A=[]
   iを０からN-１まで繰り返す
     A[i] = 0
　ここまで。
   Aを戻す。
ここまで

●(Aを）行列表示とは
   iを０からAの表行数-１まで繰り返す
     jを０からAの表列数-1まで繰り返す
      「{A[i][j]}」と「　」を連続無改行表示
     ここまで
　　「」を表示
　ここまで
ここまで。

●(xを)ベクトル表示とは
　　iを0からxの配列要素数-1まで繰り返す。
       「{x[i]}」を表示
　　ここまで
ここまで。