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

6 views (last 30 days)
Mehdi on 23 Sep 2022
Edited: Mehdi on 29 Sep 2022
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)];

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.
Torsten on 29 Sep 2022
Ok, then take the way that best fits your needs.

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
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)];

### Categories

Find more on Programming in Help Center and File Exchange

R2021b

### Community Treasure Hunt

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

Start Hunting!