0
\$\begingroup\$

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
\$\endgroup\$
  • \$\begingroup\$ 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. \$\endgroup\$ – Emily L. Apr 28 '15 at 7:10
  • \$\begingroup\$ 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. \$\endgroup\$ – Vu Ngo Thanh Apr 28 '15 at 10:34
  • \$\begingroup\$ No, it's not very fast for one element, simply you don't notice any slowdown on a such a small input. \$\endgroup\$ – Stefano Sanfilippo Apr 28 '15 at 10:40
  • \$\begingroup\$ 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? \$\endgroup\$ – Ext3h Apr 28 '15 at 10:43
  • \$\begingroup\$ 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 \$\endgroup\$ – Vu Ngo Thanh Apr 28 '15 at 11:02
1
\$\begingroup\$

I notice you've got a lot of common subexpressions there, particularly the exponentials. Would it help to assign them to new variables?

\$\endgroup\$
  • \$\begingroup\$ you mean that, I replace exp(L * Y) + exp(-L * Y) with cosh(L*Y) and the other. I has tried this, but when I calculate for 1 value of w such as: w=300, matlab calculate it very fast. but when I calculate it with more value of w to plot v and f into graph, such as w=(1:1:300) it take a lot of time. more than 1 hour \$\endgroup\$ – Vu Ngo Thanh Apr 28 '15 at 4:10
  • 1
    \$\begingroup\$ You'd rather assign values to common functions, eg. you calculate exp(-L * Y) at least three times \$\endgroup\$ – Voitcus Apr 28 '15 at 10:33
  • \$\begingroup\$ @Voitcus: Thank for your help, I just find the way to solved it \$\endgroup\$ – Vu Ngo Thanh Apr 29 '15 at 5:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.