How to improve the speed of this Newmark Algorithm?
1 回表示 (過去 30 日間)
古いコメントを表示
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;
3 件のコメント
Jan
2016 年 5 月 16 日
Please add some code for calling your function with realistic inputs, e.g. obtained by rand.
回答 (0 件)
参考
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!