Hi, I've also been attempting these problems from Chopra (2014). I'm having some trouble following your syntax and would need more descriptive %documentation. I have not attempted the constant stiffness (modified) Newton-Raphson method yet, but here is a code I made that works for generic Newton-Raphson. It reproduced the results from Chopra, Example 5.5 exactly.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%Newmark's Method for SDF Nonlinear Systems %%%
By: Christopher Wong %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%Method as per A.K. Chopra %%%
Example 5.5, (2014) %%%
%%%Using Newton-Raphson iterations
% Clear workspace for each new run
clear
% Establish time step parameters
T_n=1; % natural period
dt=0.1; % time step size
t=0:dt:1.0; % total length of time
n=length(t)-1; % number of time steps
TOL=1e-3; % convergence criteria
% Determine which special case to use: constant avg. vs. linear accel
if dt/T_n<=0.551 % Use linear accel. method - closest to theoretical
gamma=0.5;
beta=1/6;
else % Use constant avg. accel. method - unconditionally stable (Example 5.5)
gamma=0.5;
beta=0.25;
end
%%%Establish system properties
xi=0.05; % percentage of critical damping
omega_n=2*pi/T_n; % natural angular frequency
m=0.4559; % mass
k=18; % stiffness
c=2*xi*m*omega_n; % damping constant
u_y=2; % yield deformation
%%%Input excitation function
p(t<0.6)=50*sin(pi*t(t<0.6)/0.6);
p(t>=0.6)=0;
%%%Establish initial conditions @ i=1
u(1)=0; % displacement
v(1)=0; % velocity
f_s(1)=0; % restoring force
k_T(1)=k; % tangent stiffness
a(1)=(p(1)-c*v(1)-f_s(1))/m; % acceleration
%%%Calculate Newmark constants
A1=m/(beta*dt.^2)+gamma*c/(beta*dt);
A2=m/(beta*dt)+c*(gamma/beta-1);
A3=m*(1/(2*beta)-1)+dt*c*(gamma/(2*beta)-1);
k_hat=k+A1;
%%%Calculations for each time step, i=0,1,2,...,n
%%%Inititialize time step
for i=1:n
p_hat(i+1)=p(i+1)+A1*u(i)+A2*v(i)+A3*a(i); % Chopra eqn. 5.4.12
u(i+1)=p_hat(i+1)/k_hat; % linear displacement
k_T(i+1)=k_T(i); % tangent stifness at beginning of time step
f_s(i+1)=u(i+1)*k_T(i+1); % linear restoring force
%%%Determine if iteration is linear or nonlinear
if abs(f_s(i+1))<=abs(k*u_y) % If linear
u(i+1)=u(i+1); % keep linear value
f_s(i+1)=f_s(i+1); % keep linear value
else % If nonlinear
u(i+1)=u(i); % restore value from previous nonlinear iteration
f_s(i+1)=f_s(i); % resore value from previous nonlinear iteration
end
R_hat(i+1)=p_hat(i+1)-f_s(i+1)-A1*u(i+1); % Compute initial residual
%%%Begin Netwon-Raphson iteration
while abs(R_hat(i+1))>TOL % Terminate if converged
k_T_hat(i+1)=k_T(i+1)+A1;
du=R_hat(i+1)/k_T_hat(i+1);
u(i+1)=u(i+1)+du;
f_s(i+1)=f_s(i)+k*(u(i+1)-u(i));
%%%Determine if restoring force is yielding
if abs(f_s(i+1))>=abs(k*u_y) % If yielding
f_s(i+1)=k*u_y;
k_T(i+1)=0;
else % If elastic
k_T(i+1)=18;
end
R_hat(i+1)=p_hat(i+1)-f_s(i+1)-A1*u(i+1); % Compute new residual
end
%%%Calculations for new velocity and acceleration
v(i+1)=gamma*(u(i+1)-u(i))/(beta*dt)+v(i)*(1-gamma/beta)+dt*a(i)*(1-gamma/(2*beta));
a(i+1)=(u(i+1)-u(i))/(beta*dt.^2)-v(i)/(beta*dt)-a(i)*(1/(2*beta)-1);
end
%%%Generate Solution Table
solution=[t;p;R_hat;k_T;k_T_hat;u;f_s;v;a];
solution=solution';
colNames={'t','p','R_hat','k_T','k_T_hat','u','f_s','v','a'};
solution=array2table(solution,'VariableNames',colNames);
%%%Generate plots
figure(3)
plot(t,u)
xlabel('\itt\rm, s')
ylabel('\itu\rm, cm')
title('Displacement vs. Time Plot')
figure(4)
plot(t,p)
xlabel('\itt\rm, s')
ylabel('\itp\rm, kN')
title('Excitation Function Plot')
figure(5)
plot(u,f_s)
xlabel('\itu\rm, cm')
ylabel('\itf_s\rm, kn')
title('Elastoplastic Loop Plot')
