# Newmark's Method for Nonlinear Systems

323 views (last 30 days)
Markus Landwehr on 23 May 2017
Commented: Ya Su on 12 Apr 2021
Hello everyone,
based on the book "Dynamics of Structures" by Chopra I would like to simulate nonlinear vibrations in Matlab with the Newmark´s method for nonlinear systems. I attached the book chapter where the algorithm (modified Newton-Raphson and Newmark´s-method) are explained.
My current implementation of these algorithm:
function [u, udot, u2dot] = newmark_int_nlin(t,p,u0,udot0,m,k,c,gamma,beta,Werkstoffmodell,Solver)
%Newmark's Method for nonlinear systems
%--------------------------------------------------------------------------
% Integrates a nonlinear 1-DOF system with mass "m", spring stiffness "k", damping
% coeffiecient "c" and nonlinear force "fs", when subjected to an external load P(t).
% Returns the displacement, velocity and acceleration of the system with
% respect to an inertial frame of reference.
%
% SYNTAX
% [u, udot, u2dot] = newmark_int_nlin(t,p,u0,udot0,m,k,c,gamma,beta,Werkstoffmodell,Solver)
%
% INPUT
% [t] : Time Vector [n,1]
% [p] : Externally Applied Load [n,1]
% [u0]: Initial Position [1,1]
% [udot0]: Initial Velocity [1,1]
% [m]: System Mass [1,1]
% [k]: System Stiffness [1,1]
% [c]: System Damping [1,1]
% [gamma]: Newmark coefficient [1,1]
% [beta]: Newmark coefficient [1,1]
% [Werkstoffmodell]: material model parameters
% [Solver]: solver parameters
%
% OUTPUT
% [u]: Displacemente Response [n,1]
% [udot]: Velocity [n,1]
% [u2dot]: Acceleration [n,1]
%
% N = number of time steps
%
% The options include changing the value of the "gamma" and "beta"
% coefficient which appear in the formulation of the method. By default
% these values are set to gamma = 1/2 and alpha = 1/4.
dt = t(2) - t(1); %timestep
a = m/(beta*dt) + gamma*c/beta; %newmark coefficient
b = 0.5*m/beta + dt*(0.5*gamma/beta - 1)*c; %newmark coefficient
TOL = 1e-6; %Tolerance
j_max = 100; %max iterations
dp = diff(p); %external force
u = zeros(length(t),1); udot = u; u2dot = u;
u(1,1) = u0; %initial condition
udot(1,1) = udot0;%initial condition
u2dot(1,1) = 1/m*(p(1)-k*u0-c*udot0-Fresfun(u(1,1),t(1),Werkstoffmodell,Solver)); %initial condition
% u2dot(1) = 1/m*(p(1)-k*u0-c*udot0);
for i = 1:(length(t)-1) %for each timestep
dp_dach = dp(i) + a*udot(i,1) + b*u2dot(i,1);
% ki = (Fresfun(u(i+1,1),t(i),Werkstoffmodell,Solver)-Fresfun(u(i,1),t(i),Werkstoffmodell,Solver))/(u(i+1,1)-u(i,1));%tangent stiffness
ki = k;
ki_dach = ki + gamma/(beta*dt)*c + m/(beta*dt^2); %tangent stiffness
% modified Newton Raphson iterative procedure
% initialize data
j = 2;
u(i+1,1) = u(i,1);
fs(i,1) = Fresfun(u(i,1),t(i),Werkstoffmodell,Solver);
dR(1) = 0;
dR(2) = dp_dach;
kT = ki_dach;
% calculation for each iteration
while j < j_max
du(j) = (dR(j))/kT;
u(i+1,j) = u(i+1,j-1)+du(j);
fs(i,j) = Fresfun(u(i,j),t(i),Werkstoffmodell,Solver); %compute fs(j)
df(i,j) = fs(i,j)-fs(i,j-1)+(kT-k)*du(j);
dR(j+1) = dR(j)-df(i,j);
du_sum = sum(u,2);
if du(j)/du_sum(j) < TOL % repetition for next iteration
break;
end
j = j+1;
end
du(i) = du_sum(j);
dudot_i = gamma/(beta*dt)*du(i) - gamma/beta*udot(i) + dt*(1-0.5*gamma/beta)*u2dot(i);
du2dot_i = 1/(beta*dt^2)*du(i) - 1/(beta*dt)*udot(i) - 0.5/beta*u2dot(i);
u(i+1) = du(i) + u(i);
udot(i+1) = dudot_i + udot(i);
u2dot(i+1) = du2dot_i + u2dot(i);
end
In this state, the code is not working properly. I am having issues with the "modified newton raphson algorithm" and I think there are mistakes in my code.
I would appreciate if you could check my code and give me feedback.
Best regards Markus

Christopher Wong on 6 Jul 2018
Edited: Christopher Wong on 24 Sep 2018
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')

jiawei jiang on 20 Aug 2017
hello，I'm researching on the newmark-βmethod，have you resolve this problem ,could you send me a code of your。so appreciate。my email is : jjwsmile@163.com

victorroda on 29 Nov 2018
I have worked out this code, but I find that it does not work properly when the spring is unloading within the yield region.
Additional lines must be implemented in the yield condition to satisfy that if fS·dU<0, then the stiffness is no longer 0, but the stiffness of the linear region.
Regards,

civil s on 1 Apr 2020
hi
my code doesnt converge for negative post yield stiffness.can anyone help please? it is ok for zero and positive post yield stiffness but it does not work for negative one.
pls let me know if anyone can help so i will send the code.
thanks
Ya Su on 12 Apr 2021
Have you found the solution for negative stiffness?