FIFAワールドカップカタール2022で今使われているボールはAdidasとKinexon共同開発のUWB＋慣性センサ内蔵で位置が高精度にわかる。ICはQorvo(Decawave)のDW1000を使ってる。タグはアンテナもリファレンスそのまま、アンカーは外部PA追加みたい。

この記事を見た。

ワールドカップで審判支援「IoTボール」のスゴイ仕組み／世界のフットボールテック（１）

おお、すごいな。これはUWBの正しい使い方の気がする（AirTagでストーカー被害を増やすより、、、）。

重さとか大丈夫なのかなと素人考えだけどAdidasがその辺ちゃんと調整しているんでしょう。

ADIDAS REVEALS THE FIRST FIFA WORLD CUP™ OFFICIAL MATCH BALL FEATURING CONNECTED BALL TECHNOLOGY

こんな感じで入っている。

このボールをタグとしてグランドに複数置かれたアンカーで位置検知するということだそうだ。

https://kinexon.com/products/football/

で中身が気になるところ。FCCでKinexonを探すと、

まずタグ。ICはQorvo（Decawave）のDW1000だな。アンテナがDWM1001Cと同じだし、リファレンス通り作っている感じ。

https://www.qorvo.com/products/p/DW1000

次はアンカー。ICは同じだが、外部PAを使っているよう。SPDTでTx/Rxを切り替えている。

アプリケーションノートにそのまま乗っている。これもリファレンス通りっぽい。

どんどんスポーツもテクノロジー取り入れてるな。これはうかうかできない。

「人の体の1日の水分の出入り（代謝回転量）の予測式」をカシオの高精度計算サイトkeisan.casio.jpに自作式としてUP!サイエンスに論文が出たということで。

このまえこのニュースを見た。

人は一日に体内からどれくらい水分を失う？ 初の正確な計算式

おおそういうのがわかるのか。

筑波大学のプレスリリースはこちら：

https://www.tsukuba.ac.jp/journal/pdf/p20221125040000.pdf

サイエンスの論文：

Variation in human water turnover associated with environmental and lifestyle factors

https://www.science.org/doi/10.1126/science.abm8668

式そのものはものすごく簡単なので、ここはカシオの高精度計算サイトkeisan.casio.jpに作ってみよう。

リンクはこちら：

人の体の1日の水分の出入り（代謝回転量）の予測式 

画面はこんな感じ：

うちの祖母とかそうだったんですが、年取るとほとんど水分を取らない人がいるので

ちゃんとこれを見て水分補給してほしい。

Raspberry PI 3を数年ぶりに触る。Mathematicaをインストールすると最新版の13.1が入った！新機能で分数階微分とその微分方程式が解けるようになってるので例題通りやってみた。

Rasbberry PIを使いたいという人がいたので、私も昔使ってたよ、という話をしたが完全に何もかも忘れている。。。

なので思い出すためにちょっと触ってみる。最新の４は持ってないので3で。

まずはインストール。

あれ？昔と全然違ってImagerというのを使うようだ。

https://www.raspberrypi.com/software/

これは簡単！でもMathmaticaがプレインストールされてない、、、

ということでここからダウンロードしてインストール。

https://www.wolfram.com/raspberry-pi/index.php.ja?source=footer

おお、最新の13.1が入るんだな。これも簡単にできた。起動してみた。

13.1の最新の機能は、、、

https://www.wolfram.com/mathematica/quick-revision-history.html.ja?footer=lang

いろいろあるな。すぐ試せそうなのは

FractionalD，CaputoD，アップデートされたDSolveによる分数階微分および分数階微分方程式のサポート

かな。例題はDSolveのところにある。

https://reference.wolfram.com/language/ref/DSolve.html

やってみると、特に遅くもなくすぐ計算できた！

結構使えそうな感じ。

先延ばしの普遍的法則（課題提出、夏休みの宿題、進捗どうですか、などなど）N(t)=M/(D-t+C)-M/(D+C)を導出、カシオの高精度計算サイトkeisan.casio.jpにUP!締め切りギリギリじゃないとみんなやらない感覚にあってるな。

Fermat's Libraryで「"A universal law of procrastination（先延ばしの普遍的法則（万有法則？）」を見た。

This week's paper: "A universal law of procrastination" by T.Durakiewicz.

The author analyses how humans deliver under deadlines (e.g tax returns). Turns out this can be modeled by an hyperbolic function (better than an exponential!)https://t.co/dsZDOUl9nU pic.twitter.com/UdI606IKyq

— Fermat's Library (@fermatslibrary) November 1, 2022

元はこちら：

https://physicstoday.scitation.org/doi/10.1063/PT.3.3064

なるほど。確かに感覚と合ってる。Cの扱いがピンとこないのでちょっと計算。

提出物が最初が出てこなくて最後締め切りギリギリになってばっーと出てくる感じは∝-1/tのような形になるだろう。

ある時間tでの提出数をN(t)とする。締め切りt=Dになってようやく全提出N(D)=Mになる。もちろんN(0)=0。てことでこんな形になるはず。

\[ N(t) = \frac{a}{b - t }- \frac{a}{b} \]

ちょっと書き直す。

DからちょっとあとCたってこのN(t）はt=b=D+Cで発散するとする

\[ \lim_{t \to D+C} N(t) = \infty \]

なので

\[ N(t) = \frac{a}{D - t +C}- \frac{a}{D+C} \]

ここでN(D)=Mなので

\[ \frac{1}{C}- \frac{1}{D+C} = \frac{M}{a} \]

未知数が2つあるが簡単のために（or Durakiewiczさんが実際に受け取った申請データから) a=Mとする。すると

\[ \frac{D}{C(D+C)}=1 \]

\[ C^{2} + D C -D = 0 \]

\[ C=\frac{-D\pm\sqrt{D^{2}+4D}}{2}　\]

Cはもちろん正の値なのでちょっと書き直して

\[ C=\frac{1}{2}(\sqrt{D}\sqrt{D+4}-D)　\]

と Durakiewiczさんの結果が再現できた。

ということでこれをカシオの高精度計算サイトkeisan.casio.jpに作ってみた。こちら：

先延ばしの普遍的法則（提出、宿題、進捗など）

画面：

説明：

夏休みの宿題、提出する課題、申請書類、研究の進捗、などなど人間は締め切りギリギリまで先延ばしをします。
Tomasz Durakiewiczさんが実例からそれら全部に当てはまる理論式、A universal law of procrastinationを提案しました。Fermat's Libraryでも紹介されました。

計算結果：

Visual C# (C_sharp)の数値計算ライブラリ MathNET Numericsを使う(10) 数値積分としてガウス・クロンロッド積分公式と二重指数関数型積分公式を試す。

今回は数値積分やってみよう。

これは例題がちゃんと公式サイトに書かれているので簡単。

https://numerics.mathdotnet.com/Integration.html

数値積分の方法としては

DoubleExponentialTransformation
GaussKronrodRule
GaussLegendreRule
NewtonCotesTrapeziumRule
SimpsonRule

が選べる。今回はガウス・クロンロッド（オーダーいろいろ変える）と二重指数関数型積分公式をやってみよう。

積分は

∫4/(1+x^2)dx (積分範囲[0,1]）
と

∫1/√(1-x^2)dx (積分範囲[-1,1]）

とする。結果はこちら。

ガウス・クロンロッドの方は予想通りの精度だが、二重指数関数型積分公式が悪いな。

なんで？Githubのコード見てみよう。

今回のコードはこちら。めちゃくちゃ簡単。

テキストでもコードを書いておきます。

using System;
using MathNet.Numerics.Integration;

namespace MathNetIntegral
{
    class Program
    {
        static void Main(string[] args)
        {
            double ApproximatePi;
            double err;

            Console.WriteLine("ガウス・クロンロッドの積分公式で∫4/(1+x^2)dx (積分範囲[0,1]）の計算");
            for (int n = 2; n <= 20; n++)
            {
                ApproximatePi =  GaussKronrodRule.Integrate(x => 4.0 / (1.0 + x * x), 0.0, 1.0, out _, out _, order : n);
                err = ApproximatePi - Math.PI;
                Console.WriteLine("次数" + n.ToString() + ":  " + ApproximatePi.ToString() +"    誤差: " +  err.ToString());
            }
            Console.WriteLine();

            Console.WriteLine("二重指数型数値積分で∫1/√(1-x^2)dx (積分範囲[-1,1]）の計算");
            double integrate = DoubleExponentialTransformation.Integrate(x => 1.0 / Math.Sqrt(1.0 - x * x), -1.0, 1.0, 1.0e-16);
            Console.WriteLine(integrate.ToString());

        }
    }
}

 

過去のもの：

Visual C# (C_sharp)の数学ライブラリ Math.NET Numericsを使う(1) 複素行列を定義して一次方程式や逆行列、行列式などを計算する。

Visual C# (C_sharp)の数学ライブラリ Math.NET Numericsを使う(2) 補間を行う（Interpolate) リニア、３次スプライン、有理関数などいろいろ使える。

Visual C# (C_sharp)の数学ライブラリ Math.NET Numericsを使う(3) 高速フーリエ変換（FFT）を実行する。FourierOptionsにMatlabとNumerical Recipesがあるのが意外。

Visual C# (C_sharp)の数学ライブラリ Math.NET Numericsを使う(4) 多項式フィッティングをして、Array.ConvertAllで一括でフィッティングデータを得る。

Visual C# (C_sharp)の数学ライブラリ Math.NET Numericsを使う(5) 常微分方程式の数値解法、4段4次のルンゲクッタ法がRungeKutta.FourthOrderの一文でできる。ローレンツ方程式を例としてやってみる

Visual C# (C_sharp)の数値計算ライブラリ MathNET Numericsを使う(6)　OptimizationのNelder-Mead SimplexでRosenbrock関数（５パラメータ）を最小になる点を探す。

Visual C# (C_sharp)の数値計算ライブラリ MathNET Numericsを使う(7)　OptimizationのLevenberg-Marquardt法（LevenbergMarquardtMinimizer）で非線形最小二乗法（回帰）でNISTの例題Rat43を計算する。

Visual C# (C_sharp)の数値計算ライブラリ MathNET Numericsを使う(8)　特異値分解（SVD）、主成分分析（PCA)を計算してみる（ちょうど奥村先生が記事を出されてたので）

Visual C# (C_sharp)の数値計算ライブラリ MathNET Numericsを使う(9)　いろんな確率分布の乱数（メルセンヌツイスタがベース）をヒストグラムにして描く。とりあえず正規分布とガンマ分布で。

ハロウィンということで昔Scratchで作ったハロウィンのかぼちゃランタンの三体問題の8の字解を再掲。シンプレクティックオイラー法で計算してます。

Scratchのリンクはこちら：

Jack O' Lantern's Three body problem

動画にしたもの。クリックで始まります。シンプレクティックオイラー法は1次なんでちょっと8の字からずれていっているけれど。

シンプレクティックオイラー法についてはこちらを。

https://www.research.kobe-u.ac.jp/csi-viz/members/kageyama/lectures/H27_latter/Analytical_Mechanics/note_160121a.pdf

日本語プログラミング言語なでしこで数値計算（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]}」を表示
　　ここまで
ここまで。

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

さて1回目と2回目はこういうのをやってみた。

日本語プログラミング言語なでしこで数値計算（１）　まずは何はともあれ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に規格化）
# β：感染伝達率, γ：回復率

で、これは西浦博さんの「感染症疫学のためのデータ分析入門」を参考にした。

ではプログラムのリンクはこちら：

感染症流行のSIRモデルをルンゲクッタ法で計算

結果の例はこちら。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座標変換]まで線描画

ここまで

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

さて前回はルンゲクッタ法でローレンツ方程式を解いて図示してみた。

日本語プログラミング言語なでしこで数値計算（１）　まずは何はともあれ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最大値)/ 画面高 を戻すこと

ここまで

さて、次は何しようかな…

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

最近プログラム記事のネタ切れ感が半端ない、、、Visual C#を使い始めたがまあ仕事で使っている部分もあるのであまりかけない。

流行りのRUSTなどより、使っている人がとても少なそうな言語でいろいろ遊ぼうと思ったり。

そこで、日本語でプログラミングできる「なでしこ」の存在を思い出した。

https://nadesi.com/top/

おお、いつの間にかブラウザベースでプログラミングができるようになっているのか。これはちょうどいい。

しかも保存して公開できるし。誰もやってなさそうな数値計算ネタをやっていこう。

まずは何はともあれローレンツ方程式だ。

リンクはこちら：

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座標変換]まで線描画
ここまで