前々回はピタゴラスの三体問題という全員がどこかに飛んで行ってしまうものをやってみた。
今回は周期解。実はもう数百万とか解が見つかっているそうだが、まずは一番簡単で有名な8の字解。
前々回と初期値と質量が違うだけなのでコードは省略。こちらに初期値と質量の値があります。
Python+Scipyでルンゲクッタ8次のDOP853(Dormand Prince)を使う（その7） 三体問題の周期解いろいろ（１） 8の字解
そして結果の動画。
お互い追いかけあっているようで面白い。
次は別の周期解でやってみる。
今回は千葉逸人さんが中二病用語と言っていたこれをやってみる。
詳しくはスカラーペディアに。
この式を計算してます。
早速動画から。
ぐるぐる回っていたと思ったら突然広がって戻るのが面白い。スカラーペディアに出ていた図ともよく合ってる。
ソースコードはこちら。
using UnityEngine;
using System;
using System.Collections.Generic;
using Microsoft.Research.Oslo;
using System.Linq;
public class MathNET05 : MonoBehaviour
{
List<float> x0 = new List<float>();
List<float> x1 = new List<float>();
List<float> x2 = new List<float>();
int count;
int n;
void Start()
{
double t0 = 0.0;
double tmax = 1000.0;
n = 20000;
double dt = tmax / n;
var sol = Ode.RK547M(
0,
new Vector(0.012277918, -2.356078578, 0.018241293),
(t, x) => Bluesky(t, x),
new Options
{
AbsoluteTolerance = 1e-12,
RelativeTolerance = 1e-12
});
var points = sol.SolveFromToStep(t0, tmax, dt).ToArray();
n = points.Length;
count = 0;
foreach (var s in points)
{
x0.Add(Convert.ToSingle(s.X[0]));
x1.Add(Convert.ToSingle(s.X[1]));
x2.Add(Convert.ToSingle(s.X[2]));
}
}
void Update()
{
if (count < n)
{
transform.position = new Vector3(x2[count], x1[count], x0[count]);
}
count += 5;
}
Vector Bluesky(double t, Vector x)
{
double[] xdot = new double[3];
double myu = 0.456;
double eps = 0.0357;
xdot[0] = x[0] * (2.0 + myu - 10.0 * (x[0] * x[0] + x[1] * x[1]))
+ x[2] * x[2] + x[1] * x[1] + 2.0 * x[1];
xdot[1] = -x[2] * x[2] * x[2] - (1 + x[1]) * (x[2] * x[2] + x[1] * x[1] + 2.0 * x[1])
- 4.0 * x[0] + myu * x[1];
xdot[2] = (1 + x[1]) * (x[2] * x[2]) + x[0] * x[0] - eps;
return new Vector(xdot);
}
}
|
Google ColabでJuliaが使えるようになっていたということで普及しそう。しかし私は全然いじったことがない…
この機会にちょっと触ってみよう。とりあえず公式サイトのマニュアルをペラペラと読んで5%くらいは分かった。
任意精度演算がデフォルトで可能なのを知ったので、最近うちのブログによくアクセスがあるタッパーの自己言及式をやってみる。big(数）としただけで使えるのは便利。ただ途中で改行するとエラーになったので（括弧の中なのに？）一行が長い…
コード：
using Plots
setprecision(BigFloat, 2000)
k = big(960939379918958884971672962127852754715004339660129306651505519271702802395266424689642842174350718121267153782770623355993237280874144307891325963941337723487857735749823926629715517173716995165232890538221612403238855866184013235585136048828693337902491454229288667081096184496091705183454067827731551705405381627380967602565625016981482083418783163849115590225610003652351370343874461848378737238198224849863465033159410054974700593138339226497249461751545728366702369745461014655997933798537483143786841806593422227898388722980000748404719)
nx = 106
ny = 17
z = zeros(Int8, ny, nx + 1)
for i in 1:nx+1
for j in 1:ny
x = i - 1
y = k + j - 1
f = floor(floor(y / 17.0) * 2^(-17.0 * floor(x) - (floor(y) % 17)) % 2)
if f > 1/2
z[ny + 1 - j, nx + 2 - i] = 0
else
z[ny + 1 - j, nx + 2 - i] = 1
end
end
end
heatmap(z, size=(1300,200))
|
結果：
一瞬でできるな。
(過去のもの)
タッパーの自己言及式をPari/GPとExcelで描いてみる。
さて今回は非線形最小二乗法。
こちらを移植したもの。
using UnityEngine;
using MathNet.Numerics;
using MathNet.Numerics.LinearAlgebra;
using MathNet.Numerics.Optimization;
using System;
public class MathNETController : MonoBehaviour
{
public GameObject PointsPrefab;
GameObject line;
Vector<double> x;
Vector<double> y;
float scale;
int n;
int count;
// Start is called once before the first execution of Update after the MonoBehaviour is created
void Start()
{
var obj = ObjectiveFunction.NonlinearModel(OptimizeFunction, Gradient, Xdata, Ydata);
var solver = new LevenbergMarquardtMinimizer();
var result = solver.FindMinimum(obj, InitialValue1);
x = Vector<double>.Build.DenseOfArray(Generate.LinearSpaced(1000, 1.0, 15.0));
y = OptimizeFunction(result.MinimizingPoint, x);
line = GameObject.Find("Line");
scale = 50.0f;
n = x.Count;
count = 0;
for (int i = 0; i < Xdata.Count; i++)
{
GameObject point = Instantiate(PointsPrefab);
point.transform.position = new Vector3((float)Xdata[i], Convert.ToSingle(Ydata[i]/scale), 0f);
}
}
// Update is called once per frame
void Update()
{
if (count < n) {
line.transform.position = new Vector3((float)x[count], Convert.ToSingle(y[count] / scale), 0f);
count++;
}
}
Vector<double> Xdata = Vector<double>.Build.DenseOfArray(new double[] {
1.00, 2.00, 3.00, 4.00, 5.00, 6.00, 7.00, 8.00, 9.00, 10.00,
11.00, 12.00, 13.00, 14.00, 15.00
});
Vector<double> Ydata = Vector<double>.Build.DenseOfArray(new double[] {
16.08, 33.83, 65.80, 97.20, 191.55, 326.20, 386.87, 520.53, 590.03, 651.92,
724.93, 699.56, 689.96, 637.56, 717.41
});
Vector<double> InitialValue1 = Vector<double>.Build.DenseOfArray(new double[] { 100, 10, 1, 1 });
Vector<double> OptimizeFunction(Vector<double> p, Vector<double> x)
{
var y = Vector<double>.Build.Dense(x.Count);
for (int i = 0; i < x.Count; i++)
{
y[i] = p[0] / Math.Pow(1.0 + Math.Exp(p[1] - p[2] * x[i]), 1.0 / p[3]);
}
return y;
}
private Matrix<double> Gradient(Vector<double> p, Vector<double> x)
{
var prime = Matrix<double>.Build.Dense(x.Count, p.Count);
for (int i = 0; i < x.Count; i++)
{
prime[i, 0] = 1.0 / Math.Pow(1.0 + Math.Exp(p[1] - p[2] * x[i]), 1.0 / p[3]);
prime[i, 1] = -p[0] * Math.Exp(p[1] - p[2] * x[i]) /
(p[3] * Math.Pow(1.0 + Math.Exp(p[1] - p[2] * x[i]), 1.0 / p[3] + 1.0));
prime[i, 2] = p[0] * x[i] * Math.Exp(p[1] - p[2] * x[i]) /
(p[3] * Math.Pow(1.0 + Math.Exp(p[1] - p[2] * x[i]), 1.0 / p[3] + 1.0));
prime[i, 3] = Math.Log(1.0 + Math.Exp(p[1] - p[2] * x[i])) * p[0] /
(p[3] * p[3] * Math.Pow(1.0 + Math.Exp(p[1] - p[2] * x[i]), 1.0 / p[3]));
}
return prime;
}
}
|
結果の動画。計算結果を玉が軌跡を描きながら動く…
今回はWordで試してみる。
立ち上げたときからCopilotで下書きが前提になっていた。
猫のイラストを頼むとDesignerで描いてくれた。
そしてコラムも。
文章自体は挿入できるが、レイアウトが乱れているのは直せない（直接文書は編集できない）と言われた。
それが一番ストレスたまるところなんだが残念。
以前こういうのをやってみた。
この時は自分のノートPCのCPUしか使えなかったのでCNNがとても遅い…Google Colabでは無償でもGPU、TPUが使えるので比較してみる。
ランタイムのタイプを変更で選べる。無償だと2つしか選べないが…
まずはおまじない。
!pip install scikit-rf
from google.colab import drive
drive.mount('/content/drive')
path = "/content/drive/MyDrive/Colab Notebooks/"
|
次に設定とモデルの構築。今回はjaxを使う。データは以前作ったもの（上のリンク参照）。
import numpy as np
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
import os
os.environ["KERAS_BACKEND"] = "jax"
import keras
data_label = np.load(path + "data_label.npz")
data = data_label["data"].reshape(-1,200,5,1)
label = data_label["label"]
x_train, x_test, y_train, y_test = train_test_split(data, label, test_size=0.3, random_state=0)
# Functional APIでCNNを設定
inputs = keras.Input(shape=(200, 5, 1))
x = keras.layers.Conv2D(64, kernel_size=(10, 2), activation="relu")(inputs)
x = keras.layers.Conv2D(64, kernel_size=(10, 4), activation="relu")(x)
x = keras.layers.Flatten()(x)
outputs = keras.layers.Dense(6)(x)
# モデルの設定
model = keras.Model(inputs=inputs, outputs=outputs)
model.compile(loss = 'mean_squared_error' ,optimizer=keras.optimizers.Adam())
model.summary()
|
ここで学習を10epochやってみると…
batch_size = 32
epochs = 300
keras.utils.set_random_seed(1)
history = model.fit(
x_train,
y_train,
batch_size=batch_size,
verbose=0,
epochs=epochs,
validation_split=0.15,
)
y_pred = model.predict(x_test)
metric = keras.metrics.R2Score()
metric.update_state(y_test, y_pred)
result = metric.result()
print(result)
error = np.abs((y_test - y_pred)/y_test*100)
print(error.mean(axis=0))
|
結果の比較。
| 種類 | かかった時間 |
| CPU | 6分38 秒 |
| GPU | 16 秒 |
| TPU | 25 秒 |
GPUが爆速だった！これは自分のPCでやるよりずっと速い。こういう検討する時はGoogle Colab使う方向で行こう。
このポスト見て、
ｴｯﾎｴｯﾎ
— うお座 (@zlatanera17) February 24, 2025
ママに夜ご飯いらないって伝えないと
ｴｯﾎｴｯﾎ pic.twitter.com/8dshR7nJgR
さっそくいらすとやさんがイラストにしていたので、
走るフクロウのヒナのイラストhttps://t.co/wXg8Y1ue9I
— いらすとや (@irasutoya) February 28, 2025
さっそくUnityを使って動画にしてみた。
動きがいまいちか…
今回はＳパラメータの表示。dB表示は簡単だが、スミスチャートはどうする？と思ったらもうScottPlotに用意されていた。
https://scottplot.net/cookbook/5.0/SmithChart/
これなら簡単だ。ソースコードはこんな感じで（Touchstone読むところは面倒なので手抜き…）。
|
namespace PlotSmith
{
public partial class Form1 : Form
{
public Form1()
{
InitializeComponent();
radioButton1.Checked = true;
var axis = formsPlot1.Plot.Axes;
axis.Bottom.Label.Text = "Frequency[MHz]";
axis.Bottom.Label.FontSize = 21;
axis.Bottom.TickLabelStyle.FontSize = 21;
axis.Left.Label.Text = "S parameter[dB]";
axis.Left.Label.FontSize = 21;
axis.Left.TickLabelStyle.FontSize = 21;
}
private void button1_Click(object sender, EventArgs e)
{
List<double> freq = new List<double>();
List<List<double>> spara = new List<List<double>>();
StreamReader sr = new StreamReader("dea165550bt-2322a1-h.s2p");
while (!sr.EndOfStream)
{
string? line = sr.ReadLine();
if (line != null && line[0] != '!' && line[0] != '#')
{
string[] data = line.Split(' ');
freq.Add(double.Parse(data[0]));
List<double> list = new List<double>();
list.Add(double.Parse(data[1]));
list.Add(double.Parse(data[2]));
list.Add(double.Parse(data[3]));
list.Add(double.Parse(data[4]));
spara.Add(list);
}
}
int n = spara.Count;
if (radioButton1.Checked)
{
formsPlot1.Plot.Clear();
var fr = new double[n];
var s11 = new double[n];
var s21 = new double[n];
for (int i = 0; i < n; i++)
{
fr[i] = freq[i];
s11[i] = spara[i][0];
s21[i] = spara[i][2];
}
var s11plot = formsPlot1.Plot.Add.ScatterLine(fr, s11);
s11plot.LegendText = "S11";
s11plot.LineWidth = 5;
var s21plot = formsPlot1.Plot.Add.ScatterLine(fr, s21);
s21plot.LegendText = "S21";
s21plot.LineWidth = 5;
formsPlot1.Plot.Legend.FontSize = 21;
formsPlot1.Plot.ShowLegend();
formsPlot1.Plot.Axes.AutoScale();
formsPlot1.Refresh();
}
if (radioButton2.Checked)
{
formsPlot1.Plot.Clear();
var smith = formsPlot1.Plot.Add.SmithChartAxis();
var s11re = new double[n];
var s11im = new double[n];
for (int i = 0;i < n; i++)
{
var db = spara[i][0];
var phase = spara[i][1];
var mag = Math.Pow(10.0, db / 20.0);
s11re[i] = mag * Math.Cos(phase * Math.PI / 180.0);
s11im[i] = mag * Math.Sin(phase * Math.PI / 180.0);
}
var sline = formsPlot1.Plot.Add.ScatterLine(s11re, s11im);
sline.LineWidth = 5;
formsPlot1.Plot.Axes.AutoScale();
formsPlot1.Refresh();
}
}
}
}
|
実行結果の動画。dBとスミスチャートを切り替えられるようにした。
UnityでMath.NET Numericsという数値計算ライブラリを使うシリーズをやっているが、ルンゲクッタ法が4段4次のものしかないのでピタゴラスの三体問題のような精度が問題になるようなものはできず、Pythonで計算したものをCSVファイルにして読み込むということしかできなかった。
しかし、MicrosoftがOpen Solving Library for ODEs(OSLO)を出していることを知ったのでこれならできるな、とやってみた。
AssetsフォルダにMicrosoft.Research.Oslo.dllを置くと自動的に参照してくれた。
ソースコードは以下の通り。
using UnityEngine;
using System;
using System.Collections.Generic;
using Microsoft.Research.Oslo;
using System.Linq;
public class CsvReader : MonoBehaviour
{
List<float> x1 = new List<float>();
List<float> x2 = new List<float>();
List<float> x3 = new List<float>();
List<float> y1 = new List<float>();
List<float> y2 = new List<float>();
List<float> y3 = new List<float>();
int n;
int count;
GameObject Ball1;
GameObject Ball2;
GameObject Ball3;
// Start is called once before the first execution of Update after the MonoBehaviour is created
void Start()
{
Ball1 = GameObject.Find("Ball1");
Ball2 = GameObject.Find("Ball2");
Ball3 = GameObject.Find("Ball3");
var sol = Ode.RK547M(
0,
new Vector(1, 0, 3, 0, -2, 0, -1, 0, 1, 0, -1, 0),
(t, x) => threebody(t, x),
new Options
{
AbsoluteTolerance = 1e-12,
RelativeTolerance = 1e-12
});
var points = sol.SolveFromToStep(0, 100, 0.01).ToArray();
n = points.Length;
count = 0;
foreach (var s in points)
{
x1.Add(Convert.ToSingle(s.X[0]));
y1.Add(Convert.ToSingle(s.X[2]));
x2.Add(Convert.ToSingle(s.X[4]));
y2.Add(Convert.ToSingle(s.X[6]));
x3.Add(Convert.ToSingle(s.X[8]));
y3.Add(Convert.ToSingle(s.X[10]));
}
}
Vector threebody(double t, Vector x)
{
var f = new double[12];
var m1 = 3.0;
var m2 = 4.0;
var m3 = 5.0;
var qx1 = x[0];
var vx1 = x[1];
var qy1 = x[2];
var vy1 = x[3];
var qx2 = x[4];
var vx2 = x[5];
var qy2 = x[6];
var vy2 = x[7];
var qx3 = x[8];
var vx3 = x[9];
var qy3 = x[10];
var vy3 = x[11];
var Ra = Math.Pow(Math.Sqrt(Math.Pow((qx2 - qx1), 2.0) + Math.Pow((qy2 - qy1), 2)), 3);
var Rb = Math.Pow(Math.Sqrt(Math.Pow((qx3 - qx1), 2.0) + Math.Pow((qy3 - qy1), 2)), 3);
f[0] = vx1;
f[1] = m2 * (qx2 - qx1) / Ra + m3 * (qx3 - qx1) / Rb;
f[2] = vy1;
f[3] = m2 * (qy2 - qy1) / Ra + m3 * (qy3 - qy1) / Rb;
Ra = Math.Pow(Math.Sqrt(Math.Pow((qx1 - qx2), 2.0) + Math.Pow((qy1 - qy2), 2)), 3);
Rb = Math.Pow(Math.Sqrt(Math.Pow((qx3 - qx2), 2.0) + Math.Pow((qy3 - qy2), 2)), 3);
f[4] = vx2;
f[5] = m1 * (qx1 - qx2) / Ra + m3 * (qx3 - qx2) / Rb;
f[6] = vy2;
f[7] = m1 * (qy1 - qy2) / Ra + m3 * (qy3 - qy2) / Rb;
Ra = Math.Pow(Math.Sqrt(Math.Pow((qx1 - qx3), 2.0) + Math.Pow((qy1 - qy3), 2)), 3);
Rb = Math.Pow(Math.Sqrt(Math.Pow((qx2 - qx3), 2.0) + Math.Pow((qy2 - qy3), 2)), 3); ;
f[8] = vx3;
f[9] = m1 * (qx1 - qx3) / Ra + m2 * (qx2 - qx3) / Rb;
f[10] = vy3;
f[11] = m1 * (qy1 - qy3) / Ra + m2 * (qy2 - qy3) / Rb;
return f;
}
// Update is called once per frame
void Update()
{
if (count > n)
{
count = 0;
}
Ball1.transform.position = new Vector3(x1[count], y1[count], 0f);
Ball2.transform.position = new Vector3(x2[count], y2[count], 0f);
Ball3.transform.position = new Vector3(x3[count], y3[count], 0f);
count++;
}
}
|
結果をカメラを回しながら動画にしたもの。せっかくなんで背景に空もつけた。
ちゃんと計算してそうで一安心。
Visual C#で数値計算する時はMath.NET Numericsを使うのが便利だが、常微分方程式が4段4次のルンゲクッタ法しかないので硬い問題などがほぼ無理。そこでもうちょっといろいろあるのはないかなと思って探すとこれがあった。
適応型のルンゲクッタ、Dormand＆PrinceのRK547MとGear の後退差分式(BDF)が使える。まずはダウンロードして、、、
あれ？dllがない？自分でビルドしろということか。早速やってみると（古いバージョンなので）いろいろ怒られながらとりあえずできた。
使い方はとても簡単なので、早速やってみる。4段4次のルンゲクッタで計算するのがとても難しいピタゴラスの三体問題にしよう。
Python+Scipyでルンゲクッタ8次のDOP853(Dormand Prince)を使う（その5） ピタゴラスの三体問題を計算する。rtolとatolを設定しないと無茶苦茶になる。
Visual C#のコードはこんな感じで。RK547Mを使った。
using Microsoft.Research.Oslo;
namespace PlotRungeKutta02
{
public partial class Form1 : Form
{
public Form1()
{
InitializeComponent();
var axis = formsPlot1.Plot.Axes;
axis.Bottom.Label.Text = "x";
axis.Bottom.TickLabelStyle.FontSize = 21;
axis.Left.Label.Text = "y";
axis.Left.TickLabelStyle.FontSize = 21;
axis.Title.Label.Text = "3 body probrem";
axis.Title.Label.FontSize = 32;
axis.SetLimits(-6, 6, -6, 6);
formsPlot1.Refresh();
}
private void button1_Click(object sender, EventArgs e)
{
var sol = Ode.RK547M(
0,
new Vector(1, 0, 3, 0, -2, 0, -1, 0, 1, 0, -1, 0),
(t, x) => threebody(t, x),
new Options
{
AbsoluteTolerance = 1e-12,
RelativeTolerance = 1e-12
});
var points = sol.SolveFromToStep(0, 100, 0.01).ToArray();
var x1 = new double[points.Length];
var y1 = new double[points.Length];
var x2 = new double[points.Length];
var y2 = new double[points.Length];
var x3 = new double[points.Length];
var y3 = new double[points.Length];
int i = 0;
foreach (var s in points)
{
x1[i] = s.X[0];
y1[i] = s.X[2];
x2[i] = s.X[4];
y2[i] = s.X[6];
x3[i] = s.X[8];
y3[i] = s.X[10];
i++;
}
var p1 = formsPlot1.Plot.Add.ScatterLine(x1, y1);
p1.LineWidth = 5;
var p2 = formsPlot1.Plot.Add.ScatterLine(x2, y2);
p2.LineWidth = 5;
var p3 = formsPlot1.Plot.Add.ScatterLine(x3, y3);
p3.LineWidth = 5;
formsPlot1.Refresh();
}
Vector threebody(double t, Vector x)
{
var f = new double[12];
var m1 = 3.0;
var m2 = 4.0;
var m3 = 5.0;
var qx1 = x[0];
var vx1 = x[1];
var qy1 = x[2];
var vy1 = x[3];
var qx2 = x[4];
var vx2 = x[5];
var qy2 = x[6];
var vy2 = x[7];
var qx3 = x[8];
var vx3 = x[9];
var qy3 = x[10];
var vy3 = x[11];
var Ra = Math.Pow(Math.Sqrt(Math.Pow((qx2 - qx1), 2.0) + Math.Pow((qy2 - qy1), 2)), 3);
var Rb = Math.Pow(Math.Sqrt(Math.Pow((qx3 - qx1), 2.0) + Math.Pow((qy3 - qy1), 2)), 3);
f[0] = vx1;
f[1] = m2 * (qx2 - qx1) / Ra + m3 * (qx3 - qx1) / Rb;
f[2] = vy1;
f[3] = m2 * (qy2 - qy1) / Ra + m3 * (qy3 - qy1) / Rb;
Ra = Math.Pow(Math.Sqrt(Math.Pow((qx1 - qx2), 2.0) + Math.Pow((qy1 - qy2), 2)), 3);
Rb = Math.Pow(Math.Sqrt(Math.Pow((qx3 - qx2), 2.0) + Math.Pow((qy3 - qy2), 2)), 3);
f[4] = vx2;
f[5] = m1 * (qx1 - qx2) / Ra + m3 * (qx3 - qx2) / Rb;
f[6] = vy2;
f[7] = m1 * (qy1 - qy2) / Ra + m3 * (qy3 - qy2) / Rb;
Ra = Math.Pow(Math.Sqrt(Math.Pow((qx1 - qx3), 2.0) + Math.Pow((qy1 - qy3), 2)), 3);
Rb = Math.Pow(Math.Sqrt(Math.Pow((qx2 - qx3), 2.0) + Math.Pow((qy2 - qy3), 2)), 3); ;
f[8] = vx3;
f[9] = m1 * (qx1 - qx3) / Ra + m2 * (qx2 - qx3) / Rb;
f[10] = vy3;
f[11] = m1 * (qy1 - qy3) / Ra + m2 * (qy2 - qy3) / Rb;
return f;
}
}
}
|
結果はこちら。
ちゃんと計算できてるっぽい。Ode.RK547Mの代わりにOde.GearBDFにしても同様に計算できた。これは使えそう。
最近のコメント