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

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

Answer 1

[Walter Roberson]

                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.

[Mehdi]

what do you mean by temporary variable?

[Walter Roberson]

A temporary variable is a variable that is used for only a short period of time, and whose value does not need to be retained after its immediate use.

You also do not need to recalculate the legendre expressions for all of the different r values.

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');
    legi = zeros(1, II+1);
    legj = zeros(1, JJ+1);
    for i = 1 : II+1; legi = legendreP(i, eta__2); end
    for j = 1 : JJ+1; legj = legendreP(j, eta__2); end
    for r=1:M
        for i=1:II+1
            for j=1:JJ+1                
                legij = legi(i) * legj(j);
                Wxy2(r) = W(i, j, 2, r) * legij * q(r, 1) + Wxy2(r);
                Wxy3(r) = W(i, j, 3, r) * legij * 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)];

[Mehdi]

What about using parallel programming to speed up the calculation?(thread, task)

[Walter Roberson]

You could reduce the legendre calculation further by calculating once to the maximum of II+1 and JJ+1 putting the result into a single vector, and then indexing the vector inside the loops, legij=leg(i)*leg(j)

[Mehdi]

I did so, but it has not considerable effect on calculation speed. I am thinking about using parallel programming to speed up the calculation (thread, task). Let me know if somebody have opinion in this case.

[Torsten]

Let me know if somebody have opinion in this case.

If speed matters so much, why do you use symbolic instead of numerical computation ? This is the wrong approach.

[Mehdi]

can you explain a little more?

[Torsten]

I gave you numerical code:

https://de.mathworks.com/matlabcentral/answers/1807210-any-comment-to-speed-up-the-sum-calculation-of-a-function-which-includes-legendre-polynomials?s_tid=srchtitle

Why didn't you continue this way ?

[Mehdi]

In this way how to deal with compute\ing double integral Qn__2 ?

[Torsten]

You write a function (not a function handle) in which you calculate the value of the integrand given coordinates x and y.

This function is then called from integral2.

[Mehdi]

I faced problem, can you explain it on a small example.

Note that in my orginal problem Hvs2 is changing with time although here in my example I set it not to change for simplicity.

[Torsten]

I dont understand why you write Wxy2(r) in the expression

Qn__2 = [vpaintegral(vpaintegral(Wxy2(r)*heaviside(-Hvs2)*(abs(Wxy2-Wxy3)'),zeta__2,-1,1),eta__2,-1,1)];

r is the last value of the preceeding loop over r, thus M, in this case.

So you want to compute

Qn__2 = vpaintegral(vpaintegral(Wxy2(M)*heaviside(-Hvs2)*(abs(Wxy2-Wxy3)'),zeta__2,-1,1),eta__2,-1,1);

here ?

And you want to get an Mx1 vector as result ?

Or should the double integral be taken for r=1,...,M and you expect an MxM matrix as result ?

[Mehdi]

yes Qn__2 will give an Mx1 vector as result.

Qn__2 = [vpaintegral(vpaintegral(Wxy2(r)*heaviside(-Hvs2)*(abs(Wxy2-Wxy3)'),zeta__2,-1,1),eta__2,-1,1)];

is correct. This is why I computed Wxy2(r) in previous loops.

[Torsten]

But r = M in the integral.

[Mehdi]

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

[Torsten]

Not possible since Wxy2 and Wxy3 are not yet complete to be used within the integral.

Further Qn__2 is an Mx1 vector ; you can't save it in a scalar Qn__2(r).

Or do you mean

Qn__2 (r)= [vpaintegral(vpaintegral(Wxy2(r)*heaviside(-Hvs2)*abs(Wxy2(r)-Wxy3(r)),zeta__2,-1,1),eta__2,-1,1)];

[Mehdi]

Sorry! bad correction. Eddited again. Now it seems fine.

[Torsten]

for r=1:M
    wxy2=Wxy2(s)+wxy2;
    wxy3=Wxy3(s)+wxy3;   
end

A loop over r and Wxy2 and Wxy3 indexed by s which is undefined ?

Take the time and properly write down what you want to calculate (maybe in a mathematical notation, if code is too difficult). Then come back here again.

[Mehdi]

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]

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

[Mehdi]

As I pointed out in my previous post, Hvs2 is changing (is computed from another complex procedures) , although in my example for simplicty I set it a constant function. So lets modify your suggested code in a way that it works for new Hvs2s automatically.

Maybe there must be a way to set it outside of the function and import it to the function by an argument! something like bellow

Hvs2=@(x,y)((5070602400912917605986812821504.*(x +...

value = arrayfun(@(r)integral2(@(x,y)fun(x,y,r,II,JJ,M,W,q),-1,1,-1,1),1:M);

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

But faced error!

[Torsten]

It gives values in the order of 1e161. I think you should first check whether it's calculated correctly.

And it's not constant - it depends on x and y.

On what else should it depend ?

[Mehdi]

Ok, I will check its correctness ASAP.

Yes Hvs2 is not constant - it depends on x and y, but each time new function of x and y.

[Torsten]

Yes Hvs2 is not constant - it depends on x and y, but each time new function of x and y.

This won't work since integration is deterministic, not stochastic. The function to be integrated must remain unchanged during the integration process.

Or what do you expect as result of the integration if the integrand changes with each call to it ? A random variable ?

[Mehdi]

there must be a way to determine Hvs2 outside of the function and pass it to the function by an argument! something like bellow

Hvs2=@(x,y) sin(x*y)

value = arrayfun(@(r)integral2(@(x,y)fun(x,y,r,II,JJ,M,W,q),-1,1,-1,1),1:M);

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

Is not it possible? if not, how to deal with it?

[Torsten]

Yes, that's possible, but - according to what you wrote - I suspected that you wanted to change Hvs2 during the integration.

Hvs2=@(x,y) sin(x.*y)
value = arrayfun(@(r)integral2(@(x,y)fun(x,y,r,II,JJ,M,W,q,Hvs2),-1,1,-1,1),1:M);

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

[Mehdi]

I think there is a problem in your suggested codes, since getting any result even for small values of II, JJ, M.

clc
clear
II=1;JJ=1;M=2;
W=rand(II+1,JJ+1,3,M);
q=rand(M,1);
Hvs2=@(x,y)sin(x.*y);
value = arrayfun(@(r)integral2(@(x,y)fun(x,y,r,II,JJ,M,W,q,Hvs2),-1,1,-1,1),1:M);
function values = fun(x,y,r,II,JJ,M,W,q,Hvs2_fun)
Hvs2 = Hvs2_fun(x,y);
   for ii = 1:II+1
    part1(ii,:,:) = legendreP(ii, x);
  end
  for jj= 1:JJ+1
    part2(jj,:,:) = legendreP(jj, y);
  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
        Wxy2(s,:,:) = W(ii, jj, 2, s)*part1(ii,:,:).*part2(jj,:,:)*q(s, 1) + Wxy2(s,:,:);
        Wxy3(s,:,:) = W(ii, jj, 3, s)*part1(ii,:,:).*part2(jj,:,:)*q(s, 1) + Wxy3(s,:,:);
      end
    end
  end
  values = squeeze(Wxy2(r,:,:)).*heaviside(-Hvs2).*squeeze(abs(sum(Wxy2-Wxy3,1)));
end

[Torsten]

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
  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

[Mehdi]

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!

[Torsten]

The last code is for plotting the function to see where there might be problems in integration. To plot the function over [-1 1]x[-1 1], the inputs x and y to "fun" are matrices and the result "values" is a matrix of the same size as x and y.

The integration code before gives an 1xM vector coming from the line

value = arrayfun(@(r)integral2(@(x,y)fun(x,y,r,II,JJ,M,W,q,Hvs2),-1,1,-1,1),1:M);

"fun" is called by "integral2" also with matrices for x and y as inputs, and the output variable "values" is also a matrix of the same size as x and y. But the result from "integral2" is a scalar for each value of r. Since it is called for r=1:M, the result is an 1xM vector.

[Mehdi]

Running has not finished after waiting long time. !

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

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]

Thanks for your comments. As you suggested the only option is to use double trapz to approximate integrals.

We switched from symbolic to numerical computation to speed up the calculations, but no success at all.

[Torsten]

The problem is not the speed to get the function values in "fun", but to get a result from integral2.

Even this simple example needs more than 50 seconds to run.

II = 8;
JJ = 5;
tic
value = integral2(@(x,y)fun(x,y,II,JJ),-1,1,-1,1)

value = 1.9077e-13

toc

Elapsed time is 54.946624 seconds.

function values = fun(x,y,II,JJ)
  n1 = size(x,1);
  n2 = size(x,2);
  n = n1*n2;
  x = reshape(x,n,1);
  y = reshape(y,n,1);

  part1 = zeros(n,II+1);
  part2 = zeros(n,JJ+1);
  part = zeros(n,II+1,JJ+1);
  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
  values = zeros(n,1);
  for ii=1:II+1
    for jj=1:JJ+1
        values = values + part(:,ii,jj);
    end
  end
  values = reshape(values,[n1 n2]);
end

[Mehdi]

the symbolic one takes around 3 seconds to run.

[Torsten]

That's surprising. Could you include the code for testing ?

[Mehdi]

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]

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

Answer 2

[James Tursa]

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]

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