How to make this code more efficient?

So I'm a bit of a noob with matlab and I'm not very efficient with my code writing, below is one of my scripts that gets the job done however it just seems...long and not very efficient at all which is an issue because I'll have multiple scripts that will follow this bad habit until a new method comes forward, can anyway help me condense this code?

%constants
ConrodLength=0.1375; %m
Stroke=0.082; %m
Throw=0.041; %m
PistonMass=0.28866; %Kg
CrankAngle=(-360:1:359)'; %Degrees
CrankAngleRad=deg2rad(CrankAngle); %Radians
%Angular velocity
w=[104.7198;130.8997;157.0796;183.2596;209.4395;235.6194;261.7994;314.1593;366.5191;418.879;471.2389;523.5988;575.9587;575.9587;628.3185;680.6784];

%Inertia
for i=CrankAngleRad
        acc1=Throw*(w(1)^2)*cos(i);
        acc2=(Throw/ConrodLength)*Throw*(w(1)^2)*cos(2*i);
        acct=acc1+acc2;
        Inertia1000=PistonMass*acct;
        figure
        subplot(1,1,1)
        plot(CrankAngle,Inertia1000)
        title('Inertia/Crank Angle')
        xlabel('Crank Angle (degrees)')
        ylabel('Inertia (Kg rad/s^2)')
        hold on
end

for i=CrankAngleRad
        acc1=Throw*(w(2)^2)*cos(i);
        acc2=(Throw/ConrodLength)*Throw*(w(2)^2)*cos(2*i);
        acct=acc1+acc2;
        Inertia1250=PistonMass*acct;
        plot(CrankAngle,Inertia1250)
        hold on
end
for i=CrankAngleRad
        acc1=Throw*(w(3)^2)*cos(i);
        acc2=(Throw/ConrodLength)*Throw*(w(3)^2)*cos(2*i);
        acct=acc1+acc2;
        Inertia1500=PistonMass*acct;
        plot(CrankAngle,Inertia1500)
        hold on
end
for i=CrankAngleRad
        acc1=Throw*(w(4)^2)*cos(i);
        acc2=(Throw/ConrodLength)*Throw*(w(4)^2)*cos(2*i);
        acct=acc1+acc2;
        Inertia1750=PistonMass*acct;
        plot(CrankAngle,Inertia1750)
        hold on
end
for i=CrankAngleRad
        acc1=Throw*(w(5)^2)*cos(i);
        acc2=(Throw/ConrodLength)*Throw*(w(5)^2)*cos(2*i);
        acct=acc1+acc2;
        Inertia2000=PistonMass*acct;
        plot(CrankAngle,Inertia2000)
        hold on
end
for i=CrankAngleRad
        acc1=Throw*(w(6)^2)*cos(i);
        acc2=(Throw/ConrodLength)*Throw*(w(6)^2)*cos(2*i);
        acct=acc1+acc2;
        Inertia2250=PistonMass*acct;
        plot(CrankAngle,Inertia2250)
        hold on
end
for i=CrankAngleRad
        acc1=Throw*(w(7)^2)*cos(i);
        acc2=(Throw/ConrodLength)*Throw*(w(7)^2)*cos(2*i);
        acct=acc1+acc2;
        Inertia2500=PistonMass*acct;
        plot(CrankAngle,Inertia2500)
        hold on
end
for i=CrankAngleRad
        acc1=Throw*(w(8)^2)*cos(i);
        acc2=(Throw/ConrodLength)*Throw*(w(8)^2)*cos(2*i);
        acct=acc1+acc2;
        Inertia3000=PistonMass*acct;
        plot(CrankAngle,Inertia3000)
        hold on
end
for i=CrankAngleRad
        acc1=Throw*(w(9)^2)*cos(i);
        acc2=(Throw/ConrodLength)*Throw*(w(9)^2)*cos(2*i);
        acct=acc1+acc2;
        Inertia3500=PistonMass*acct;
        plot(CrankAngle,Inertia3500)
        hold on
end
for i=CrankAngleRad
        acc1=Throw*(w(10)^2)*cos(i);
        acc2=(Throw/ConrodLength)*Throw*(w(10)^2)*cos(2*i);
        acct=acc1+acc2;
        Inertia4000=PistonMass*acct;
        plot(CrankAngle,Inertia4000)
        hold on
end
for i=CrankAngleRad
        acc1=Throw*(w(11)^2)*cos(i);
        acc2=(Throw/ConrodLength)*Throw*(w(11)^2)*cos(2*i);
        acct=acc1+acc2;
        Inertia4500=PistonMass*acct;
        plot(CrankAngle,Inertia4500)
        hold on
end
for i=CrankAngleRad
        acc1=Throw*(w(12)^2)*cos(i);
        acc2=(Throw/ConrodLength)*Throw*(w(12)^2)*cos(2*i);
        acct=acc1+acc2;
        Inertia5000=PistonMass*acct;
        plot(CrankAngle,Inertia5000)
        hold on
end
for i=CrankAngleRad
        acc1=Throw*(w(13)^2)*cos(i);
        acc2=(Throw/ConrodLength)*Throw*(w(13)^2)*cos(2*i);
        acct=acc1+acc2;
        Inertia5500=PistonMass*acct;
        plot(CrankAngle,Inertia5500)
        hold on
end
for i=CrankAngleRad
        acc1=Throw*(w(14)^2)*cos(i);
        acc2=(Throw/ConrodLength)*Throw*(w(14)^2)*cos(2*i);
        acct=acc1+acc2;
        Inertia6000=PistonMass*acct;
        plot(CrankAngle,Inertia6000)
        hold on
end
for i=CrankAngleRad
        acc1=Throw*(w(15)^2)*cos(i);
        acc2=(Throw/ConrodLength)*Throw*(w(15)^2)*cos(2*i);
        acct=acc1+acc2;
        Inertia6500=PistonMass*acct;
        plot(CrankAngle,Inertia6500)
        hold on
end
hold off

Answer 1

[Star Strider]

You can take advantage of the bsxfun function to do all the matrix multiplications at once. Then do the scalar multiplications and plot the results.

See if this does what you want:

%constants
ConrodLength=0.1375; %m
Stroke=0.082; %m
Throw=0.041; %m
PistonMass=0.28866; %Kg
CrankAngle=(-360:1:359)'; %Degrees
CrankAngleRad=deg2rad(CrankAngle); %Radians
%Angular velocity
w=[104.7198;130.8997;157.0796;183.2596;209.4395;235.6194;261.7994;314.1593;366.5191;418.879;471.2389;523.5988;575.9587;575.9587;628.3185;680.6784];

acc1 = Throw * bsxfun(@times, w.^2, cos(CrankAngleRad)');
acc2 = Throw^2/ConrodLength * bsxfun(@times, w.^2, cos(2*CrankAngleRad)');
acct = acc1 + acc2;
Inertia = PistonMass * acct;

figure(1)
plot(CrankAngle, Inertia)
title('Inertia/Crank Angle')
xlabel('Crank Angle (degrees)')
ylabel('Inertia (Kg rad/s^2)')

I didn’t time it, although I’d be surprised if it wasn’t significantly faster. The plot looks the same to me as your original code produced.

—————

EDIT — Thinking about this further, even the bsxfun calls aren’t necessary. Straightforward vector multiplication creates the required matrices:

acc1 = Throw * (w.^2 * cos(CrankAngleRad)');
acc2 = Throw^2/ConrodLength * (w.^2 * cos(2*CrankAngleRad)');

This is likely even faster.

Answer 2

[Jonathan Chin]

   acct1=Throw*bsxfun(@times,cos(CrankAngleRad)*ones(size(w)).',w.^2.');
   acct2= Throw/ConrodLength*Throw*bsxfun(@times,cos(2*CrankAngleRad)*ones(size(w)).',w.^2.');
   acct=acct1+acct2;
   Inertia=PistonMass*acct;

   plot(CrankAngleRad,Inertia(:,1:15))

%%old code
 for i=CrankAngleRad % <- this  is a column vector
 acct1=Throw*w(1)^2*cos(i)
 %rest of code
 end

I replaced cos(i) with cos(CrankAngleRad) because these are the same thing due to CrankAngleRad being a column vector and not a row vector when you assign it in the for loop

I use bsxfun to do all your multiplication at once instead of having to create a line for each

cos(CrankAngleRad)*ones(size(w)).'

this creates a 720 x 16 array. This is essentially the same as cos(CrankAngleRad) but repeated 16 times and concatenated. This is done so you can do an element wise multiplication with w

bsxfun(@times,cos(CrankAngleRad)*ones(size(w)).',w.^2.');

Element wise multiply with w squared.

Then i just repeat this concept for acct2 and you have all the values you need to calculate Inertia for all 16 values of w

plotting is done by plotting the Inertia matrix against your angles