Last active
January 22, 2018 23:39
-
-
Save ti-nspire/cc77304e7a1a35708b8dcab65c4c337f to your computer and use it in GitHub Desktop.
Python による常微分方程式の数値解法 / 補外法
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# Python による常微分方程式の数値解法 / 補外法 | |
class Extrapo: | |
def __init__(self, funcs, t0, inits, h, numOfDiv=1): | |
self.funcs = funcs | |
self.t0 = t0 | |
self.inits = inits | |
self.h = h / numOfDiv | |
self.numOfDiv = numOfDiv | |
def __midPoints(self, funcs, t0, inits, H, S=7): | |
dim = len(funcs) | |
RANGE = 2**S+2 | |
listRange = [None for i in range(RANGE)] | |
temp = [None for i in range(dim)] | |
listS = [None for i in range(S)] | |
t = listRange.copy() | |
t[0] = t0 | |
x = [[] for i in range(dim)] | |
X = x.copy() | |
for i in range(dim): | |
x[i] = listRange.copy() | |
x[i][0] = inits[i] | |
X[i] = listS.copy() | |
for s in range(S+1): | |
N = 2**s | |
hs = H/N | |
t[1] = t[0] + hs | |
for i in range(dim): | |
x[i][1] = x[i][0] + hs * funcs[i](t[0], *inits) | |
for n in range(0, N): | |
t[n+2] = t[n] + 2 * hs | |
for i in range(dim): | |
temp[i] = x[i][n+1] | |
for i in range(dim): | |
x[i][n+2] = x[i][n] + 2 * hs * funcs[i](t[n+1], *temp) | |
for i in range(dim): | |
X[i][s-1] = (1/4) * (x[i][N-1] + 2 * x[i][N] + x[i][N+1]) | |
return X | |
def __neville(self, listOfLists): | |
dim = len(listOfLists) | |
S = len(listOfLists[0]) | |
output = [None for i in range(dim)] | |
x = [v for v in listOfLists] # x = [[xList],[yList],[List]] | |
for i in range(dim): | |
x[i].insert(0, None) # x = [[0,xList],[0,yList],[0,List]] # 添字をテキストと揃える手段として、使わない [0] 番目の要素を設けておく。 | |
table = [[[None for i in range(S+1)] for i in range(S+1)] for i in range(dim)] # 添字をテキストと揃える手段として、使わない [0] 行目、[0] 列目を設けておく。 | |
for i in range(1,S+1): # テーブルの [1] 列目に点列を並べる。 | |
for j in range(1,S+1): | |
for k in range(dim): | |
table[k][i][j] = x[k][i] | |
for i in range(2,S+1): | |
for j in range(2,i+1): | |
for k in range(dim): | |
table[k][i][j] = table[k][i][j-1]+(table[k][i][j-1]-table[k][i-1][j-1])/(4**(j-1)-1) | |
output[k] = table[k][S][S] | |
return output | |
def update(self): | |
for i in range(self.numOfDiv): | |
self.inits = self.__neville(self.__midPoints(self.funcs, self.t0, self.inits, self.h)) | |
self.t0 = self.t0 + self.h | |
return self | |
def print(self): | |
print(self.t0, self.inits) | |
return self | |
######## | |
# test # | |
######## | |
if __name__ == "__main__": | |
# 解くべき聯立微分方程式を定義する。リストで括っておく。 | |
def xDot(t, x, y): # x'(t) = y(t) | |
return y | |
def yDot(t, x, y): # y'(t) = t - x(t) | |
return t - x | |
funcs = [xDot, yDot] | |
# 独立変数の開始値と終了値とを指定する。 | |
t0 = 0 | |
tMax = 7.5 | |
# 従属変数の初期値を指定する。リストで括っておく。 | |
x0 = 0 | |
y0 = 0 | |
inits = [x0, y0] | |
# 刻み幅を指定する。2 の整数乗分の 1 にすることが望ましい。 | |
h = 1/2**2 | |
# 1 刻みの内部分割数を指定する場合は 2 の整数乗にすることが望ましい。 | |
# numOfDiv = 2**4 | |
# 1 刻みだけ計算する函数を実体化して、 | |
sol = Extrapo(funcs, t0, inits, h) | |
# 初期値を更新しながら必要な回数だけ実行を繰り返す。 | |
while sol.t0 < tMax: | |
sol.update().print() |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment