« 松屋で「チーズボロネーゼコンボ牛めし」(ご飯大盛)をいただく。なんか懐かしい、、、と思ったら子供のころ実家でミートソースが余ったときにご飯にかけて食べた味そのものだ。 | トップページ | 「赤ずきん、旅の途中で死体と出会う。」を読んだ。赤ずきんちゃんが探偵役となっていろいろな童話の世界で起きる殺人事件を解決するお話。しかし旅の目的は…? »

2022年10月21日 (金)

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

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

日本語プログラミング言語なでしこで数値計算(1) まずは何はともあれ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

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

Nadeshiko_lu01

ソースコード:

Nadeshiko_lu02

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

 


#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を0からN-1まで繰り返す
     A[i] = []
     jを0からM-1まで繰り返す
   A[i][j]=0
  ここまで。
 ここまで。
   Aを戻す。
ここまで

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

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

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



« 松屋で「チーズボロネーゼコンボ牛めし」(ご飯大盛)をいただく。なんか懐かしい、、、と思ったら子供のころ実家でミートソースが余ったときにご飯にかけて食べた味そのものだ。 | トップページ | 「赤ずきん、旅の途中で死体と出会う。」を読んだ。赤ずきんちゃんが探偵役となっていろいろな童話の世界で起きる殺人事件を解決するお話。しかし旅の目的は…? »

パソコン・インターネット」カテゴリの記事

学問・資格」カテゴリの記事

日記・コラム・つぶやき」カテゴリの記事

コメント

コメントを書く

(ウェブ上には掲載しません)

« 松屋で「チーズボロネーゼコンボ牛めし」(ご飯大盛)をいただく。なんか懐かしい、、、と思ったら子供のころ実家でミートソースが余ったときにご飯にかけて食べた味そのものだ。 | トップページ | 「赤ずきん、旅の途中で死体と出会う。」を読んだ。赤ずきんちゃんが探偵役となっていろいろな童話の世界で起きる殺人事件を解決するお話。しかし旅の目的は…? »

最近の記事

最近のコメント

2022年12月
        1 2 3
4 5 6 7 8 9 10
11 12 13 14 15 16 17
18 19 20 21 22 23 24
25 26 27 28 29 30 31
フォト
無料ブログはココログ