%-----------------------------------------------------------------
% Script 'Discretization Methods', Example 4.7
% (c) 12/2004 by A. Schmidt, IAM, Uni Stuttgart
%-----------------------------------------------------------------
clear all

% Variables
c = 1.0;                            % spring stiffness
m = 1.0;                            % mass
f = 0.0;                            % external force
u0 = 1.0;                           % initial condition displacement
v0 = 0.0;                           % initial condition velocity
a0 = 0.0;                           % initial condition acceleration
t_max = 100.0;                      % time under consideration
dt = 0.1;                           % time step size
n_a = 500;                          % number of sample points
                                    % (analytical solution)
delta = 0.5;                        % Newmark parameter #1
alpha = 0.25;                       % Newmark parameter #2

% Parameters
w = sqrt(c/m);                      % angular frequency

%-----------------------------------------------------------------
% Calculation / time integration using Newmark's method
t(1) = 0.0;                         % first value (time)
u(1) = u0;                          % initial value (displacement)
vel = v0;                           % needed in first time integration step
acc = a0;                           % needed in first time integration step

for i=1:1:t_max/dt
  t(i+1) = i*dt;                    % actual time
                                    % actual disp from Eq. (4.252) with Eq.(4.257)
  u(i+1) = (f+m*(u(i)/alpha/dt^2+vel/alpha/dt+acc*(1.0/2.0/alpha-1.0)))/(m/alpha/dt^2+c);
  acc_new = (u(i+1)-u(i))/alpha/dt^2-vel/alpha/dt-acc*(1.0/2.0/alpha-1.0);      % Eq. (4.257)
  vel_new = vel+((1.0-delta)*acc+delta*acc_new)*dt;                             % Eq. (4.253)

  vel = vel_new;                    % update of velocities for next inc.
  acc = acc_new;                    % update of accelerations for next inc.
end

%-----------------------------------------------------------------
% Analytical solution
t_a(1) = 0.0;                       % first value (time)
u_a(1) = u0;                        % initial value (displacement)
for i=1:1:n_a
  t_a(i+1) = i*t_max/n_a;
  u_a(i+1) = u0*cos(w*t_a(i+1));    % Eq. (3.117)
end

%-----------------------------------------------------------------
% Output / plot
figure(1)
plot(t,u,'r',t_a,u_a,'k')
xlabel('Time / s')
ylabel('Displacement / m')
title('red: Newmark`s method      black: analytical')