Having problem with a infinite double sum

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.

[Lucas Resende]

I just need help putting that double sum into code, I can figure out from that on

[Torsten]

Just out of curiosity: Which PDE does the above sum solve ?

[Lucas Resende]

Momentum in the x axis of a thin square plate.

Answer 1

[Walter Roberson]

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)

[Lucas Resende]

Thanks so much! I'll try it like that!

[Torsten]

@Walter Roberson

Does your code try to find a symbolic expression for the double sum ?

Does vpasum do all the numerics for the user - thus decide up to which indices to sum such that an acceptable precision is reached ?

[Walter Roberson]

@Torsten

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.

vpasum() does somehow figure out when to stop; see https://www.mathworks.com/help/symbolic/vpasum.html#mw_bd9ec10e-e8f5-4483-a536-d7e7a3a5ec26

[Torsten]

Thank you for the explanation.

Answer 2

[Torsten]

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