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)
switch acce
case 'average'
beta = 1/4; gamma = 1/2;
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;
