How to improve the speed of this Newmark Algorithm?

Hello guys,

I wrote this Newmark algorithm based on someone's post to solve a dynamic problem. My main program calls it for iterations. Profiler told me that this Newmark algorithm ran 25334 times and costed most of the computing time. Can somebody help me to see if there are anything I could do to improve the speed of this function?

I wrote it with superscript _r and Phi so that if I have a reduced system, I could use the same function to solve the problem.

Thanks!

A bit more explanation: for input, Phi is the reduced basis, if not computing a reduced system, one can set Phi to identity matrix with same size of M_r. M_r, C_r, K_r are mass, damping, stiffness matrix, respectively. F_r is the force vector in time. acce defines the coefficients beta and gamma, I usually use beta = 1/4, gamma = 1/2 for numerical stability. dT is length of time steps, maxT is the maximum time for the problem, U0 and V0 are initial displacement and velocity, respectively.

function [U_r, V_r, A_r, U, V, A, t, time_step_NO] = NewmarkBetaReducedMethod...
    (Phi, M_r, C_r, K_r, F_r, acce, dT, maxT, U0, V0)
% Solve dynamic problem with FORCE IN TIME.

switch acce
    case 'average'
         beta = 1/4; gamma = 1/2; % al = alpha
    case 'linear'
         beta = 1/6; gamma = 1/2;
end

Time step and initial conditions U0, V0.

t = 0:dT:(maxT);
time_step_NO = maxT/dT;

U_r = zeros(length(K_r), length(t));
U_r(:, 1) = U_r(:, 1)+U0;
V_r = zeros(length(K_r), length(t));
V_r(:, 1) = V_r(:, 1)+V0;
A_r = zeros(length(K_r), length(t));
A_r(:, 1) = A_r(:, 1)+M_r\(F_r(:, 1)-C_r*V_r(:, 1)-K_r*U_r(:, 1));

Coefficients

a1 = gamma/(beta*dT);
a2 = 1/(beta*dT^2);
a3 = 1/(beta*dT);
a4 = gamma/beta;
a5 = 1/(2*beta);
a6 = dT*(gamma/(2*beta)-1);
Khat = K_r+a1*C_r+a2*M_r;
a = a3*M_r+a4*C_r; 
b = a5*M_r+a6*C_r;
for i_nm = 1:length(t)-1

      dFhat = F_r(:, i_nm+1)-F_r(:, i_nm)+a*V_r(:, i_nm)+b*A_r(:, i_nm);
      dU_r = Khat\dFhat;
      dV_r = a1*dU_r-a4*V_r(:, i_nm)-a6*A_r(:, i_nm);
      dA_r = a2*dU_r-a3*V_r(:, i_nm)-a5*A_r(:, i_nm);
      U_r(:, i_nm+1) = U_r(:, i_nm)+dU_r;
      V_r(:, i_nm+1) = V_r(:, i_nm)+dV_r;
      A_r(:, i_nm+1) = A_r(:, i_nm)+dA_r;

end
A_r = A_r(:, 1:size(A_r, 2));
V_r = V_r(:, 1:size(V_r, 2));
U_r = U_r(:, 1:size(U_r, 2));
A = Phi*A_r;
V = Phi*V_r;
U = Phi*U_r;

[Brattv]

Have you tried to implement by doing vector/matrix multiplications. This is often faster compared with loops. I do not have Matlab available on this computer, and do not know the variable sizes in your program - but i would have attempted something like

% Calculating parameters to dFhat
F_rDiff = F_r(:,2:end) .- F_r(:,1:end-1);

dFhat = F_rDiff.-(a*V_r).-(b*A_r);

This is kind of a trail and error process, but it could potentially decrease the content of your for loop - which could result in a lower run-time.

[Jan]

Please add some code for calling your function with realistic inputs, e.g. obtained by rand.

[Xiaohan Du]

Hi Vegard,

Thanks for the reply!

I guess you mean put these lines outside the loop? I'm afraid this cannot be done because this is time integration method, each current step calls the result from the previous step. Therefore columns of dFhat changes every iteration and cannot be computed in one go or doing matrix multiplications.