%% simple harmonic motion

clear; clc; close all;
hold on;
grid on;
axis equal;

%% 参数设置
K=10; %弹簧劲度系数 N/m
m=0.2; %震子质量 Kg

t_total = 10;       % 总时间 (s)
dt = 0.0001;          % 时间步长 (s)

%% 数值模拟（欧拉法）
% 初始化数组
n_steps = round(t_total / dt);
t = linspace(0, t_total, n_steps);
x=[];
v=[];
a=[];

x(1)=-0.2; %从负半轴释放
v(1)=0;
a(1)=-x(1)*K;

for i = 1:n_steps-1

  a(i+1)=-x(i)*K;
  v(i+1)=v(i) + a(i) * dt;
  x(i+1)=x(i)+v(i) *dt;

end

subplot(2,2,1)

plot(t, x, 'b-', 'LineWidth', 1);
title('t-x');
grid on
axis equal

subplot(2,2,2)
plot(t, v, 'r-', 'LineWidth', 1);
title('t-v');
grid on
axis equal

subplot(2,2,3)
plot(t, a, 'g-', 'LineWidth', 1);
title('t-a');
grid on
% axis equal