%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ファンデルポールの方程式，自励発振のみ
% 初期条件の違い
% x(0) = 0.2の場合
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear ;

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

figure

subplot(1, 2, 1)
plot(t, x(:,1))
xlabel('t')
ylabel('x')
title('van der Pol Equation, \epsilon = 0.15, x0 = 0.2, Time history')

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








%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ファンデルポールの方程式，自励発振のみ
% 初期条件の違い
% x(0) = 3の場合
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear ;

tspan = [0, 80] ;
x0 = [3; 0] ;
epsilon = 0.15 ;
omega0 = 1 ;
A = 0 ;
omega = 0 ;
ode = @(t, x) vanderpol(t, x, epsilon, omega0, A, omega) ;
[t, x] = ode45(ode, tspan, x0) ;

figure

subplot(1, 2, 1)
plot(t, x(:,1))
xlabel('t')
ylabel('x')
title('van der Pol Equation, \epsilon = 0.15, x0 = 3, Time history')

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








%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ファンデルポールの方程式，自励発振のみ
% 非線形項の大きさの違い
% ε = 1の場合
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear ;

tspan = [0, 20] ;
x0 = [2; 0] ;
epsilon = 1 ;
omega0 = 1 ;
A = 0 ;
omega = 0 ;
ode = @(t, x) vanderpol(t, x, epsilon, omega0, A, omega) ;
[t, x] = ode45(ode, tspan, x0) ;

figure

subplot(1, 2, 1)
plot(t, x(:,1))
xlabel('t')
ylabel('x')
title('van der Pol Equation, \epsilon = 1, Time history')

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








%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ファンデルポールの方程式，自励発振のみ
% 非線形項の大きさの違い
% ε = 1000の場合(スティッフな場合)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear ;

tspan = [0, 3000] ;
x0 = [2; 0] ;
epsilon = 1000 ;
omega0 = 1 ;
A = 0 ;
omega = 1 ;
ode = @(t, x) vanderpol(t, x, epsilon, omega0, A, omega) ;
% 1次常微分方程式(ODE)の初期値問題(スティッフな問題)を，数値微分法(NDF)に基づく可変次数のソルバーで解く
[t, x] = ode15s(ode, tspan, x0) ;

figure

subplot(1, 2, 1)
plot(t, x(:,1))
xlabel('t')
ylabel('x')
title('van der Pol Equation, \epsilon = 1000, Time history')

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








%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ファンデルポールの方程式，自励発振と強制振動
% ω = 1の場合
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear ;

tspan = [0, 200] ;
x0 = [0; 0] ;
epsilon = 0.1 ;
omega0 = 1 ;
A = 0.1 ;
omega = 1 ;
ode = @(t, x) vanderpol(t, x, epsilon, omega0, A, omega) ;
[t, x] = ode45(ode, tspan, x0) ;

figure

subplot(1, 2, 1)
plot(t, x(:,1))
xlabel('t')
ylabel('x')
title('van der Pol Equation, \epsilon = 0.1, \omega = 1, Time history')

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