« 2022年9月 | トップページ | 2022年11月 »
Scratchのリンクはこちら:
Jack O' Lantern's Three body problem
動画にしたもの。クリックで始まります。シンプレクティックオイラー法は1次なんでちょっと8の字からずれていっているけれど。
シンプレクティックオイラー法についてはこちらを。
ドラマも始まってますが、medium読んでない方、ドラマの最終回を見てない方はこれは結末に触れられているのでまだ読まない方がいいです。mediumの衝撃にはどうしても比べてしまいますが、これも面白かった。
1本目の生者の言伝も倒叙ものですが、男子中学生を翻弄する翡翠さんが面白い。。。がそんな単純な結末ではもちろんなかった。
2本目の覗き窓の死角は、本格ミステリ好きの女性カメラマンでもある初めてできた友達がある殺人事件の容疑者になり、しかもそのアリバイは翡翠自身が証言することに(殺害時刻には翡翠が撮影されていた)。
トリック以上に、友達を疑うことになる翡翠さんが今までにない苦悩を。。。真さんも大ピンチに。トリックもかなり意外なもの。あれはちょっと出来すぎですが、、、
しかし真さんの「はわわてへぺろこっつんこ」で解決する、という言葉は笑った。
ドラマもまだ見てない(録画済み、まとめてみる)ですが、おそらく原作読んでない人はラストで衝撃を受けるだろう。
久しぶりに京都国立近代美術館へ。市民からの寄贈で始まっているというのがすごいな。
内容もバラエティに富んでいて面白かった。リキテンスタインの絵にちょっと衝撃を受けたり。
ロシアの美術も面白かった。もちろんピカソも。あと同時開催のコレクション展も面白い。特に日本の作品。
このハシビロコウだけは写真OKでした。
まずは各国のリニアスケール。日本、よく見ると変曲点があって増加に転じだしている。
あとフランス、ドイツも増えてる。イギリスは、、、もう統計取ってないのだろう。中国もゼロコロナ&どこまで正確に出してるかわからないしてるわりには増えてる。
次は各国のログログプロット。もうイギリス(何度も言いますがもうちゃんと計測されてない)と日本がほぼ同等になった。韓国は超えてるけど。
次は日本のワクチン接種者数累計と累積陽性者数。ワクチンが増えないなあ。私も4回目接種したけど。
韓国、日本、アメリカの陽性者数と人口比。日本はもう17.5%になった。
最後は日本の詳細ログログプロット。もう3年たつので経過1000日では足りなくなってきた。。。
さて、(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]}」を表示
ここまで
ここまで。
さて1回目と2回目はこういうのをやってみた。
日本語プログラミング言語なでしこで数値計算(1) まずは何はともあれ4段4次のルンゲクッタ法でローレンツ方程式を計算してみる。
日本語プログラミング言語なでしこで数値計算(2) カオス(ローレンツ方程式)の次はフラクタルだということでマンデルブロ集合を描いてみる。
今回はせっかくルンゲクッタ法の計算ルーチンを作ったので、このコロナ禍で有名になったSIRモデル
https://ja.wikipedia.org/wiki/SIR%E3%83%A2%E3%83%87%E3%83%AB
を計算してみよう。
説明はプログラムにも書いていますが
# コロナウイルスのような感染症流行についての微分方程式
# SIRモデル
# dS/dt = -βSI
# dI/dt = βSI-γI
# dR/dt = γI
# を4段4次のルンゲクッタ法で計算する。
# S: まだ未感染で免疫がない感受性人口
# I: 感染している感染性人口
# R: 回復した治癒人口
# N:全人口=S+I+R (1に規格化)
# β:感染伝達率, γ:回復率
で、これは西浦博さんの「感染症疫学のためのデータ分析入門」を参考にした。
ではプログラムのリンクはこちら:
結果の例はこちら。S,I,Rの初期値やβ、γのパラメータをいろいろいじって遊べる。
リストはこんな感じ。
テキストデータでのソースコード:
# コロナウイルスのような感染症流行についての微分方程式
# SIRモデル
# dS/dt = -βSI
# dI/dt = βSI-γI
# dR/dt = γI
# を4段4次のルンゲクッタ法で計算する。
# S: まだ未感染で免疫がない感受性人口
# I: 感染している感染性人口
# R: 回復した治癒人口
# N:全人口=S+I+R (1に規格化)
# β:感染伝達率, γ:回復率
全描画クリア。
# 画面設定
画面幅=描画中キャンバスの「width」をDOM属性取得。
画面高=描画中キャンバスの「height」をDOM属性取得。
「18px sans-serif」に描画フォント設定。
[80,20]に「SIRモデル(ルンゲクッタ法で計算)」を文字描画。
[80,40]に「S:赤, I:青, R:緑」を文字描画。
時間に0を代入する。
時間ステップに0.01を代入する。
最大時間に140を代入する。
x最大値に最大時間を代入する。
x最小値に-10を代入する。
y最大値に1.2を代入する。
y最小値に-0.2を代入する。
x最小値と0からx最大値と0まで黒色で1で直線プロット。
0とy最小値から0とy最大値まで黒色で1で直線プロット。
#初期値
S0に0.9999を代入する。
I0に0.0001を代入する。
R0に0を代入する。
ベータに0.3を代入する。
ガンマに0.1を代入する。
時間が最大時間以下の間
#ルンゲクッタ法のメイン計算
KS1にS0とI0とR0の関数fSを代入する。
KI1にS0とI0とR0の関数fIを代入する。
KR1にS0とI0とR0の関数fRを代入する。
KS2に(S0+時間ステップ×KS1/2)と(I0+時間ステップ×KI1/2)と(R0+時間ステップ×KR1/2)の関数fSを代入する。
KI2に(S0+時間ステップ×KS1/2)と(I0+時間ステップ×KI1/2)と(R0+時間ステップ×KR1/2)の関数fIを代入する。
KR2に(S0+時間ステップ×KS1/2)と(I0+時間ステップ×KI1/2)と(R0+時間ステップ×KR1/2)の関数fRを代入する。
KS3に(S0+時間ステップ×KS2/2)と(I0+時間ステップ×KI2/2)と(R0+時間ステップ×KR2/2)の関数fSを代入する。
KI3に(S0+時間ステップ×KS2/2)と(I0+時間ステップ×KI2/2)と(R0+時間ステップ×KR2/2)の関数fIを代入する。
KR3に(S0+時間ステップ×KS2/2)と(I0+時間ステップ×KI2/2)と(R0+時間ステップ×KR2/2)の関数fRを代入する。
KS4に(S0+時間ステップ×KS3)と(I0+時間ステップ×KI3)と(R0+時間ステップ×KR3)の関数fSを代入する。
KI4に(S0+時間ステップ×KS3)と(I0+時間ステップ×KI3)と(R0+時間ステップ×KR3)の関数fIを代入する。
KR4に(S0+時間ステップ×KS3)と(I0+時間ステップ×KI3)と(R0+時間ステップ×KR3)の関数fRを代入する。
S1 にS0 + 時間ステップ×(KS1+2*KS2+2*KS3+KS4)/6を代入する。
I1 にI0 + 時間ステップ×(KI1+2*KI2+2*KI3+KI4)/6を代入する。
R1 にR0 + 時間ステップ×(KR1+2*KR2+2*KR3+KR4)/6を代入する。
時間とS0から(時間+時間ステップ)とS1まで赤色で5で直線プロット。
時間とI0から(時間+時間ステップ)とI1まで青色で5で直線プロット。
時間とR0から(時間+時間ステップ)とR1まで緑色で5で直線プロット。
S0=S1
I0=I1
R0=R1
時間を時間ステップだけ増やす。
ここまで。
#SIRモデル
●(SとIとRの)関数fSとは
-ベータ×S×Iを戻すこと
ここまで
●(SとIとRの)関数fIとは
ベータ×S×I - ガンマ×Iを戻すこと
ここまで
●(SとIとRの)関数fRとは
ガンマ×Iを戻すこと
ここまで
●(xの)x座標変換とは
それ=画面幅 * (x - x最小値)/(x最大値 - x最小値)
ここまで
●(yの)y座標変換とは
それ=画面高 * (y - y最大値)/(y最小値 - y最大値)
ここまで
●(x1とy1からx2とy2までcでwで)直線プロットとは
cに線色設定
wに線太設定
[x1のx座標変換, y1のy座標変換]から[x2のx座標変換, y2のy座標変換]まで線描画
ここまで
さて前回はルンゲクッタ法でローレンツ方程式を解いて図示してみた。
日本語プログラミング言語なでしこで数値計算(1) まずは何はともあれ4段4次のルンゲクッタ法でローレンツ方程式を計算してみる。
次は、、、まあローレンツ方程式に次いで有名なのはマンデルブロ集合じゃないだろか、ということでやってみた。
リンクはこちら:
こんな感じで描ける。
プログラムリストはこんな感じ。
リストをテキストでも書いておきます。
#マンデルブロ集合を描くプログラム
全描画クリア。
# 画面設定
画面幅=描画中キャンバスの「width」をDOM属性取得。
画面高=描画中キャンバスの「height」をDOM属性取得。
x最大値に1を代入する。
x最小値に-2.5を代入する。
y最大値に1.5を代入する。
y最小値に-1.5を代入する。
最大繰り返し回数に255を代入する。
iを0から画面幅まで1ずつ増やし繰り返す
jを0から画面高まで1ずつ増やし繰り返す
#複素数でz(n+1) = z(n) + cを反復計算し、収束するかどうかの判定をする。
#z(0)=0とする。c = cx + i*cy
cx=iのx座標変換
cy=jのy座標変換
繰り返し回数に0を代入する。
xに0を代入する。
yに0を代入する。
(x^2+y^2 < 4 かつ 繰り返し回数 < 最大繰り返し回数) の間
x1 に xを代入する。
y1 に yを代入する。
#z(n)^2 + c = (x(n)+i*y(n))^2 + cx + i*cy = x(n)^2 -y(n)^2 + cx + i*( 2*x(n)*y(n) + cy)
x に x1 ^ 2 - y1 ^ 2 + cx を代入する。
y に 2 × x1 × y1 + cy を代入する。
繰り返し回数を1増やす
ここまで。
1に線太設定
#色は適当なので、もっときれいな設定をしたほうがいいです。
色 = RGB(INT(LOG(繰り返し回数)/LOG(最大繰り返し回数)*最大繰り返し回数),INT(LOG(繰り返し回数)/LOG(最大繰り返し回数)*最大繰り返し回数),INT(LOG(繰り返し回数)/LOG(最大繰り返し回数)*最大繰り返し回数))
色に線色設定。
[i, j, 1, 1]へ四角描画。
ここまで。
ここまで。
●(iの)x座標変換とは
x最小値 + i * (x最大値- x最小値)/ 画面幅を戻すこと
ここまで
●(jの)y座標変換とは
y最大値 + j * (y最小値 - y最大値)/ 画面高 を戻すこと
ここまで
さて、次は何しようかな…
今月のIEEE Microwave Magazineの特集はマイクロ波磁気工学。
https://ieeexplore.ieee.org/xpl/mostRecentIssue.jsp?punumber=6668
アイソレータ、サーキュレータ、チューナブルフィルタなどいろいろな用途はあるが、私も昔こういうのに携わっていたので懐かしい!と思ったのはLandau-Lifshitz-Gilbert(LLG)方程式!あの教科書のランダウ・リフシッツです。減衰に関する項をギルバートさんが入れたという。
そしてMicrowave Journalの特集はパッシブと制御回路。
https://www.microwavejournal.com/publications/1
ミリ波のPCBに関する記事も面白いが、iPhone14が衛星緊急通信に使うGlobalstarという会社の記事が出てた!なんてタイムリー。
エリクソンの技術資料。E-band(70/80GHz)が大事とのこと。ほんとかな。。。
アナログデバイセズのプレスリリース。
個人的に注目しているKumuの記事。
リコーの360°カメラがu-bloxのGNSSモジュールを使う。
最近プログラム記事のネタ切れ感が半端ない、、、Visual C#を使い始めたがまあ仕事で使っている部分もあるのであまりかけない。
流行りのRUSTなどより、使っている人がとても少なそうな言語でいろいろ遊ぼうと思ったり。
そこで、日本語でプログラミングできる「なでしこ」の存在を思い出した。
おお、いつの間にかブラウザベースでプログラミングができるようになっているのか。これはちょうどいい。
しかも保存して公開できるし。誰もやってなさそうな数値計算ネタをやっていこう。
まずは何はともあれローレンツ方程式だ。
リンクはこちら:
https://n3s.nadesi.com/index.php?page=1875&action=show
こんな感じで日本語でプログラムして、
実行すると、、、おお!おなじみのグラフが描けた。
普通に何でもプログラムできそうだ(まあJavascriptに変換されているそうなのでJavascriptで出来ることは何でもできそう)。
ちょっと遊んで行ってみよう(続く)。マンデルブロ集合とかSIRもでるとかどうだろう。
テキストでのプログラム:
# ローレンツ方程式
# dx/dt に-σ(x-y)
# dy/dt に-xz + rz -y
# dz/dt にxy - bz
# を4段4次のルンゲクッタ法で計算する。
全描画クリア。
# 画面設定
画面幅=描画中キャンバスの「width」をDOM属性取得。
画面高=描画中キャンバスの「height」をDOM属性取得。
「24px sans-serif」に描画フォント設定。
[25,20]に「ローレンツ方程式(ルンゲクッタ法で計算)」を文字描画。
x最大値に25を代入する。
x最小値に-25を代入する。
y最大値に50を代入する。
y最小値に-30を代入する。
x最小値と0からx最大値と0まで黒色で直線プロット。
0とy最小値から0とy最大値まで黒色で直線プロット。
時間に0を代入する。
時間ステップに0.01を代入する。
#初期値
X0に0.1を代入する。
Y0に0.1を代入する。
Z0に0.1を代入する。
10000回繰り返す
#ルンゲクッタ法のメイン計算
KX1にX0とY0とZ0の関数fXを代入する。
KY1にX0とY0とZ0の関数fYを代入する。
KZ1にX0とY0とZ0の関数fZを代入する。
KX2に(X0+時間ステップ×KX1/2)と(Y0+時間ステップ×KY1/2)と(Z0+時間ステップ×KZ1/2)の関数fXを代入する。
KY2に(X0+時間ステップ×KX1/2)と(Y0+時間ステップ×KY1/2)と(Z0+時間ステップ×KZ1/2)の関数fYを代入する。
KZ2に(X0+時間ステップ×KX1/2)と(Y0+時間ステップ×KY1/2)と(Z0+時間ステップ×KZ1/2)の関数fZを代入する。
KX3に(X0+時間ステップ×KX2/2)と(Y0+時間ステップ×KY2/2)と(Z0+時間ステップ×KZ2/2)の関数fXを代入する。
KY3に(X0+時間ステップ×KX2/2)と(Y0+時間ステップ×KY2/2)と(Z0+時間ステップ×KZ2/2)の関数fYを代入する。
KZ3に(X0+時間ステップ×KX2/2)と(Y0+時間ステップ×KY2/2)と(Z0+時間ステップ×KZ2/2)の関数fZを代入する。
KX4に(X0+時間ステップ×KX3)と(Y0+時間ステップ×KY3)と(Z0+時間ステップ×KZ3)の関数fXを代入する。
KY4に(X0+時間ステップ×KX3)と(Y0+時間ステップ×KY3)と(Z0+時間ステップ×KZ3)の関数fYを代入する。
KZ4に(X0+時間ステップ×KX3)と(Y0+時間ステップ×KY3)と(Z0+時間ステップ×KZ3)の関数fZを代入する。
X1 にX0 + 時間ステップ×(KX1+2*KX2+2*KX3+KX4)/6を代入する。
Y1 にY0 + 時間ステップ×(KY1+2*KY2+2*KY3+KY4)/6を代入する。
Z1 にZ0 + 時間ステップ×(KZ1+2*KZ2+2*KZ3+KZ4)/6を代入する。
X0とY0からX1とY1まで赤色で直線プロット。
X0とZ0からX1とZ1まで青色で直線プロット。
X0=X1
Y0=Y1
Z0=Z1
ここまで。
#ローレンツ方程式
●(XとYとZの)関数fXとは
-10*(X-Y)を戻すこと
ここまで
●(XとYとZの)関数fYとは
-Y-X*Z+26*Xを戻すこと
ここまで
●(XとYとZの)関数fZとは
X*Y - 8*Z/3を戻すこと
ここまで
●(xの)x座標変換とは
それ=画面幅 * (x - x最小値)/(x最大値 - x最小値)
ここまで
●(yの)y座標変換とは
それ=画面高 * (y - y最大値)/(y最小値 - y最大値)
ここまで
●(x1とy1からx2とy2までcで)直線プロットとは
cに線色設定
[x1のx座標変換, y1のy座標変換]から[x2のx座標変換, y2のy座標変換]まで線描画
ここまで
続きを読む "日本語プログラミング言語なでしこで数値計算(1) まずは何はともあれ4段4次のルンゲクッタ法でローレンツ方程式を計算してみる。" »
サイコロの絵から(ようやく、、、)あ、これ絵は挿絵じゃないんだ!ここもふくめて小説になってるんだと気づく(遅すぎ)。
そこから再度絵を見直していってさらに面白くなってきた。絵がなかったらクライマックスのシーンのイメージが全然違ったものになったと思うしこれはとてもよかった!
あらすじは「製菓会社の広報部で働く岸は、商品への異物混入トラブルに見舞われる。疲労困憊する岸だったが、ある議員の登場で状況が変わる。その議員はある夢について語り始める。。。」
というもの。現実と別の世界が夢で繋がる、というとよくある異世界ものみたいですが伊坂幸太郎さんが書くともっとひねりがきいていて
めちゃくちゃ面白い。またハシビロコウが出てきますがそれがイラストがとてもいいんですよ。
そしてラストは本当に大ピンチ。そのクライマックスもイラストをずっと見てるとさらに感慨深くなる。
本当に面白い作品なのでお勧めです。
今年のノーベル物理学賞はアラン・アスペ、ジョン・クラウザー、アントン・ツァイリンガーさんの3人が「量子もつれ、ベルの不等式の破れの確立、量子情報科学の開拓」 で受賞!
https://www.nobelprize.org/prizes/physics/2022/summary/
教科書に出てくるような話で、これは私も初めて聞いたときすごいなと思った。
で、毎年やっている素粒子・宇宙と、物性・基礎論がどう推移しているかというグラフ書いてみた。
昔は一年おきだったがどうなるか、というのがだいぶ前に気になったので。
結果はこちら。
なんとなく周期が1年おきから2年おきに変わったような気が、、、
最近、Visual C#で数値計算シリーズをやっている。でスミスチャートを描くライブラリも作ろうかと思ったが、、、
あのバックグラウンドにひかれているたくさんの円(R一定とX一定)の式なんだっけ、、、と思うが計算するのがめんどくさい(平方完成するだけだが、、、)。
こういうのは数式処理ソフトに計算してもらおう、ということでSympyの出番。
まあそんなに難しくないのでリストだけ。ただ、r=一定のuの項だけがうまく簡略化できなかった。。。
うまいやり方知っている人は教えて!
陽性者数のカウントには厚生労働省のサイトを使っているのですが、9/26に全数把握止めたのを受けて9/27からフォーマットが変わった。
前はこんな感じでした。
それはそれとして、まずは各国のリニアスケール。日本がいうほど減ってない。韓国もか。他の国はあんまり増えてなさそうだがどれだけ管理してるのか(イギリスはほぼあきらめているのかな。。。)
次は各国のログログプロット。韓国がイギリスを抜いているが、これもどこまで信用できるやら。
日本のワクチン接種回数と陽性者数のログログプロット。それでもじわじわワクチン接種回数は増えてるかな。
日本・韓国・アメリカの人口に対する総陽性者プロット。日本ももう17%に達している。
最後は日本の詳細ログログプロット。何度も言っているが現状、全然減ってなくて過去の中では増加が激しいレベル。
最近のコメント