# Having problem with a infinite double sum

2 views (last 30 days)
Lucas Resende on 4 Apr 2022
Commented: Torsten on 5 Apr 2022
I'm having a problem trying to code this infine double sum
A, b, x, y and q0 are known values. My problem is specificaly tryind to transform that double sum in a code.
##### 3 CommentsShow 1 older commentHide 1 older comment
Torsten on 4 Apr 2022
Just out of curiosity: Which PDE does the above sum solve ?
Lucas Resende on 5 Apr 2022
Momentum in the x axis of a thin square plate.

Walter Roberson on 5 Apr 2022
Edited: Walter Roberson on 5 Apr 2022
You might be able to simplify the output if you can put constraints on the values.
This takes a while... It might possibly take less time with specific numeric values.
With specific numeric values for everything you could potentially use vpasum()
syms a b M N m n q0 x y nu
Pi = sym(pi);
inner = ((m^2/a^2) + nu*(n^2/b^2)) / (m*n*((m^2/a^2) + (n^2/b^2))^2) * sin(m*Pi*x/a) * sin(n*pi*y/b)
innerMN = subs(inner, {m, n}, {2*M-1, 2*N-1})
M_x = 16*q0/Pi^4 * symsum( symsum(innerMN, N, 1, inf), M, 1, inf)
Walter Roberson on 5 Apr 2022
Yes, my code does try to find the symbolic expression for the double sum. In practice, after thinking a fair while, it returns
(16*q0*symsum((sin((x*pi*(2*M - 1))/a)*symsum((sin((y*pi*(2*N - 1))/b)*((2*M - 1)^2/a^2 + (nu*(2*N - 1)^2)/b^2))/((2*N - 1)*((2*M - 1)^2/a^2 + (2*N - 1)^2/b^2)^2), N, 1, Inf))/(2*M - 1), M, 1, Inf))/pi^4
At the moment I do not know if it would be able to get further if it were given specific numeric values for the constants.
Torsten on 5 Apr 2022
Thank you for the explanation.

Torsten on 5 Apr 2022
Edited: Torsten on 5 Apr 2022
You might want to try a numerical solution:
a = 2.0;
b = 4.0;
nu = 20;
q0 = 1.0;
X = linspace(-2*pi,2*pi,250);
Y = linspace(-2*pi,2*pi,500);
eps = 1e-5; % precision of series evaluation
tic
for i = 1:numel(X)
for j = 1:numel(Y)
Z(j,i) = func(a,b,nu,q0,X(i),Y(j),eps);
end
i
end
toc
[XX,YY] = meshgrid(X,Y) ;
surf(XX,YY,Z)
function fvalue = func(a,b,nu,q0,x,y,eps)
total_sum = 0.0;
diagonal_sum = 1.0;
i = 1;
if abs(nu) >= 1
while abs(diagonal_sum) > eps
J = 1:i;
diagonal_sum = sum((((2*J-1)/a).^2/nu + ((2*(i-J)-1)/b).^2)./...
((2*J-1).*(2*(i-J)-1).*(((2*J-1)/a).^2 +...
((2*(i-J)-1)/b).^2).^2) .*...
sin((2*J-1)*pi*x/a).*sin((2*(i-J)-1)*pi*y/b));
total_sum = total_sum + diagonal_sum;
i = i + 1;
end
total_sum = nu*total_sum;
else
while abs(diagonal_sum) > eps
J = 1:i;
diagonal_sum = sum((((2*J-1)/a).^2 + nu*((2*(i-J)-1)/b).^2)./...
((2*J-1).*(2*(i-J)-1).*(((2*J-1)/a).^2 +...
((2*(i-J)-1)/b).^2).^2) .*...
sin((2*J-1)*pi*x/a).*sin((2*(i-J)-1)*pi*y/b));
total_sum = total_sum + diagonal_sum;
i = i + 1;
end
end
fvalue = 16*q0/pi^4*total_sum;
end