非線形自由振動 - how to code something

非線形自由振動
$\ddot x +\omega^2 x = -\varepsilon x^3$
初期条件
$x=a, \dot x = 0 (t=0)$
とする

解析解
f:id:seinzumtode:20200504150401p:plain

clear; close;

t=0:0.1:10
omega=1
e=0.1

for a=[1,5,10,11,12]

x = a*cos(omega*t)+...
    e*a^3/(32*omega^2)*(cos(3*omega*t)-cos(omega*t))+... 
    e^2*a^5./(1024*omega^4)*(cos(5*omega*t)-cos(omega*t))+...
    e^3*a^7/(32768*omega^6)*(5*cos(omega*t)-18*cos(3*omega*t)+...
    cos(7*omega*t))

plot(t,x,'DisplayName',sprintf("a=%d",a));
hold on

end

legend();

初期位置aが非線形性をもたらしている
f:id:seinzumtode:20200504150404p:plain