Python 学習:スプライン曲線

2014-12-06 :  PCクリニック
本文の前に、
-・・・ -・-
現時点での blogramのランクインカテゴリは、
6、1、0、1、1、 0、0、0、0、0(41)で、換算ポイント 74pt 。
3日前に予想したとおり、一昨日朝方「グルコサミン」が2位に、
「Firefox」も5位に、共にダウン。
そして、昨日朝方「グルコサミン」がトップ返り咲きで、今に至る。
相変わらず「C言語」「Firefox」「化学業界」「グルコサミン」のbg値が、
毎日の如く変動しているから。
・-・ - -・

さて、本文。

ふと、
”Python”で「スプライン曲線」を描くには?
と思い立った。

何はともあれ、「Python スプライン曲線」で検索してみた。

「西方見聞録」なるブログの、2010-08-19 の記事:
ROOT での Spline 補完 (python 記述)
がある。
  ROOT を使ったデータの Spline 補間についてのメモ。
  ここでは、3 次の spline 補間 TSpline3 を試している。
  ・・・・・
  3次の Spline補間(TSpline3)は、
  各データ点の間の区間に分けて多項式近似する「区分的多項式近似
  (piecewise-polynomial approximation)」の一種です。
  ・・・・・

“ROOT”って?
例の「Gohlkeさんのページ」では、
残念ながら、
  Other useful packages and applications
  not currently available on this page
に入っていた。


他の検索結果では、
「ソフトとハードの交差点」ブログの、March 08, 2013 記事:
【Python】時系列データをスプライン曲線で・・・【pandas】

“pandas”? ・・・
“<紙>さんLoG”内pandas記事


更に他の検索結果では、
「数値計算とかの備忘録(仮)」ブログの、2011-08-09 記事:
スプライン補間(3次)

ここには、投稿者の自作(?)コードが載っている。

取り敢えず、
早速、コピー&ペーストさせて頂き、
実行。

当然(?)ながら、'data.dat' ファイルが出来て、
これを“pyplot”で表示すると、・・・ブログにある図の様に出ました。

これか。


ソース・コード を眺めていたが、・・・

ブログの最初に載せてある、リンク先:
3 スプライン補間
この記事で、学習した。


結果、解説の通りにコード化。
<紙>流(所謂、我流コード)で作った。

対象データ点は、上記コードにある(X座標が、0.9 から 9.2 迄の14点)
を借用させて頂き、( xj, yj ) として、
先ずは、
補間関数の係数配列( aj, bj, cj, dj ) を計算し、
尚、ここでの変数(配列)は全てグローバルとした
その後、
( X座標を 0.005 刻みで、補間関数の値を求め、)グラフ描画する。

import numpy as np
import matplotlib.pyplot as plt

def spline(): # xj, yj ---> aj, bj, cj, dj
M=np.zeros([N+1,N+1]); v=np.zeros(N+1); u=np.zeros(N+1)
hj=np.zeros(N+1); dyj=np.zeros(N+1)
for j in range(N):
hj[j]=xj[j+1]-xj[j]; dyj[j]=yj[j+1]-yj[j]
for j in range(1,N):
M[j,j]=2*(hj[j-1]+hj[j]) # 左辺:マトリックス
M[j,j-1]=hj[j-1] #
M[j,j+1]=hj[j] #
v[j]=6*(dyj[j]/hj[j]-dyj[j-1]/hj[j-1]) # 右辺:ベクトル
u[1:N] = np.linalg.solve( M[1:N,1:N], v[1:N] ) # ======= 方程式求解
for j in range(N):
aj[j]=(u[j+1]-u[j]) / (6*hj[j])
bj[j]=u[j]/2
cj[j]=dyj[j] / hj[j] - hj[j] * (2*u[j]+u[j+1]) /6
dj[j]=yj[j]

def sval(x):
j=0
while x<xj[j] or x>xj[j+1]: j+=1
dx=x-xj[j]
return aj[j]*dx**3+bj[j]*dx**2+cj[j]*dx+dj[j]

# =================================================================
xj = [ 0.9, 1.3, 1.9, 2.1, 2.6, 3.0, 3.9, 4.4, 4.7, 5.0, 6.0, 7.0, 8.0, 9.2 ]
yj = [ 1.3, 1.5, 1.85, 2.1, 2.6, 2.7, 2.4, 2.15, 2.05, 2.1, 2.25, 2.3, 2.25, 1.95 ]

N = len(xj)-1
aj=np.zeros(N); bj=np.zeros(N); cj=np.zeros(N); dj=np.zeros(N)
spline() # ========== 3次スプライン曲線計算 ( xj, yj ---> aj, bj, cj, dj )

X=[]; Y=[]
for x in np.arange( xj[0], xj[-1], 0.005 ):
X.append(x); Y.append(sval(x))
plt.plot( xj, yj, 'ro' ) # 'o'(circle) 赤色:オリジナル
plt.plot( X, Y, 'b-', lw=2 ) # '-'(line) 青色:スプライン曲線
plt.show()



これで、自己満足。メデタシ、めでたし。


尚、上記「3 スプライン補間」のページは、
「計算機応用」学科の教材の一部でした。

他の参考書では、
3次スプライン補間法」<<「数値解析入門I」
もあった。


本日はここまで。


見ていただいた序でとは厚かましい限りですが、
お帰りに投票して頂けるとなお嬉しいです。 ⇒ blogram投票ボタン


141030

関連記事
スポンサーサイト

コメントの投稿

管理者にだけ表示を許可する

おきてがみ/blogram
blogram投票ボタン



おきてがみ

最新記事
カレンダー
05 | 2017/06 | 07
- - - - 1 2 3
4 5 6 7 8 9 10
11 12 13 14 15 16 17
18 19 20 21 22 23 24
25 26 27 28 29 30 -
月別アーカイブ
カテゴリ
最新コメント
検索フォーム
リンク
プロフィール

<紙>

Author:<紙>
ようこそ。
「パソコンヲタクの雑記帳」
もろもろなことを綴っています。
パソコン ヲタクってねくら?
画像は kami でなく kani です。

カウンター(fc2、i2i) /Google Analytics


i2i(from 2010-08-24)
Total =
Today  =  
Yesterday=
アンチエイジング

Google Analytics
ブックマーク