私のpython練習のために、「三次元スプライン補間」を書いてみた。

どうも、yuukivelです。荷重試験も無事終わり、設計としては試験飛行待つばかりとなりました。
今年度、私が設計した機体は、時間部門の高速機の中ではかなりスパンが大きい機体です。
某SNS内でも、短スパンが良いのか長スパンが良いのか、問題提起をしてくれたようで、私としては何よりです。
今年度のうちの機体は、まさに「今を流れる小型TT機へのアンチテーゼ」なのです。

さて、今回もまた数学的（？）なお話。
前回の記事にて、octaveから離れた他の言語に移ると書いた。その言語とは「python」である。
GUIを実装するという観点から、「C++」か「python」に悩んでいた。
「C++」の方も「liboctave」というoctaveをそのままC++でコンパイルしました、微分代数方程式でも何でもござれという凄まじいライブラリを見つけて、狂喜乱舞しこちらに傾きかけたのだが、多くの方に話を伺ったところ「C++は暗黒」という言葉をいただいた。そしてoctaveを進めて下さったお師匠様(勝手に弟子入り)にもpythonの方が生産性が高いと言うことで、結局pythonにした。
とりあえず科学計算系はnumpy scipy scikit matplotlib といったライブラリをそろえれば実現できるし、そのうちoctaveで書いた設計プログラムをGUI化して移植するのも可能そう。

で、pythonに一通り慣れるために「三次元スプライン補間」のコードを書いてみた。
プログラムを始めるときには強烈な目標がないと、モチベーションが続かないので、何かしら課題を課してやることにしている。
VBAの時は渦法のプロペラ設計プログラム(結局未完)、octaveの時は主翼設計プログラムみたいな感じ。

三次元スプライン補間は昔記事に書いたように、離散データである空力解析データを如何に設計に活かすか。という問題で大きな役割を果たす。
スプライン補間が出来るか出来ないかで、設計の自由度は大きく変わるといっていい。
今年度の私の設計のような、激しく最適化を行う計算では、空力解析データのスプライン補間が必須となる。
図面出力もoctaveの豊富な補間関数に寄るところが大きい。

ということで、コードを書いた。それがこちら

#-------------------------------------------------------------------------------
# Name:        三次元スプライン補間
# Purpose:     私のpython練習
#
# Author:      Yuukivel
#
# Created:     05/12/2012
#-------------------------------------------------------------------------------
# -*- coding: utf-8 -*-

#スプライン基準点x,yの入力　補間点のxの値val
x=[1,1.5,2,2.5,3,3.5,4,4.5,5]
y=[1,2.25,4,6.25,9,12.25,16,20.25,25]
val=3

if len(x)!=len(y):
    print('x,yの長さが違う')
    exit()

#スプライン解行列Aの作成のための値計算
h=[x[1]-x[0]]
v=[0]
u=[0]
A=[(len(x)-2)*[0]]

j=1
while j<len(x)-1:
    if j!=len(x)-2:
        A.append((len(x)-2)*[0])
    h.append(x[j+1]-x[j])
    v.append(6*((y[j+1]-y[j])/h[j]-(y[j]-y[j-1])/h[j-1]))

    j+=1

j=0
while j<len(x)-2:
    #スプライン解行列Aの作成
    if j==0:
        A[0][0]=2*(h[0]+h[1])
        A[0][1]=h[1]

    elif j==(len(x)-3):
        A[j][j-1]=h[j-2]
        A[j][j]=2*(h[j-2]+h[j-1])
    else:
        A[j][j-1]=h[j]
        A[j][j]=2*(h[j]+h[j+1])
        A[j][j+1]=h[j+1]


    j+=1

#連立方程式を解く　scipyつかったのは逆行列演算めんどうくさいから
import numpy as np
import scipy as sp

Asp=sp.mat(A)
vsp=(sp.mat(v[1:len(x)-1])).T
usp=(Asp.I)*vsp
usp=usp.tolist()

for j in range(len(x)-1):
    if j==len(x)-2:
        u.append(0)
    else:
        u.append(usp[j][0])

#三次の近似関数の係数a,b,c,dを決定
a=[(u[1]-u[0])/(6*(x[1]-x[0]))]
b=[u[0]/2]
c=[(y[1]-y[0])/(x[1]-x[0])-1/6*(x[1]-x[0])*(2*u[0]+u[1])]
d=[y[0]]
for j in range(1,len(x)-1):
    a.append((u[j+1]-u[j])/(6*(x[j+1]-x[j])))
    b.append(u[j]/2)
    c.append((y[j+1]-y[j])/(x[j+1]-x[j])-1/6*(x[j+1]-x[j])*(2*u[j]+u[j+1]))
    d.append(y[j])

#補間点の位置を特定
for plc in range(len(x)-1):
    if val<x[0]:
        break
    elif val>=x[plc]:
        if val<x[plc+1]:
            break

#解の提示
ans=a[plc]*(val-x[plc])**3+b[plc]*(val-x[plc])**2+c[plc]*(val-x[plc])+d[plc]
print(ans)

引用、改変等は自由に行って下さい。再配布の時は一声戴けると嬉しいです。

まぁ読んでいただければ分かるとおり、numpy scipyの使用は逆行列を求めるところ以外避けている。
だって、逆行列を求める掃き出し法めんどうくさいんだもん。
だって、scipyの中にspline補間入ってるんだもん。

pythonに慣れることが目的だったので、まあこんな出来かなぁと。リストの参照先の関係で、代入したところ以外の値も変わるのにはびっくりした。
正直、本気でこの内容書くなら最初からnumpy scipy使うべき。

そして思ったのが、昔あんなに苦労して挫折したスプライン補間も、大学の2年にもなってくるとすらすら理解できるもんだなぁとしみじみ。
やっぱり、大学ってすごいわ。