crossingorder

or crossing order

2021.8.13

数値データからn個のグラフを描く際に，曲線が滑らかになるようにデータを並び替えるMathematica関数

問題の説明

パラメタxに依存するn個の値 yₖ(x) を，m個のxの値 x=x₁, x₂, …, xₘ に対して計算して，xy面にプロットする場合を考える。下図は，そのようなデータ(n=3, m=16)を，(xₗ, yₖ(xₗ))と(xₗ₊₁, yₖ(xₗ₊₁))を線分で結んでプロットした例である。丸点が(xₗ, yₖ(xₗ))，青橙緑の順にn=1,2,3である。この図では，yₖ(xₗ)は値の昇順，yₖ(x)<yₖ₊₁(x) に並んでいる。

各xₗで，このようなデータyₖ(xₗ)を並べかえて，(xₗ, yₖ(xₗ))と(xₗ₊₁, yₖ(xₗ₊₁))を線分で結んだときに，次の図のように，折れ線が滑らかになるようにするには，どうしたらよいか。例えば，x=14.73付近で，青と橙の線は，上の図では交差せずに急カーブで反発しているが，下の図では，滑らかに交差している。

アルゴリズム

試行錯誤の結果，2次の補間多項式で補外を行うと，安定してそのような滑らかな結び方ができるようである。手順は次の通り：

  1. 最初の3点 (xₗ, yₖ(xₗ)), l=1,2,3 は，滑らかにつながる順に並んでいるとする。

  2. fₗₖ(x)を3点，(xₛ, yₖ(xₛ)), s=l-2, l-1, lを通るxの2次式とする。これにx=xₗ₊₁を代入した値とyₜ(xₗ₊₁)の差の二乗，(yₜ(xₗ₊₁)-fₗₖ(xₗ₊₁))², t=k, k+1, k+2, …, n を最小にするtを探索し，yₜ(xₗ₊₁)を新たにyₖ(xₗ₊₁)とする。すなわち，yₖ(xₗ₊₁)とyₜ(xₗ₊₁)の値を交換する。これを，t=1, 2, …, nの順に繰り返す。

  3. l=3, 4, …, n-1の順に2を行う。

ただし，データ間隔 xₗ₊₁-xₗ があまり広いと，補間関数の値が暴れるので，線の結び方が変な場合は，データを適宜追加する必要がある。

プログラムリスト

以下は，これをMathematicaの関数で書いたものである。

crossingorder[o_] := Module[

{n = Length[o], ni = Length[o[[1]]], ja, jn, p, f, nf = 2},

ja = Table[Table[j, {j, n}], ni];

jn = Table[j, {j, n}];

Print["data=", ni, ", function=", n, "];

Do[

Do[

f = InterpolatingPolynomial[

Table[o[[ja[[i + k, j]], i + k]], {k, -nf, 0}], o[[1, i + 1, 1]]];

p = OrderingBy[

Table[o[[jn[[k]], i + 1]], {k, j, n}], (# - f)^2 &, 1][[1]] + j - 1;

{jn[[j]], jn[[p]]} = {jn[[p]], jn[[j]]};

, {j, n}

];

ja[[i + 1]] = jn;

, {i, nf + 1, ni - 1}

];

Table[o[[ja[[i, j]], i]], {j, n}, {i, ni}]

]

引数oは，データ(xₗ, yₖ(xₗ))を表すリスト

{ {{x₁, y₁(x₁)}, {x₂, y₁(x₂)}, …, {{xₘ, y₁(xₘ)}}, {{x₁, y₂(x₁)}, {x₂, y₂(x₂)}, …, {{xₘ, y₂(xₘ)}}, …, {{x₁, yₙ(x_1)}, {x_2, yₙ(x_2)}, …, {{xₘ, yₙ(xₘ)}} }

である。

サンプルコード

oms = {

{{14.62, 0.046773}, {14.64, 0.079457}, {14.66, 0.108559}, {14.68, 0.134723}, {14.7, 0.158429}, {14.72, 0.180047}, {14.74, 0.190402}, {14.76, 0.192134}, {14.78, 0.193817}, {14.8, 0.195453}, {14.85, 0.199346}, {14.9, 0.202971}, {14.95, 0.206341}, {15., 0.209464}, {15.05, 0.212351}, {15.1, 0.215013}},

{{14.62, 0.178733}, {14.64, 0.180867}, {14.66, 0.182911}, {14.68, 0.184878}, {14.7, 0.186778}, {14.72, 0.188618}, {14.74, 0.199868}, {14.76, 0.218126}, {14.78, 0.235012}, {14.8, 0.250685}, {14.85, 0.285374}, {14.9, 0.314832}, {14.95, 0.331745}, {15., 0.335138}, {15.05, 0.337924}, {15.1, 0.340201}},

{{14.62, 0.277765}, {14.64, 0.284036}, {14.66, 0.289629}, {14.68, 0.294649}, {14.7, 0.299179}, {14.72, 0.303285}, {14.74, 0.30702}, {14.76, 0.310427}, {14.78, 0.313544}, {14.8, 0.316402}, {14.85, 0.322577}, {14.9, 0.327615}, {14.95, 0.340157}, {15., 0.362137}, {15.05, 0.381362}, {15.1, 0.398282}}

};

ListLinePlot[oms, Mesh -> Full, MeshStyle -> PointSize[Medium]]

oma = crossingorder[oms];

ListLinePlot[oma, Mesh -> Full, MeshStyle -> PointSize[Medium]]