2014年11月11日火曜日

TSplineの境界条件の設定の方法

ROOTで点を補間しようと思った時に、点のある区間はスプライン補間で行い、その外は別の関数で補間することを考えた。
その際に、端点でスプラインと任意の関数が滑らかにつながるようにしたい。
このとき、TSplineに境界条件を入れることを考える。
CernのTSpline3の説明を見ると、

TSpline3(const char* title, const TGraph* g, const char* opt = 0, Double_t valbeg = 0, Double_t valend = 0)

となっており、valbegとvalendに端点の1階微分を境界条件として入れてあげれば良さそう。
デフォルトでは0となっているoptを設定しなければならない。
このoptに何を設定するかがわかりづらかった。
どうやらこのページを見ると、optにcharで"b1e1"と入れると、1階微分の境界条件が設定されるようだ。

試しにやってみたのが以下のプログラム(pythonで動かしてみた)。
ここでは、外の点は端から2点をつないだ直線でつなぐことを考えて、境界条件として端から2点の傾きを求めた。

>>> import ROOT as pR
>>> import array
>>> a = array.array("d",[1,2,3,4,5]) 
>>> b = array.array("d",[10,3,7,4,8])
>>> g = pR.TGraph(len(a),a,b)
>>> g.SetMarkerStyle(2)
>>> s = pR.TSpline3("s",g) # 境界条件を設定しないスプライン
>>> s.SetLineColor(2) #赤色
>>> tmp = [pR.Double(),pR.Double(),pR.Double(),pR.Double()]
>>> # 左端の傾き
>>> g.GetPoint(0,tmp[0],tmp[1]) # 始点
>>> g.GetPoint(1,tmp[2],tmp[3]) # 始点のお隣
>>> derF = (tmp[3]-tmp[1])/(tmp[2]-tmp[0])
>>> # 右端の傾き
>>> g.GetPoint(g.GetN()-2,tmp[0],tmp[1]) # 終点一つ前
>>> g.GetPoint(g.GetN()-1,tmp[2],tmp[3]) # 終点
>>> derE = (tmp[3]-tmp[1])/(tmp[2]-tmp[0])
>>> s2 = pR.TSpline3("s2",g,"b1e1",derF,derE) # 境界条件を設定したスプライン
>>> s2.SetLineColor(3) # 緑色
>>> 
>>> g.Draw("AP")
>>> s.Draw("same")
>>> s2.Draw("same")