Any comment to speed up the speed of caculation of symbolic loops having Legendre polynomials?
Show older comments
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
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
Mehdi
on 23 Sep 2022
Walter Roberson
on 23 Sep 2022
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
on 24 Sep 2022
Walter Roberson
on 24 Sep 2022
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
on 24 Sep 2022
Torsten
on 24 Sep 2022
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
on 24 Sep 2022
I gave you numerical code:
Why didn't you continue this way ?
Mehdi
on 24 Sep 2022
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 ?
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
on 24 Sep 2022
Torsten
on 24 Sep 2022
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.
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
on 24 Sep 2022
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
on 24 Sep 2022
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 ?
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
on 25 Sep 2022
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
on 28 Sep 2022
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
on 28 Sep 2022
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
on 29 Sep 2022
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)
toc
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
on 29 Sep 2022
Torsten
on 29 Sep 2022
That's surprising. Could you include the code for testing ?
Torsten
on 29 Sep 2022
Ok, then take the way that best fits your needs.
More Answers (1)
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)
sprintf('%f',81129638414606681695789005144064)
sprintf('%f',324518553658426726783156020576256)
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')
1 Comment
Mehdi
on 23 Sep 2022
Categories
Find more on Programming in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
