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

Question

Mehdi 2022-9-23

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1811020-any-comment-to-speed-up-the-speed-of-caculation-of-symbolic-loops-having-legendre-polynomials

编辑： Mehdi 2022-9-29

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Walter Roberson 2022-9-23

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1811020-any-comment-to-speed-up-the-speed-of-caculation-of-symbolic-loops-having-legendre-polynomials#answer_1059585

在 MATLAB Online 中打开

                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 个评论
显示 36更早的评论隐藏 36更早的评论

Mehdi 2022-9-24

编辑：Mehdi 2022-9-24

在 MATLAB Online 中打开

sorry, now corrected.

syms eta__2 zeta__2
II=10;JJ=11;M=3;
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 s=1:M
    for i=1:II+1
        for j=1:JJ+1
            Wxy2(s) = W(i, j, 2, s)*legendreP(i, zeta__2)*legendreP(j, eta__2)*q(s, 1) + Wxy2(s);
            Wxy3(s) = W(i, j, 3, s)*legendreP(i, zeta__2)*legendreP(j, eta__2)*q(s, 1) + Wxy3(s);
       end
   end
end
for r=1:M
Qn__2(r,1) = [vpaintegral(vpaintegral(Wxy2(r)*heaviside(-Hvs2)*(abs(sum(Wxy2)-sum(Wxy3))),zeta__2,-1,1),eta__2,-1,1)];
end

Mehdi 2022-9-24

编辑：Mehdi 2022-9-24

在 MATLAB Online 中打开

I deleted them and used sum. now it is fine. Excuse again.

syms eta__2 zeta__2
II=10;JJ=11;M=3;
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 s=1:M
    for i=1:II+1
        for j=1:JJ+1
            Wxy2(s) = W(i, j, 2, s)*legendreP(i, zeta__2)*legendreP(j, eta__2)*q(s, 1) + Wxy2(s);
            Wxy3(s) = W(i, j, 3, s)*legendreP(i, zeta__2)*legendreP(j, eta__2)*q(s, 1) + Wxy3(s);
       end
   end
end
for r=1:M
Qn__2(r,1) = [vpaintegral(vpaintegral(Wxy2(r)*heaviside(-Hvs2)*(abs(sum(Wxy2)-sum(Wxy3))),zeta__2,-1,1),eta__2,-1,1)];
end

Torsten 2022-9-24

编辑：Torsten 2022-9-24

在 MATLAB Online 中打开

Check whether it's correctly implemented.

I don't know the runtime.

II=10;JJ=11;M=3;
W=rand(II+1,JJ+1,3,M);
q=rand(M,1);
Hvs2_fun = @(x,y)((5070602400912917605986812821504.*(x + 2251799813683713./2251799813685248).^2)./2356225 + (9007199254740992.*(y + 2935286035937695./18014398509481984).^2)./196937227765191 - 1).*((81129638414606681695789005144064.*(x + 9007199254732683./9007199254740992).^2)./69039481 + (576460752303423488.*(y + 3261970163074917./4503599627370496).^2)./6904142590940591 - 1).*((324518553658426726783156020576256.*(x + 140737488355209./140737488355328).^2)./231983361 + (144115188075855872.*(y - 262292457514301./562949953421312).^2)./2637878570603985 - 1).*((144115188075855872.*(x + 4028041154330599./4503599627370496).^2)./424643881623313 + y.^2 - 1).*((20282409603651670423947251286016.*(x - 4503599627213111./4503599627370496).^2)./24770038225 + (288230376151711744.*(y - 7530397878711147./9007199254740992).^2)./5204731445635785 - 1).*((324518553658426726783156020576256.*(x + 4503599627365785./4503599627370496).^2)./355058649 + (36893488147419103232.*(y + 4434826747744735/4503599627370496).^2)./8603290501959015 - 1).*((4611686018427387904.*(y + 2213733699584161./2251799813685248).^2)./1317884237102575 + (324518553658426726783156020576256.*(x - 4503599627284663./4503599627370496).^2)./117876175561 - 1).*((81129638414606681695789005144064.*(x + 9007199254735975./9007199254740992).^2)./25170289 + (576460752303423488.*(y - 4066832143866835./4503599627370496).^2)./2374649627355687 - 1);
value = arrayfun(@(r)integral2(@(x,y)fun(x,y,r,II,JJ,M,W,q,Hvs2_fun),-1,1,-1,1),1:M);
function values = fun(x,y,r,II,JJ,M,W,q,Hvs2_fun)
  Hvs2 = Hvs2_fun(x,y);
  part1 = zeros(II+1,size(x,1),size(x,2));
  part2 = zeros(JJ+1,size(x,1),size(x,2));
  part = zeros(II+1,JJ+1,size(x,1),size(x,2));
  for ii = 1:II+1
    part1(ii,:,:) = legendreP(ii, x);
  end
  for jj= 1:JJ+1
    part2(jj,:,:) = legendreP(jj, y);
  end
  for ii = 1:II+1
    for jj = 1:JJ+1
      part(ii,jj,:,:) = part1(ii,:,:).*part2(jj,:,:);
    end
  end
  Wxy2 = zeros(M,size(x,1),size(x,2));
  Wxy3 = zeros(size(Wxy2));
  for s=1:M
    for ii=1:II+1
      for jj=1:JJ+1
        partij = reshape(part(ii,jj,:,:),1,size(x,1),size(x,2)); 
        Wxy2(s,:,:) = W(ii, jj, 2, s)*partij*q(s, 1) + Wxy2(s,:,:);
        Wxy3(s,:,:) = W(ii, jj, 3, s)*partij*q(s, 1) + Wxy3(s,:,:);
      end
    end
  end
  values = squeeze(Wxy2(r,:,:)).*heaviside(-Hvs2).*squeeze(abs(sum(Wxy2-Wxy3,1)));
end

Torsten 2022-9-25

编辑：Torsten 2022-9-25

在 MATLAB Online 中打开

I don't see an error in the code (that I tried to optimize further a bit).

Note that heaviside introduces a sharp discontinuity. You must be aware that this might turn integration extremely difficult or even impossible.

II=10;JJ=11;M=3;

W=rand(II+1,JJ+1,3,M);

q=rand(M,1);

Hvs2_fun = @(x,y)sin(x.*y);

r = 2;

x = -1:0.1:1;

y = -1:0.1:1;

[X,Y] = meshgrid(x,y);

Z = fun(X,Y,r,II,JJ,M,W,q,Hvs2_fun);

surf(X,Y,Z)

function values = fun(x,y,r,II,JJ,M,W,q,Hvs2_fun)

Hvs2 = Hvs2_fun(x,y);

part1 = zeros(II+1,size(x,1),size(x,2));

part2 = zeros(JJ+1,size(x,1),size(x,2));

part = zeros(II+1,JJ+1,size(x,1),size(x,2));

for ii = 1:II+1

part1(ii,:,:) = legendreP(ii, x);

end

for jj= 1:JJ+1

part2(jj,:,:) = legendreP(jj, y);

end

for ii = 1:II+1

for jj = 1:JJ+1

part(ii,jj,:,:) = part1(ii,:,:).*part2(jj,:,:);

end

Wxy2 = zeros(M,size(x,1),size(x,2));

Wxy3 = zeros(size(Wxy2));

for s=1:M

for ii=1:II+1

for jj=1:JJ+1

partij = reshape(part(ii,jj,:,:),1,size(x,1),size(x,2));

Wxy2(s,:,:) = W(ii, jj, 2, s)*partij*q(s, 1) + Wxy2(s,:,:);

Wxy3(s,:,:) = W(ii, jj, 3, s)*partij*q(s, 1) + Wxy3(s,:,:);

end

values = squeeze(Wxy2(r,:,:)).*heaviside(-Hvs2).*squeeze(abs(sum(Wxy2-Wxy3,1)));

end

Torsten 2022-9-28

编辑：Torsten 2022-9-29

在 MATLAB Online 中打开

Why do you stress your comments by exclamation marks ? Shall I feel responsible for that your assignment does not make progress ?

I don't know the exact reason why integral2 needs so long for computation. I told you that "heaviside" and "abs" are poison for every integrator.

If you are satisfied with results without error estimates, choose

x = -1:0.1:1;

y = -1:0.1:1;

evaluate the function on this mesh by calling "fun" and use trapz twice on the matrix of function values returned. This will give you an approximation of the double integral.

I think there is still a mistake in your codes, since I am looking for a Mx1 vector as result, but your codes result in 21*21!

I could reproduce this behaviour. Apparently, the "squeeze" commands to calculate "values" change dimensions not as wanted if x and y are vector inputs. Changing the last line from

values = squeeze(Wxy2(r,:,:)).*heaviside(-Hvs2).*squeeze(abs(sum(Wxy2-Wxy3,1)));

to

values = reshape(Wxy2(r,:,:),[size(x,1),size(x,2)]).*heaviside(-Hvs2).*reshape(abs(sum(Wxy2-Wxy3,1)),[size(x,1) size(x,2)]);

should remove this fault.

Mehdi 2022-9-29

编辑：Mehdi 2022-9-29

在 MATLAB Online 中打开

clear
syms eta__2 zeta__2
II=1;JJ=1;M=2;
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 s=1:M
    for i=1:II+1
        for j=1:JJ+1
            Wxy2(s) = W(i, j, 2, s)*legendreP(i, zeta__2)*legendreP(j, eta__2)*q(s, 1) + Wxy2(s);
            Wxy3(s) = W(i, j, 3, s)*legendreP(i, zeta__2)*legendreP(j, eta__2)*q(s, 1) + Wxy3(s);
       end
   end
end
for r=1:M
Qn__2(r,1) = [vpaintegral(vpaintegral(Wxy2(r)*heaviside(-Hvs2)*(abs(sum(Wxy2)-sum(Wxy3))),zeta__2,-1,1),eta__2,-1,1)];
end

Torsten 2022-9-29

Ok, then take the way that best fits your needs.

请先登录，再进行评论。

Answer 2

James Tursa 2022-9-23

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1811020-any-comment-to-speed-up-the-speed-of-caculation-of-symbolic-loops-having-legendre-polynomials#answer_1059570

在 MATLAB Online 中打开

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 个评论
显示 -1更早的评论隐藏 -1更早的评论

Mehdi 2022-9-23

在 MATLAB Online 中打开

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