quadruple summation using function

If I use function to calculate the double summation for , the code is

n = 100;

a = @(p) sin(p);

sum(sum(a(0:n)' * a(0:n)))

However, if I want to calculate the quadruple summation , how to modify the code?

[Dyuman Joshi]

Similar to my answer in a previous question of yours, use ndgrid

%Random value for N for example
N = 10;

kvec = 0:1:N; %1 as increment is not necessary

[p,q,r,s] = ndgrid(kvec);
arr = sin(p).*sin(q).*sin(r).*sin(s);

s = sum(arr,'all')

s = 3.9659

[Matt J]

But note that this approach will be quite slow:

%Random value for N for example
N = 100;

kvec = 0:1:N; %1 as increment is not necessary

tic
 [p,q,r,s] = ndgrid(kvec);
 arr = sin(p).*sin(q).*sin(r).*sin(s);
 s1 = sum(arr,'all');
toc

Elapsed time is 3.951835 seconds.

tic
 s2=sum(sin(kvec))^4;
toc

Elapsed time is 0.002075 seconds.

s1,s2

s1 = 2.6155e-04

s2 = 2.6155e-04

[Matt J]

You should be exploiting the separability of the summation,

and likewise for the quadruple sum.

[Dyuman Joshi]

Good point, @Matt J

Answer 1

[Matt J]

n = 100;
tic;
sum(sin(0:n))^4;
toc

Elapsed time is 0.002355 seconds.

Answer 2

[Steven Lord]

You could use implicit expansion to avoid having to create quite so many large arrays. It's not as fast as the case that exploits the separability, but it is significantly faster than the original general approach.

%Random value for N for example
N = 100;

kvec = 0:1:N; %1 as increment is not necessary

tic
 [p,q,r,s] = ndgrid(kvec);
 arr = sin(p).*sin(q).*sin(r).*sin(s);
 s1 = sum(arr,'all');
toc

Elapsed time is 3.802691 seconds.

tic
 s2=sum(sin(kvec))^4;
toc

Elapsed time is 0.002488 seconds.

% Use implicit expansion
tic
n = numel(kvec);
sinK = sin(kvec);
p = reshape(sinK, n, 1, 1, 1); % unnecessary in this case, but useful for generality
q = reshape(sinK, 1, n, 1, 1);
r = reshape(sinK, 1, 1, n, 1);
s = reshape(sinK, 1, 1, 1, n);
s3 = sum(p.*q.*r.*s, 'all');
toc

Elapsed time is 0.332431 seconds.

The trailing 1's in the reshape calls aren't really necessary, but they do make the pattern of sizes quite easy to see.

format longg
results = [s1, s1-s2, s1-s3; s2-s1, s2, s2-s3; s3-s1, s3-s2, s3]

results = 3×3

1.0e+00 * 0.000261548679431158 -2.30172814402047e-13 0 2.30172814402047e-13 0.00026154867966133 2.30172814402047e-13 0 -2.30172814402047e-13 0.000261548679431158

Those diagonal elements are in pretty good agreement, and all the off-diagonal elements (the differences between the approaches) are all quite small in magnitude,

[Matt J]

You can simplify this by downloading ndgridVecs from the File Exchange,

https://www.mathworks.com/matlabcentral/fileexchange/74956-ndgridvecs

tic
 [p,q,r,s] = ndgridVecs(sin(0:100));
 result = sum( p.*q.*r.*s ,  'all');
toc

Elapsed time is 0.268219 seconds.