日本語プログラミング言語なでしこで数値計算(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
こんな感じで計算できます。
ソースコード:
テキストでもコード書いておきます。そのうち行列の演算(四則演算、逆行列、行列式)とかまとめてプラグインにします。
#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]}」を表示
ここまで
ここまで。
« 松屋で「チーズボロネーゼコンボ牛めし」(ご飯大盛)をいただく。なんか懐かしい、、、と思ったら子供のころ実家でミートソースが余ったときにご飯にかけて食べた味そのものだ。 | トップページ | 「赤ずきん、旅の途中で死体と出会う。」を読んだ。赤ずきんちゃんが探偵役となっていろいろな童話の世界で起きる殺人事件を解決するお話。しかし旅の目的は…? »
「パソコン・インターネット」カテゴリの記事
- Interface2025年8月号Pythonで体験!はじめての暗号を買った。上杉暗号からRSA、AES、DHなど、特に楕円曲線暗号についてはコードも実際に動かすところまで詳しくかかれていた。耐量子暗号や聞いたことなかったY-00暗号や関数型暗号も記載。(2025.07.10)
- Gemini CLIが使えるようになっていたので早速VSCodeのターミナルから使って、JavaScriptで連立一次方程式を計算するコードを書いてもらった。普通にガウスの消去法で計算するhtmlを作ってくれた。(2025.06.27)
- Google ColabのJulia言語で1次元のGray-Scottモデル(∂u/∂t=u²v-(F+k)u+Du∂²u/∂x²,∂v/∂t=-u²v+F(1-v)+Dv∂²v/∂x²)を計算してパルスが次々分裂する様子を見る。空間6次の差分、時間8次のルンゲクッタ法で計算。(2025.07.08)
「学問・資格」カテゴリの記事
- 高周波・RFニュース 2025年7月13日 Pythonの高周波ライブラリscikit-rfがv1.8.0に、SamsungがGalaxy Z Fold7など発表→QualcommがSnapdragon 8 Eliteが使われていると発表、NGMNが基地局アンテナの推奨事項をまとめる、STMicroとMetalenzがメタサーフェス光学のライセンス締結(2025.07.14)
- Interface2025年8月号Pythonで体験!はじめての暗号を買った。上杉暗号からRSA、AES、DHなど、特に楕円曲線暗号についてはコードも実際に動かすところまで詳しくかかれていた。耐量子暗号や聞いたことなかったY-00暗号や関数型暗号も記載。(2025.07.10)
- 高周波・RFニュース 2025年7月8日 NordicとSercommのセルラーIoTモジュール、iFixitがFairphone 6を分解、スコアは10/10、RCR wireless newsのウェビナー2件(6GとIndustry4.0)、SEMCOが高耐圧C0G MLCCを車載急速充電に提案(2025.07.09)
- 高周波・RFニュース 2025年7月2日 5G Americasが6Gに向けセンシングと通信ホワイトペーパー発行、KYOCERA AVXが3dBハイブリッドカプラ発表、TDKが車載薄膜インダクタ発表、Nordicが1次電池向けPMIC発表、ローデ・シュワルツの6GとAI/ML解説記事(2025.07.02)
- 高周波・RFニュース 2025年6月30日 QualcommがAIを用いた6Rxアンテナ解説、Next G Allianceと日本のXGMFが5G,6Gで協力、5G Americasが25Q1で5G加入者増加と発表、TechInsigtsがHuawei Pura 80 Pro+分解、Qorvoが5-7GHzのWi-Fi 7 FEM発表(2025.06.30)
「日記・コラム・つぶやき」カテゴリの記事
- 高周波・RFニュース 2025年7月13日 Pythonの高周波ライブラリscikit-rfがv1.8.0に、SamsungがGalaxy Z Fold7など発表→QualcommがSnapdragon 8 Eliteが使われていると発表、NGMNが基地局アンテナの推奨事項をまとめる、STMicroとMetalenzがメタサーフェス光学のライセンス締結(2025.07.14)
- 高周波・RFニュース 2025年7月8日 NordicとSercommのセルラーIoTモジュール、iFixitがFairphone 6を分解、スコアは10/10、RCR wireless newsのウェビナー2件(6GとIndustry4.0)、SEMCOが高耐圧C0G MLCCを車載急速充電に提案(2025.07.09)
- 高周波・RFニュース 2025年7月2日 5G Americasが6Gに向けセンシングと通信ホワイトペーパー発行、KYOCERA AVXが3dBハイブリッドカプラ発表、TDKが車載薄膜インダクタ発表、Nordicが1次電池向けPMIC発表、ローデ・シュワルツの6GとAI/ML解説記事(2025.07.02)
- 高周波・RFニュース 2025年6月30日 QualcommがAIを用いた6Rxアンテナ解説、Next G Allianceと日本のXGMFが5G,6Gで協力、5G Americasが25Q1で5G加入者増加と発表、TechInsigtsがHuawei Pura 80 Pro+分解、Qorvoが5-7GHzのWi-Fi 7 FEM発表(2025.06.30)
- 高周波・RFニュース 2025年6月27日 Qualcommが6Gに向けての3GPPリリース20解説、TDKが100V1608サイズ1μFのMLCCを発表、Skyworksが低ジッタのクロックファミリーを発表、Elisa,Ericsson,MediaTekが5G SAで8Gbpsを達成(2025.06.27)
« 松屋で「チーズボロネーゼコンボ牛めし」(ご飯大盛)をいただく。なんか懐かしい、、、と思ったら子供のころ実家でミートソースが余ったときにご飯にかけて食べた味そのものだ。 | トップページ | 「赤ずきん、旅の途中で死体と出会う。」を読んだ。赤ずきんちゃんが探偵役となっていろいろな童話の世界で起きる殺人事件を解決するお話。しかし旅の目的は…? »
コメント