%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ダフィングの方程式の自由振動
% 初期条件の違いによる角振動数の変化
% x(0) = 1の場合
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear ;

% 計算する時間区間　→　tspan
tspan = [0, 30] ;
% 問題の初期値　→　x0
x0 = [1; 0] ;
% ダフィングの方程式のパラメータ
gamma = 0 ;
omega0 = 1 ;
epsilon = 0.1 ;
A = 0 ;
omega = 0 ;
% ダフィングの方程式を1次常微分方程式(ODE)として定義　→　Duffing( )
% 1次常微分方程式(ODE)の関数ハンドル　→　ode
ode = @(t, x) Duffing(t, x, gamma, omega0, epsilon, A, omega) ;
% 1次常微分方程式(ODE)の初期値問題を，陽的ルンゲ・クッタで解く
[t, x] = ode45(ode, tspan, x0) ;

figure

subplot(1, 2, 1)
plot(t, x(:,1))
xlabel('t')
ylabel('x')
title('Duffing Equation, \gamma = 0, \omega0 = 1, \epsilon = 0.1, x0 = 1, Time history')

subplot(1, 2, 2)
plot(x(:, 1), x(:, 2))
xlabel('x')
ylabel('dx/dt')
title('Phase plane plot')








%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ダフィングの方程式の自由振動
% 初期条件の違いによる角振動数の変化
% x(0) = 2の場合
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear ;

tspan = [0, 30] ;
x0 = [2; 0] ;
gamma = 0 ;
omega0 = 1 ;
epsilon = 0.1 ;
A = 0 ;
omega = 0 ;
ode = @(t, x) Duffing(t, x, gamma, omega0, epsilon, A, omega) ;
[t, x] = ode45(ode, tspan, x0) ;

figure

subplot(1, 2, 1)
plot(t, x(:,1))
xlabel('t')
ylabel('x')
title('Duffing Equation, \gamma = 0, \omega0 = 1, \epsilon = 0.1, x0 = 2, Time history')

subplot(1, 2, 2)
plot(x(:, 1), x(:, 2))
xlabel('x')
ylabel('dx/dt')
title('Phase plane plot')








%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ダフィングの方程式の強制振動
% 上田のパラメータによるカオスの発生
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear ;

tspan = [0, 100] ;
x0 = [0; 0] ;
gamma = 0.05/2 ;
omega0 = 0 ;
epsilon = 1 ;
A = 7.5 ;
omega = 1 ;
ode = @(t, x) Duffing(t, x, gamma, omega0, epsilon, A, omega) ;
[t, x] = ode45(ode, tspan, x0) ;

figure

subplot(1, 2, 1)
plot(t, x(:,1))
xlabel('t')
ylabel('x')
title('Duffing Equation, \gamma = 0.05/2, \omega0 = 0, \epsilon = 1, A = 7.5, Time history')

subplot(1, 2, 2)
plot(x(:, 1), x(:, 2))
xlabel('x')
ylabel('dx/dt')
title('Phase plane plot')








%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ダフィングの方程式の強制振動
% 位相平面図における解の収束　→　アトラクタ
% 収束していく軌道　→　吸引リミットサイクル
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear ;

tspan = [0, 100] ;
x0 = [0; 0] ;
gamma = 0.1/2 ;
omega0 = 0 ;
epsilon = 1 ;
A = 12 ;
omega = 1 ;
ode = @(t, x) Duffing(t, x, gamma, omega0, epsilon, A, omega) ;
[t, x] = ode45(ode, tspan, x0) ;

figure

subplot(1, 2, 1)
plot(t, x(:,1))
xlabel('t')
ylabel('x')
title('Duffing Equation, \gamma = 0.1/2, \omega0 = 0, \epsilon = 1, A = 12, Time history')

subplot(1, 2, 2)
plot(x(:, 1), x(:, 2))
xlabel('x')
ylabel('dx/dt')
title('Phase plane plot')
