%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% マシューの方程式の自由振動
% 係数の角振動数の違いによる変化
% ω = 1.94の場合
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear ;

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

figure

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

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








%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% マシューの方程式の自由振動
% 係数の角振動数の違いによる変化
% ω = 2の場合
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear ;

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

figure

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

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








%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% マシューの方程式の自由振動
% 係数の角振動数の違いによる変化
% ω = 2.06の場合
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear ;

tspan = [0, 340] ;
x0 = [0.1; 0] ;
gamma = 0 ;
omega0 = 1 ;
epsilon = 0.1 ;
omega = 2.06 ;
ode = @(t, x) Mathieu(t, x, gamma, omega0, epsilon, omega) ;
[t, x] = ode45(ode, tspan, x0) ;

figure

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

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








