def gevrey_tanh(T, n, sigma=sigma_tanh, K=K_tanh):
"""
Provide the flat output y(t) = phi(t), with the gevrey-order
1+1/sigma, and the derivatives up to order n.
:param t: [0, ... , t_end] (numpy array)
:param n: (integer)
:param sigma: (float)
:param K: (float)
:return: np.array([[phi], ... ,[phi^(n)]])
"""
t_init = t = np.linspace(0., T, int(0.5*10**(2+np.log10(T))))
# pop
t = np.delete(t, 0, 0)
t = np.delete(t, -1, 0)
# main
tau = t/T
a = dict()
a[0] = K*(4*tau*(1-tau))**(1-sigma)/(2*(sigma-1))
a[1] = (2*tau - 1)*(sigma-1)/(tau*(1-tau))*a[0]
for k in range(2, n+2):
a[k] = (tau*(1-tau))**-1 * ((sigma-2+k)*(2*tau-1)*a[k-1]+(k-1)*(2*sigma-4+k)*a[k-2])
yy = dict()
yy[0] = np.tanh(a[1])
if n > 0:
yy[1] = a[2]*(1-yy[0]**2)
z = dict()
z[0] = (1-yy[0]**2)
for i in range(2, n+1):
sum_yy = np.zeros(len(t))
for k in range(i):
if k == 0:
sum_z = np.zeros(len(t))
for j in range(i):
sum_z += -sm.comb(i-1, j)*yy[j]*yy[i-1-j]
z[i-1] = sum_z
sum_yy += sm.comb(i-1, k)*a[k+2]*z[i-1-k]
yy[i] = sum_yy
# push
phi = np.nan*np.zeros((n+1, len(t)+2))
for i in range(n+1):
phi_temp = 0.5*yy[i]
if i == 0:
phi_temp += 0.5
phi_temp = np.insert(phi_temp, 0, [0.], axis=0)
phi[i, :] = np.append(phi_temp, [1.])
else:
phi_temp = np.insert(phi_temp, 0, [0.], axis=0)
# attention divide by T^i
phi[i, :] = np.append(phi_temp, [0.])/T**i
return phi, t_init
评论列表
文章目录
