Excel VBAからOfficeスクリプト(TypeScript)へ乗換るための数値計算(2) 4段4次のルンゲクッタ法でローレンツ方程式を計算してみる。

このブログではかつて、Excel VBAでいろいろな数値計算をするのが看板だった。例えばこんなの。

Excelマクロ（VBA)で数値計算-複素TDGL方程式、蔵本シバシンスキー方程式、ピタゴラスの三体問題、シュレーディンガー方程式、FPU問題などなど

Dormand-Prince(ルンゲ・クッタ8次）の有名なルーチンDOP853をExcel VBAに移植

Excel VBAで複素数/FFTが使えるライブラリ

ただ、VBAからOfficeスクリプトにマイクロソフトは置き換えようとしているように見える。（Python in Excelは別用途）

で先日こんなのをやってみた。

Excelの自動化がいつの間にかVBAじゃなくてデスクトップ版でもOfficeスクリプト（TypeScript)で出来るようになっていた。Python in Excelとか触っているうちに出遅れた…そこで先日作った複素数クラスを使ってセルに入力した値を複素数計算できるようにした。 

VBAよりはるかに簡単に数値計算プログラム書けそう、ということで前回を(1)として今回の(2)ではもう何回やったかわからない4段4次のルンゲクッタ法でローレンツ方程式を計算してみよう。

プログラムはこんな感じで簡単。

function main(workbook: ExcelScript.Workbook) {

    let selectedSheet = workbook.getActiveWorksheet();

    let x : Array<number> = [0.1, 0.1, 0.1]; 

    let t: number = 0;

    const h: number = 0.01;

    for (let i = 0; i < 5000; i++) {

      selectedSheet.getCell(1 + i, 0).setValue(t);

      selectedSheet.getRangeByIndexes(1 + i, 1, 1, 3).setValues([x]);

      x = RungeKutta(t, x, h);

      t += h;

    }

    let chart = selectedSheet.addChart(ExcelScript.ChartType.xyscatterLinesNoMarkers, selectedSheet.getRange("B1:D1").getExtendedRange(ExcelScript.KeyboardDirection.down));

    chart.getTitle().setText("ローレンツ方程式");

}

function func(t : number, x : Array<number>) : Array<number> {

  let f : Array<number> = new Array(x.length);

  const p : number = 10.0;

  const r : number = 28.0;

  const b : number = 8.0/3.0;

  f[0] = -p * (x[0] - x[1]);

  f[1] = -x[0] * x[2] + r * x[0] - x[1];

  f[2] = x[0] * x[1] - b * x[2];

  return f;

}

function AddArray(a : Array<number>, b : Array<number>, c : number) : Array<number> { 

  let tmp : Array<number> = new Array(a.length);

  for (let i = 0; i < a.length; i++) {

    tmp[i] = a[i] + c * b[i];

  }

  return tmp;

}

function RungeKutta(t: number, x: Array<number>, h : number) : Array<number> {

  const k1: Array<number> = func(t, x);

  const k2: Array<number> = func(t + h / 2, AddArray(x, k1, h / 2));

  const k3: Array<number> = func(t + h / 2, AddArray(x, k2, h / 2));

  const k4: Array<number> = func(t + h, AddArray(x, k3, h));

  let x_next : Array<number> = new Array(x.length);

  for (let i = 0; i < x.length; i++) {

    x_next[i] = x[i] + h * (k1[i] + 2*k2[i] + 2*k3[i] + k4[i])/6;

  }

  return x_next;

}

実行するとすぐにこんなグラフが得られた！これは簡単。

コードエディタではこんな感じになって、補完とかしてくれる。

このシリーズもちょっとやっていこう（続く）。