% ch1example1prg1.m
g=9.8;  % 重力加速度
v=0;    % 設定初始速度條件
s=0;    % 設定初始位移條件
t=0;    % 設定起始時間
dt=0.1; % 設置計算步長
N=20;   % 設置仿真遞推次數(shù). 仿真時間等于N與dt的乘積
for k=1:N
  v=v+g*dt;                 % 計算新時刻的速度
  s(k+1)=s(k)+v*dt;         % 新位移
  t(k+1)=t(k)+dt;           % 時間更新
end
% 理論計算, 以便與仿真結果對照
t_theory=0:0.01:N*dt;       % 設置解析計算的時間點
v_theory=g*t_theory;        % 解析計算的瞬時速度
s_theory=1/2*g*t_theory.^2; % 解析計算的瞬時位移
% 作圖: 仿真結果與解析結果對比
t=0:dt:N*dt;
plot(t,s,'o', t_theory,s_theory, '-');
xlabel('時間 t'); ylabel('位移 s');
legend('仿真結果','理論結果');