%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 抵抗力が粘性減衰による場合の減衰振動方程式
% 制動係数の大きさの違い
% k = 0.5の場合(減衰が弱い場合)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear ;

tspan = [0, 10] ;
x0 = [0.1; 0] ;
k = 0.5 ;
omega0 = sqrt(40) ;
ode = @(t, x) VdVibration(t, x, k, omega0) ;
[t, x] = ode45(ode, tspan, x0) ;

figure

subplot(1, 2, 1)
plot(t, x(:,1))
xlabel('t')
ylabel('x')
title('Viscous Damped Vibration Equation, k = 0.5, x0 = 0.1, Time history')

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








%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 抵抗力が粘性減衰による場合の減衰振動方程式
% 制動係数の大きさの違い
% k = 2*ω0の場合(臨界制動の場合)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear ;

tspan = [0, 10] ;
x0 = [0.1; 0] ;
k = 2*sqrt(40) ;
omega0 = sqrt(40) ;
ode = @(t, x) VdVibration(t, x, k, omega0) ;
[t, x] = ode45(ode, tspan, x0) ;

figure

subplot(1, 2, 1)
plot(t, x(:,1))
xlabel('t')
ylabel('x')
title('Viscous Damped Vibration Equation, k = 12.649, x0 = 0.1, Time history')

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








%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 抵抗力が粘性減衰による場合の減衰振動方程式
% 制動係数の大きさの違い
% k = 20の場合(減衰が強い場合)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear ;

tspan = [0, 10] ;
x0 = [0.1; 0] ;
k = 20 ;
omega0 = sqrt(40) ;
ode = @(t, x) VdVibration(t, x, k, omega0) ;
[t, x] = ode45(ode, tspan, x0) ;

figure

subplot(1, 2, 1)
plot(t, x(:,1))
xlabel('t')
ylabel('x')
title('Viscous Damped Vibration Equation, k = 20, x0 = 0.1, Time history')

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