# Plotting a trigonometric function appears to be too slow

This calculates the mobility and it's derivative of one element. If I calculate with 1 value of w, it's very fast.

%**setting the equation to calculate**
syms m w k e mi c w t L z ;
a = (((m * w ^ 2 - k) ^ 2 + (c * w) ^ 2) / (e * mi) ^ 2) ^ (1 / 8);
B = 1 / 4 * atan(-c * w / (m * w ^ 2 - k)) + 3.14159 / 4;
X = a * cos(B);
Y = a * sin(B);
ap = 1 / (8 * (e * mi) ^ .25) * (4 * m * w * (m * w ^ 2 - k) + 2 * c ^ 2 * w) * ((m * w ^ 2 - k) ^ 2 + (c * w) ^ 2) ^ (-7 / 8);
bp = 1 / 4 * (c * k + c * m * w ^ 2) / ((c * w) ^ 2 + (m * w ^ 2 - k) ^ 2);
xp = ap * cos(B) - a * bp * sin(B);
yp = ap * sin(B) + a * bp * cos(B);
sab = cos(L * X) / 2 * (exp(L * Y) + exp(-L * Y));
sa = cos(L * Y) / 2 * (exp(L * X) + exp(-L * X));
sbb = sin(L * X) / 2 * (exp(L * Y) - exp(-L * Y));
sb = sin(L * Y) / 2 * (exp(L * X) - exp(-L * X));
sabp = -L * xp / 2 * sin(L * X) * (exp(L * Y) + exp(-L * Y)) + L * yp / 2 * cos(L * X) * (exp(L * Y) - exp(-L * Y));
sap = -L * yp / 2 * sin(L * Y) * (exp(L * X) + exp(-L * X)) + L * xp / 2 * cos(L * Y) * (exp(L * X) - exp(-L * X));
sbbp = L * xp / 2 * cos(L * X) * (exp(L * Y) - exp(-L * Y)) + L * yp / 2 * sin(L * X) * (exp(L * Y) + exp(-L * Y));
sbp = L * yp / 2 * cos(L * Y) * (exp(L * X) - exp(-L * X)) + L * xp / 2 * sin(L * Y) * (exp(L * X) + exp(-L * X));
re1 = -cos(L * Y) / 2 * (exp(L * X) - exp(-L * X)) + sin(L * X) / 2 * (exp(L * Y) + exp(-L * Y));
re4 = cos(z * Y) / 2 * (exp(z * X) - exp(-z * X)) + sin(z * X) / 2 * (exp(z * Y) + exp(-z * Y));
re2 = cos(z * Y) / 2 * (exp(z * X) + exp(-z * X)) + cos(z * X) / 2 * (exp(z * Y) + exp(-z * Y));
re3 = cos(L * Y) / 2 * (exp(L * X) + exp(-L * X)) - cos(L * X) / 2 * (exp(L * Y) + exp(-L * Y));
im1 = -sin(L * Y) / 2 * (exp(L * X) + exp(-L * X)) + cos(L * X) / 2 * (exp(L * Y) - exp(-L * Y));
im4 = sin(z * Y) / 2 * (exp(z * X) + exp(-z * X)) + cos(z * X) / 2 * (exp(z * Y) - exp(-z * Y));
im2 = sin(z * Y) / 2 * (exp(z * X) - exp(-z * X)) - sin(z * X) / 2 * (exp(z * Y) - exp(-z * Y));
im3 = sin(L * Y) / 2 * (exp(L * X) - exp(-L * X)) + sin(L * X) / 2 * (exp(L * Y) - exp(-L * Y));
re1p = L / 2 * cos(L * X) * xp * (exp(L * Y) + exp(-L * Y)) + L / 2 * sin(L * X) * yp * (exp(L * Y) - exp(-L * Y)) + L * yp / 2 * sin(L * Y) * (exp(L * X) - exp(-L * X)) - L / 2 * xp * cos(L * Y) * (exp(L * X) + exp(-L * X));
re2p = z / 2 * cos(z * Y) * xp * (exp(z * X) + exp(-z * X)) - z / 2 * sin(z * Y) * yp * (exp(z * X) + exp(-z * X)) + z * yp / 2 * cos(z * X) * (exp(z * Y) - exp(-z * Y)) - z / 2 * xp * sin(z * X) * (exp(z * Y) + exp(-z * Y));
re3p = L / 2 * cos(L * Y) * xp * (exp(L * X) - exp(-L * X)) - L / 2 * sin(L * Y) * yp * (exp(L * X) + exp(-L * X)) + L * xp / 2 * sin(L * X) * (exp(L * Y) + exp(-L * Y)) - L / 2 * yp * cos(L * X) * (exp(L * Y) - exp(-L * Y));
re4p = z / 2 * cos(z * X) * xp * (exp(z * Y) + exp(-z * Y)) + z / 2 * sin(z * X) * yp * (exp(z * Y) - exp(-z * Y)) + z * xp / 2 * cos(z * Y) * (exp(z * X) + exp(-z * X)) - z / 2 * yp * sin(z * Y) * (exp(z * X) - exp(-z * X));
im1p = -L / 2 * sin(L * X) * xp * (exp(L * Y) - exp(-L * Y)) + L / 2 * cos(L * X) * yp * (exp(L * Y) + exp(-L * Y)) - L * yp / 2 * cos(L * Y) * (exp(L * X) + exp(-L * X)) - L / 2 * xp * sin(L * Y) * (exp(L * X) - exp(-L * X));
im2p = z / 2 * cos(z * Y) * yp * (exp(z * X) - exp(-z * X)) + z / 2 * sin(z * Y) * xp * (exp(z * X) + exp(-z * X)) - z * xp / 2 * cos(z * X) * (exp(z * Y) - exp(-z * Y)) - z / 2 * yp * sin(z * X) * (exp(z * Y) + exp(-z * Y));
im3p = L / 2 * cos(L * Y) * yp * (exp(L * X) - exp(-L * X)) + L / 2 * sin(L * Y) * xp * (exp(L * X) + exp(-L * X)) + L * xp / 2 * cos(L * X) * (exp(L * Y) - exp(-L * Y)) + L / 2 * yp * sin(L * X) * (exp(L * Y) + exp(-L * Y));
im4p = -z / 2 * sin(z * X) * xp * (exp(z * Y) - exp(-z * Y)) + z / 2 * cos(z * X) * yp * (exp(z * Y) + exp(-z * Y)) + z * yp / 2 * cos(z * Y) * (exp(z * X) + exp(-z * X)) + z / 2 * xp * sin(z * Y) * (exp(z * X) - exp(-z * X));
re5 = re1 * re2 - im1 * im2;
re6 = re3 * re4 - im3 * im4;
im5 = re1 * im2 + im1 * re2;
im6 = re3 * im4 + im3 * re4;
re5p = re1p * re2 + re1 * re2p - im1p * im2 - im1 * im2p;
re6p = re3p * re4 + re3 * re4p - im3p * im4 - im3 * im4p;
im5p = re1p * im2 + re1 * im2p + im1p * re2 + im1 * re2p;
im6p = re3p * im4 + re3 * im4p + im3p * re4 + im3 * re4p;
b1 = ((re5 + re6) ^ 2 + (im5 + im6) ^ 2) ^ .5;
b2 = ((1 - sa * sab - sb * sbb) ^ 2 + (sb * sab - sa * sbb) ^ 2) ^ .5;
b1p = 2 * (re5 + re6) * (re5p + re6p) + 2 * (im5 + im6) * (im5p + im6p);
b2p = 2 * (1 - sa * sab - sb * sbb) * (-sab * sap - sa * sabp - sbp * sbb - sb * sbbp) + 2 * (sb * sab - sa * sbb) * (sbp * sab + sb * sabp - sap * sbb - sa * sbbp);
v= w * b1 / 2 / e / mi / a ^ 3 / b2;
f = 1 / 2 / e / mi / a ^ 4 / b2 ^ 2 * ((b1 + w * b1p / 2 / b1) * (a * b2) - (3 * ap * b2 + a * b2p / 2 / b2) * (w * b1));

**%Input parameter to calculate**
e=2.7*10^10;
L=21;
mi=0.001332;
z=5;
m=500;
k=1*10^7;
c=30000;
w=300;

%**calculate the v and f**
resultv = double(subs(v));
resultf = double(subs(f));


However, If I want to calculate velocity and accelerator with a long range of w, such as w(1:1:300), it take a lot of time, more than 1 hour. The code I add is:

w=linspace(0,1,300);
resultv(w) = double(subs(v));
resultf(w) = double(subs(f));
subplot(2,1,1);
plot(w,resultv)
subplot(2,1,2);
plot(w,resultf)


Can anyone tell me why? And how can I fix it to increase its speed?

Second way:

e=2.7*10^10;
L=21;
mi=0.001332;
z=5;
m=500;
k=1*10^7;
c=30000;
for w=1:300
a = (((m * w ^ 2 - k) ^ 2 + (c * w) ^ 2) / (e * mi) ^ 2) ^ (1 / 8);
B = 1 / 4 * atan(-c * w / (m * w ^ 2 - k)) + 3.14159 / 4;
X = a * cos(B);
Y = a * sin(B);
ap = 1 / (8 * (e * mi) ^ .25) * (4 * m * w * (m * w ^ 2 - k) + 2 * c ^ 2 * w) * ((m * w ^ 2 - k) ^ 2 + (c * w) ^ 2) ^ (-7 / 8);
bp = 1 / 4 * (c * k + c * m * w ^ 2) / ((c * w) ^ 2 + (m * w ^ 2 - k) ^ 2);
xp = ap * cos(B) - a * bp * sin(B);
yp = ap * sin(B) + a * bp * cos(B);
sab = cos(L * X) / 2 * (cosh(L*Y));
sa = cos(L * Y) / 2 * (cosh(L*X));
sbb = sin(L * X) / 2 * (sinh(L*Y));
sb = sin(L * Y) / 2 * (sinh(L*X));
sabp = -L * xp / 2 * sin(L * X) * (cosh(L*Y)) + L * yp / 2 * cos(L * X) * (sinh(L*Y));
sap = -L * yp / 2 * sin(L * Y) * (cosh(L*X)) + L * xp / 2 * cos(L * Y) * (sinh(L*X));
sbbp = L * xp / 2 * cos(L * X) * (sinh(L*Y)) + L * yp / 2 * sin(L * X) * (cosh(L*Y));
sbp = L * yp / 2 * cos(L * Y) * (sinh(L*X)) + L * xp / 2 * sin(L * Y) * (cosh(L*X));
re1 = -cos(L * Y) / 2 * (sinh(L*X)) + sin(L * X) / 2 * (cosh(L*Y));
re4 = cos(z * Y) / 2 * (sinh(z*X)) + sin(z * X) / 2 * (cosh(z*Y));
re2 = cos(z * Y) / 2 * (cosh(z*X)) + cos(z * X) / 2 * (cosh(z*Y));
re3 = cos(L * Y) / 2 * (cosh(L*X)) - cos(L * X) / 2 * (cosh(L*Y));
im1 = -sin(L * Y) / 2 * (cosh(L*X)) + cos(L * X) / 2 * (sinh(L*Y));
im4 = sin(z * Y) / 2 * (cosh(z*X)) + cos(z * X) / 2 * (sinh(z*Y));
im2 = sin(z * Y) / 2 * (sinh(z*X)) - sin(z * X) / 2 * (sinh(z*Y));
im3 = sin(L * Y) / 2 * (sinh(L*X)) + sin(L * X) / 2 * (sinh(L*Y));
re1p = L / 2 * cos(L * X) * xp * (cosh(L*Y)) + L / 2 * sin(L * X) * yp * (sinh(L*Y)) + L * yp / 2 * sin(L * Y) * (sinh(L*X)) - L / 2 * xp * cos(L * Y) * (cosh(L*X));
re2p = z / 2 * cos(z * Y) * xp * (cosh(z*X)) - z / 2 * sin(z * Y) * yp * (cosh(z*X)) + z * yp / 2 * cos(z * X) * (sinh(z*Y)) - z / 2 * xp * sin(z * X) * (cosh(z*Y));
re3p = L / 2 * cos(L * Y) * xp * (sinh(L*X)) - L / 2 * sin(L * Y) * yp * (cosh(L*X)) + L * xp / 2 * sin(L * X) * (cosh(L*Y)) - L / 2 * yp * cos(L * X) * (sinh(L*Y));
re4p = z / 2 * cos(z * X) * xp * (cosh(z*Y)) + z / 2 * sin(z * X) * yp * (sinh(z*Y)) + z * xp / 2 * cos(z * Y) * (cosh(z*X)) - z / 2 * yp * sin(z * Y) * (sinh(z*X));
im1p = -L / 2 * sin(L * X) * xp * (sinh(L*Y)) + L / 2 * cos(L * X) * yp * (cosh(L*Y)) - L * yp / 2 * cos(L * Y) * (cosh(L*X)) - L / 2 * xp * sin(L * Y) * (sinh(L*X));
im2p = z / 2 * cos(z * Y) * yp * (sinh(z*X)) + z / 2 * sin(z * Y) * xp * (cosh(z*X)) - z * xp / 2 * cos(z * X) * (sinh(z*Y)) - z / 2 * yp * sin(z * X) * (cosh(z*Y));
im3p = L / 2 * cos(L * Y) * yp * (sinh(L*X)) + L / 2 * sin(L * Y) * xp * (cosh(L*X)) + L * xp / 2 * cos(L * X) * (sinh(L*Y)) + L / 2 * yp * sin(L * X) * (cosh(L*Y));
im4p = -z / 2 * sin(z * X) * xp * (sinh(z*Y)) + z / 2 * cos(z * X) * yp * (cosh(z*Y)) + z * yp / 2 * cos(z * Y) * (cosh(z*X)) + z / 2 * xp * sin(z * Y) * (sinh(z*X));

re5 = re1 * re2 - im1 * im2;
re6 = re3 * re4 - im3 * im4;
im5 = re1 * im2 + im1 * re2;
im6 = re3 * im4 + im3 * re4;
re5p = re1p * re2 + re1 * re2p - im1p * im2 - im1 * im2p;
re6p = re3p * re4 + re3 * re4p - im3p * im4 - im3 * im4p;
im5p = re1p * im2 + re1 * im2p + im1p * re2 + im1 * re2p;
im6p = re3p * im4 + re3 * im4p + im3p * re4 + im3 * re4p;
b1 = ((re5 + re6) ^ 2 + (im5 + im6) ^ 2) ^ .5;
b2 = ((1 - sa * sab - sb * sbb) ^ 2 + (sb * sab - sa * sbb) ^ 2) ^ .5;
b1p = 2 * (re5 + re6) * (re5p + re6p) + 2 * (im5 + im6) * (im5p + im6p);

b2p = 2 * (1 - sa * sab - sb * sbb) * (-sab * sap - sa * sabp - sbp * sbb - sb * sbbp) + 2 * (sb * sab - sa * sbb) * (sbp * sab + sb * sabp - sap * sbb - sa * sbbp);

v(w)= w * b1 / 2 / e / mi / a ^ 3 / b2;

f(w) = 1 / 2 / e / mi / a ^ 4 / b2 ^ 2 * ((b1 + w * b1p / 2 / b1) * (a * b2) - (3 * ap * b2 + a * b2p / 2 / b2) * (w * b1));
drawnow;
subplot(2,1,1);
plot(v)
subplot(2,1,2);
plot(f)

end

• Please post some context. What exactly is it calculating, we cannot offer you a improvement unless we can be sure we are calculating the right thing. – Emily L. Apr 28 '15 at 7:10
• This code is used for calculate the mobility and it's derivative of one element in the ground. a,s,sab,sabp,re,im... is cutted from the long equation to solved it better and easier. and this long equation is difficult to express here. but all the equation in my code solved the requirement correctly. However when I need it to calculate a list of "resultv" and "resultf " respectively with one value of "w" from 1-300 then plot it into the grahp, it take a lot of time to calculate (more than 1 hour). This is my problem. – Vu Ngo Thanh Apr 28 '15 at 10:34
• No, it's not very fast for one element, simply you don't notice any slowdown on a such a small input. – Stefano Sanfilippo Apr 28 '15 at 10:40
• How did you "measure" that is is fast for a single parameter set? When I'm not mistaken, then feeding in constants when calling subs will result in the equation being computed as a constant expression - which means it is largely evaluated at compile time. Have you tried to look at how long the term gets when w isn't set yet? – Ext3h Apr 28 '15 at 10:43
• Mr Stefano & Mr @Ext3h : I just find the other way to solved it. this is much more faster. Someone can explain the reason why this way is faster. update in the question – Vu Ngo Thanh Apr 28 '15 at 11:02

• You'd rather assign values to common functions, eg. you calculate exp(-L * Y) at least three times – Voitcus Apr 28 '15 at 10:33