This notebook uses the newer solver. This solver minimizes frequency domain error. hb_freq can also ignore the constant term ($\omega = 0$) in the solution process. The error in Fourier Series of the state-equation calculated state derivative as compared to that obtained by taking the derivative of the input state. Any variety of time points may be used to ensure substantial averaging over a single cycle.
In [100]:
%matplotlib inline
%load_ext autoreload
%autoreload 2
import scipy as sp
import numpy as np
import matplotlib.pyplot as plt
import mousai as ms
from scipy import pi, sin
In [101]:
# Test that all is working.
# f_tol adjusts accuracy. This is smaller than reasonable, but illustrative of usage.
t, x, e, amps, phases = ms.hb_freq(ms.duff_osc, np.array([[0,1,-1]]), omega = .7, f_tol = 1e-8)
print('Equation errors (should be zero): ', e)
print('Constant term of FFT of signal should be zero: ', ms.fftp.fft(x)[0,0])
In [102]:
# Using more harmonics.
t, x, e, amps, phases = ms.hb_freq(ms.duff_osc, x0 = np.array([[0,1,-1]]), omega = .1, num_harmonics= 1)
print('Equation errors (should be zero): ', e)
print('Constant term of FFT of signal should be zero: ', ms.fftp.fft(x)[0,0])
In [103]:
np.abs(e)
Out[103]:
Sometimes we can improve just by restarting from the prior end point. Sometimes, we just think it's improved.
In [104]:
t, x, e, amps, phases = ms.hb_freq(ms.duff_osc, x0 = x, omega = 0.1, num_harmonics= 7)
print('Errors: ', e)
print('Constant term of FFT of signal should be zero: ', ms.fftp.fft(x)[0,0])
In [105]:
# Let's get a smoother response
time, xc = ms.time_history(t,x)
plt.plot(time,xc.T,t,x.T,'*')
plt.grid(True)
print('The average for this problem is known to be zero, we got', sp.average(x))
In [106]:
def duff_osc2(x, v, params):
omega = params['omega']
t = params['cur_time']
return np.array([[-x-.01*x**3-.01*v+1*sin(omega*t)]])
In [107]:
t, x, e, amps, phases = ms.hb_freq(duff_osc2, np.array([[0,1,-1]]), omega = 0.8, num_harmonics=7)
print(amps, x, e)
print('Constant term of FFT of signal should be zero: ', ms.fftp.fft(x)[0,0])
time, xc = ms.time_history(t,x)
plt.plot(time, xc.T, t, x.T, '*')
plt.grid(True)
In [108]:
omega = np.linspace(0.1,3,200)+1/200
amp = np.zeros_like(omega)
x = np.array([[0,-1,1]])
for i, freq in enumerate(omega):
#print(i,freq,x)
try:
t, x, e, amps, phases = ms.hb_freq(duff_osc2, x, omega = freq, num_harmonics = 1)# , callback = resid)
#print(freq, amps, e)
amp[i]=amps[0]
except:
amp[i] = sp.nan
print(np.hstack((omega.reshape(-1,1), amp.reshape(-1,1))))
plt.plot(omega, amp)
Out[108]:
In [110]:
t, x, e, amps, phases = ms.hb_freq(duff_osc2, np.array([[0,1,-1]]), omega = 1.1, num_harmonics=1)
print(' amps = {}\n x = {}\n e = {}\n phases = {}'.format(amps, x, e, phases))
print('Constant term of FFT of signal should be zero: ', ms.fftp.fft(x)[0,0])
time, xc = ms.time_history(t,x)
plt.plot(time, xc.T, t, x.T, '*')
plt.grid(True)
In [111]:
phases
Out[111]:
In [112]:
omega = sp.linspace(0.1,3,90)+1/200
amp = sp.zeros_like(omega)
x = np.array([[0,-1,1]])
for i, freq in enumerate(omega):
#print(i,freq,x)
#print(sp.average(x))
x = x-sp.average(x)
try:
t, x, e, amps, phases = ms.hb_freq(duff_osc2, x, freq, num_harmonics=1)#, callback = resid)
amp[i]=amps[0]
except:
amp[i] = sp.nan
plt.plot(omega, amp)
Out[112]:
In [113]:
omegal = sp.arange(3,.03,-1/200)+1/200
ampl = sp.zeros_like(omegal)
x = np.array([[0,-1,1]])
for i, freq in enumerate(omegal):
# Here we try to obtain solutions, but if they don't work,
# we ignore them by inserting `np.nan` values.
x = x-sp.average(x)
try:
t, x, e, amps, phases = ms.hb_freq(duff_osc2, x, freq, num_harmonics=1, f_tol = 1e-6)#, callback = resid)
ampl[i]=amps[0]
except:
ampl[i] = sp.nan
plt.plot(omegal, ampl)
Out[113]:
In [114]:
plt.plot(omegal,ampl)
plt.plot(omega,amp)
#plt.axis([0,3, 0, 10.5])
Out[114]:
In [115]:
from scipy.optimize import newton_krylov
In [116]:
def duff_amp_resid(a):
return (mu**2+(sigma-3/8*alpha/omega_0*a**2)**2)*a**2-(k**2)/4/omega_0**2
In [117]:
mu = 0.05 # damping
k = 1 # excitation amplitude
sigma = -0.9 #detuning
omega_0 = 1 # driving frequency
alpha = 0.1 # cubic coefficient
In [118]:
newton_krylov(duff_amp_resid,-.1)
Out[118]:
In [119]:
sigmas = sp.linspace(-1,3,200)
amplitudes = sp.zeros_like(sigmas)
x = newton_krylov(duff_amp_resid,1)
for i, sigma in enumerate(sigmas):
try:
amplitudes[i] = newton_krylov(duff_amp_resid,x)
x = amplitudes[i]
except:
amplitudes[i] = newton_krylov(duff_amp_resid,0)
x = amplitudes[i]
plt.plot(sigmas,amplitudes)
Out[119]:
In [120]:
sigmas = sp.linspace(-1,3,200)
sigmasr = sigmas[::-1]
amplitudesr = sp.zeros_like(sigmas)
x = newton_krylov(duff_amp_resid,3)
for i, sigma in enumerate(sigmasr):
try:
amplitudesr[i] = newton_krylov(duff_amp_resid,x)
x = amplitudesr[i]
except:
amplitudesr[i] = sp.nan#newton_krylov(duff_amp_resid,0)
x = amplitudesr[i]
plt.plot(sigmasr,amplitudesr)
Out[120]:
In [121]:
plt.plot(sigmasr,amplitudesr)
plt.plot(sigmas,amplitudes)
Out[121]:
In [122]:
def duff_osc2(x, v, params):
omega = params['omega']
t = params['cur_time']
return np.array([[-x-.1*x**3-.1*v+1*sin(omega*t)]])
_,_,_,a,_ = ms.hb_freq(duff_osc2, np.array([[0,1,-1]]), 0.7, num_harmonics=1)
print(a)
In [123]:
_,_,_,a,_ = ms.hb_freq(lambda x,v, params:np.array([[-x-.1*x**3-.1*v+1*sin(0.7*params['cur_time'])]]), np.array([[0,1,-1]]), .7, num_harmonics=1)
a
Out[123]:
Two things to note:
n by 1 array.