%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 抵抗力が運動摩擦係数μによる場合の減衰振動方程式
% 初期条件の違い
% x(0) = 0.1の場合
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear ;

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

figure

subplot(1, 2, 1)
plot(t, x(:,1))
xlabel('t')
ylabel('x')
title('Dry Solid Friction Vibration Equation, \mu = 0.01, x0 = 0.1, Time history')

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








%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 抵抗力が運動摩擦係数μによる場合の減衰振動方程式
% 初期条件の違い
% x(0) = -0.1の場合
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear ;

tspan = [0, 15] ;
x0 = [-0.1; 0] ;
mu = 0.01 ;
omega0 = sqrt(40) ;
ode = @(t, x) MuVibration(t, x, mu, omega0) ;
[t, x] = ode45(ode, tspan, x0) ;

figure

subplot(1, 2, 1)
plot(t, x(:,1))
xlabel('t')
ylabel('x')
title('Dry Solid Friction Vibration Equation, \mu = 0.01, x0 = -0.1, Time history')

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








%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 抵抗力が運動摩擦係数μによる場合の減衰振動方程式
% 運動摩擦係数の大きさの違い
% μ = 0.001の場合
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear ;

tspan = [0, 120] ;
x0 = [0.1; 0] ;
mu = 0.001 ;
omega0 = sqrt(40) ;
ode = @(t, x) MuVibration(t, x, mu, omega0) ;
[t, x] = ode45(ode, tspan, x0) ;

figure

subplot(1, 2, 1)
plot(t, x(:,1))
xlabel('t')
ylabel('x')
title('Dry Solid Friction Vibration Equation, \mu = 0.001, x0 = 0.1, Time history')

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








%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 抵抗力が運動摩擦係数μによる場合の減衰振動方程式
% 運動摩擦係数の大きさの違い
% μ = 0.1の場合
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear ;

tspan = [0, 2] ;
x0 = [0.1; 0] ;
mu = 0.1 ;
omega0 = sqrt(40) ;
ode = @(t, x) MuVibration(t, x, mu, omega0) ;
[t, x] = ode45(ode, tspan, x0) ;

figure

subplot(1, 2, 1)
plot(t, x(:,1))
xlabel('t')
ylabel('x')
title('Dry Solid Friction Vibration Equation, \mu = 0.1, x0 = 0.1, Time history')

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