The following code seems to be taking forever to run. I don't get the reasons behind it though. Any suggestions to make it run quicker are highly appreciated

 %Prob 1.b - Distributed Loading.
n=6; %No. of mesh elements considered
L=1; %Length of total beam
l=L/n;%Length of each element.
c=0.1*L; h=0.1*L;b=1.875; zeta=0.7341;  
Rho=1000; a=c*h; E=12*(10^9); I=c*(h^3)/12;k_spring = E*I/(h*L*L);
m=(Rho*a*l/420)*[156 22*l 54 -13*l; 22*l 4*l*l 13*l -3*l*l; 54 -13*l 156 -22*l; -13*l -3*l*l -22*l 4*l*l];
k=(E*I/(l^3))*[12 6*l -12 6*l; 6*l 4*l*l -6*l 2*l*l; -12 -6*l 12 -6*l; 6*l 2*l*l -6*l 4*l*l];
global A
D=0.1*eye(2*n+2);
M=zeros(2*n+2);
K=zeros(2*n+2);
for i=1:n
    p=2*i-2;
    for j=1:4
        for s=1:4
            M(p+j,p+s)=M(p+j,p+s)+m(j,s); %Global mass and stiffness matrices.
            K(p+j,p+s)=K(p+j,p+s)+k(j,s);
        end
    end
end
M_sub=M(3:length(M),3:length(M));
K_sub=K(3:length(K),3:length(K)); %deleting the first two rows and columns to account for fixed end.
D_sub=D(3:length(D),3:length(D));
q=length(K_sub);
K_sub(q-1,q-1)=K_sub(q-1,q-1)+k_spring;%adding the spring stiffness at the boundary
A = [zeros(q) eye(q); -M_sub\K_sub -M_sub\D_sub];
tspan = linspace(0,20,21);
syms eta
si(1)=1-3*(eta^2)+2*(eta^3); %shape functions used
si(2)=(eta-2*(eta^2)+(eta^3))*l;
si(3)=3*(eta^2)-2*(eta^3);
si(4)=(-(eta^2)+(eta^3))*l;
%g=zeros(n,1); %distributed forcing load per element
for i=1:n
    g(i)=(cos(b*l*(eta+i-1))-cosh(b*l*(eta+i-1))-zeta*(sin(b*l*(eta+i-1))-sinh(b*l*(eta+i-1))));
end
for i=1:4
    for j=1:n
        F(i,j)=int(g(j)*si(i)*l,eta,0,1);
    end
end
F_single(1:2,1)=F(1:2,1);
for i=1:5
    F_single(2*i+1,1)=F(3,i)+F(1,i+1);
    F_single(2*i+2,1)=F(4,i)+F(2,i+1);
end
F_single(13:14,1)=F(3:4,6);
F_sub=F_single(3:14,1);
[t,Q]=ode45(@(t,Q) ganti(t,Q,M_sub,F_sub,A,n,l),tspan,zeros(24,1));
save('prob1b2r')

And the function is

 function dQdt=ganti(t,Q,M_sub,F_sub,A,n,l)
omega=2363.28;
dQdt = zeros(24,1);
F_final=F_sub*sin(omega*t);
U = [zeros(12,1);M_sub\F_final];
for i=1:24
dQdt(i)=A(i,:)*[Q(1);Q(2);Q(3);Q(4);Q(5);Q(6);Q(7);Q(8);Q(9);Q(10);Q(11);Q(12);Q(13);Q(14);Q(15);Q(16);Q(17);Q(18);Q(19);Q(20);Q(21);Q(22);Q(23);Q(24)]+U(i);
end

end

Answer 1

[Walter Roberson]

You do not initialize F, so it is growing on the fly (but it is not large.) More importantly, the right hand side is int() which is for finding closed form integrals, giving exact results such as

(82126976*exp(-5/16))/390625 - (18263296*cos(5/16))/78125 + (9189504*exp(5/16))/390625 + (47337472*sin(5/16))/390625 - 2447/3125

This leads to F_sub being symbolic. You then pass F_sub into your ode function, so much of the math in the ode ends up being done in symbolic form, before ending up getting converted to numeric form.

Considering your A values are pure numeric, you cannot justify calculating with F_sub in symbolic form: you should be passing double(F_sub) to the ode. Or, more likely, instead of using int() you should be using vpaintegral() and assigning into F that has been initialized to double precision.

Note: there is no obvious reason why you made A global.

[Sanjit Ganti]

Thank you very much Walter Roberson. Using double(F_sub) has reduced the run time more than vpaintegral. Thank you for the answer.

Regarding declaring A as global, that was being used for some other code previously and I just forgot to change it in this. Anyways, would declaring A as global have impact on the run time as well?

[Walter Roberson]

I am not sure exactly how MATLAB treats a global that is passed as a parameter. I think it becomes ordinary within the routine it is passed to, in which case performance within that routine should not be affected, but... hmmm... I'm not completely sure how it would work.

Generally speaking globals are the slowest form of variables to look up.

Okay, my tests show that if you alter a global that has also been passed as a parameter, that in the places that know it only as a parameter, the value is not changed. This means that the slowdown of using globals only occurs in the places that it is declared as global (or a shared variable.)