Any comment to speed up the speed of caculation of symbolic loops having Legendre polynomials?

14 views (last 30 days)
syms eta__2 zeta__2
II=12;JJ=11;M=22;
Hvs2 = ((5070602400912917605986812821504*(zeta__2 + 2251799813683713/2251799813685248)^2)/2356225 + (9007199254740992*(eta__2 + 2935286035937695/18014398509481984)^2)/196937227765191 - 1)*((81129638414606681695789005144064*(zeta__2 + 9007199254732683/9007199254740992)^2)/69039481 + (576460752303423488*(eta__2 + 3261970163074917/4503599627370496)^2)/6904142590940591 - 1)*((324518553658426726783156020576256*(zeta__2 + 140737488355209/140737488355328)^2)/231983361 + (144115188075855872*(eta__2 - 262292457514301/562949953421312)^2)/2637878570603985 - 1)*((144115188075855872*(zeta__2 + 4028041154330599/4503599627370496)^2)/424643881623313 + eta__2^2 - 1)*((20282409603651670423947251286016*(zeta__2 - 4503599627213111/4503599627370496)^2)/24770038225 + (288230376151711744*(eta__2 - 7530397878711147/9007199254740992)^2)/5204731445635785 - 1)*((324518553658426726783156020576256*(zeta__2 + 4503599627365785/4503599627370496)^2)/355058649 + (36893488147419103232*(eta__2 + 4434826747744735/4503599627370496)^2)/8603290501959015 - 1)*((4611686018427387904*(eta__2 + 2213733699584161/2251799813685248)^2)/1317884237102575 + (324518553658426726783156020576256*(zeta__2 - 4503599627284663/4503599627370496)^2)/117876175561 - 1)*((81129638414606681695789005144064*(zeta__2 + 9007199254735975/9007199254740992)^2)/25170289 + (576460752303423488*(eta__2 - 4066832143866835/4503599627370496)^2)/2374649627355687 - 1);
W=rand(II+1,JJ+1,3,M);
q=rand(M,1);
Wxy2 = sym('Wxy2',[1 M]);
Wxy3 = sym('Wxy3',[1 M]);
Wxy2(1:M) = sym('0');
Wxy3(1:M) = sym('0');
for r=1:M
for i=1:II+1
for j=1:JJ+1
Wxy2(r) = W(i, j, 2, r)*legendreP(i, zeta__2)*legendreP(j, eta__2)*q(r, 1) + Wxy2(r);
Wxy3(r) = W(i, j, 3, r)*legendreP(i, zeta__2)*legendreP(j, eta__2)*q(r, 1) + Wxy3(r);
end
end
end
Qn__2 = [vpaintegral(vpaintegral(Wxy2(r)*heaviside(-Hvs2)*(abs(Wxy2-Wxy3)'),zeta__2,-1,1),eta__2,-1,1)];

Accepted Answer

Walter Roberson
Walter Roberson on 23 Sep 2022
Wxy2(r) = W(i, j, 2, r)*legendreP(i, zeta__2)*legendreP(j, eta__2)*q(r, 1) + Wxy2(r);
Wxy3(r) = W(i, j, 3, r)*legendreP(i, zeta__2)*legendreP(j, eta__2)*q(r, 1) + Wxy3(r);
You are calculating the exact same legendre on both lines. Calculate the product into a temporary variable and use the temporary variable in both lines.
  38 Comments

Sign in to comment.

More Answers (1)

James Tursa
James Tursa on 23 Sep 2022
The Symbolic Toolbox is going to be slower than IEEE floating point ... that's just something you have to accept. And if you need to have those integer numbers represented exactly you should probably create them as symbolic integers, not double precision integers. E.g., your values with more than 15 decimal digits seem to be exactly representable:
sprintf('%f',5070602400912917605986812821504)
ans = '5070602400912917605986812821504.000000'
sprintf('%f',81129638414606681695789005144064)
ans = '81129638414606681695789005144064.000000'
sprintf('%f',324518553658426726783156020576256)
ans = '324518553658426726783156020576256.000000'
So I am guessing these came from some calculation that ensures this, but to guarantee this in general you would need to do something like this instead:
sym('5070602400912917605986812821504')
ans = 
5070602400912917605986812821504
  1 Comment
Mehdi
Mehdi on 23 Sep 2022
I think the problem is on loops rather than those symbolic numeric problems.
syms eta__2 zeta__2
II=10;JJ=11;M=3;
Hvs2 = sym('5070602400912917605986812821504')*(zeta__2);
W=rand(II+1,JJ+1,3,M);
q=rand(M,1);
Wxy2 = sym('Wxy2',[1 M]);
Wxy3 = sym('Wxy3',[1 M]);
Wxy2(1:M) = sym('0');
Wxy3(1:M) = sym('0');
for r=1:M
for i=1:II+1
for j=1:JJ+1
Wxy2(r) = W(i, j, 2, r)*legendreP(i, zeta__2)*legendreP(j, eta__2)*q(r, 1) + Wxy2(r);
Wxy3(r) = W(i, j, 3, r)*legendreP(i, zeta__2)*legendreP(j, eta__2)*q(r, 1) + Wxy3(r);
end
end
end
Qn__2 = [vpaintegral(vpaintegral(Wxy2(r)*heaviside(-Hvs2)*(abs(Wxy2-Wxy3)'),zeta__2,-1,1),eta__2,-1,1)];

Sign in to comment.

Products


Release

R2021b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!