微分方程式の数値解の計算は4次R-K法がよく用いられる。精度が高いらしい。そのアルゴリズムは以前にお勉強したので知っているしgFortran, C++, Pythonでも試した。数式通りに書くだけである。。。
それをF#で書くとどうなるか。F#の特徴は変数はimmutableであるとか回帰を使うのだとか能書きが多いのである。回帰を使って どう書くのよぉぉぉぉ・・・ よーし チャッターズに訊くである。
chatGPT3.5先生の回答である。
// 微分方程式 dy/dt = f(t, y) を 4 次 Runge-Kutta 法で解く
let rungeKutta4 f t0 y0 h n =
let rec loop t y i acc =
if i > n then List.rev acc // 計算終了時、解のリストを返す
else
// 4 次 Runge-Kutta 法の計算
let k1 = h * f t y
let k2 = h * f (t + h / 2.0) (y + k1 / 2.0)
let k3 = h * f (t + h / 2.0) (y + k2 / 2.0)
let k4 = h * f (t + h) (y + k3)
let yNext = y + (k1 + 2.0 * k2 + 2.0 * k3 + k4) / 6.0
loop (t + h) yNext (i + 1) ((t, yNext) :: acc)
loop t0 y0 0 [(t0, y0)]
// 微分方程式 dy/dt = -2y (指数関数的減衰) を解く例
let f t y = -2.0 * y
// 初期条件
let t0 = 0.0 // 初期時刻
let y0 = 1.0 // 初期値 y(0) = 1
let h = 0.1 // 時間刻み幅
let n = 20 // ステップ数
// 数値解を計算
let solution = rungeKutta4 f t0 y0 h n
// 結果を出力
solution
|> List.iter (fun (t, y) -> printfn "t = %.2f, y = %.6f" t y)
ま～なんてことでしょ。回帰ではなくて再帰でしたわ。再帰関数を多用するのがF#らしい書き方らしいのですわ。回帰というのは機械学習の回帰でしたわ。
しかし、このCodeは・・・・・C++で書いても似たようなものだがな・・・そりゃ同じアルゴリズムだからなあ(´・ω・`)
ここはあれだな、つまり、それだよ。再帰関数をいろいろ書いて慣れないことにはどこならんほーるどですわな。ま～その前に上の結果が正しいのか出力されたデータをグラフ化するなり、解は分かっているので誤差がどのぐらいかを確認したりしてからだな。。。。もちろんF#でグラフ出したりするのだよよん・・・・
グラフ描画は昨日、一つは動いたが残り2つも動かしたい次第である。ブラウザ表示のものではなくて。Python + matplotlibのように簡単に使えるようになると良きかな良きかな。
蛇足 簡単にグラフ出す方法 出力結果をコピペして秀丸でdata.txtとする。VScodeのTerminalで > wgnuplot と叩いて GnuPlotを起動して data.txtを食わせる。どの列の数字かは列を指定する。それでグラフは出る。ま～なんて簡潔な（笑）
それをDeedleとかXPlotとかをインストールして使い方を調べて、それからグラフ化するってんですって・・・・・・・・・・・・・・(´・ω・`)
posted by toinohni at 07:54| 東京 ☀|
Comment(0)
|
ソフト系雑学
|
|
