Help with a double series?

I am trying to write function for the equation, can someone tell me what i am doing wrong

clear,clc
a=input('Please enter the dimension, a\n');
b=input('Please enter the dimension, b\n');
h=input('Please enter the thickness, h\n');
E=input('Please enter the Modulus of Eleasticity, E\n');
v=input('Please enter the Poisson Ratio, v\n');

P=input('Please enter the Load, P\n');
xo=input('Please enter the Position, xo\n');
yo=input('Please enter the Position, yo\n');
rec_plate_pt(a,b,P,xo,yo,x,y);

D=(E*h^3)/(12*(1-v^2))
c=(4*P)/(pi^4*a*b*D)

function w = rec_plate_pt(a,b,P,xo,yo,x,y)

s=0;
for m=1:Inf
    for n=1:Inf
s=c*(sin((m*pi*xo)/a)*sin((n*pi*yo)/b)/((m^2/a^2+n^2/b^2).^2)*(sin((m*pi*x)/a))*(sin(n*pi*y)/b))
        
        
format long
        fprintf('%d',w);
    end
end
end

[Walter Roberson]

Those are infinite series. You should use symsum()

Answer 1

[Jyothis Gireesh]

As Walter mentioned in the previous comment the double summation may be implemented using a combination of symbolic variables and the “symsum()” in MATLAB instead of defining a function.

One possible implementation is given below

syms Wc(x,y) m n; 
Wc(x,y) = c*symsum(symsum(sin(m*pi*x0/a)*sin(n*pi*y0/b)*sin(m*pi*x/a)*sin(n*pi*y/b)/(m^2/a^2 + n^2/b^2)^2,n,1,Inf),m,1,Inf); 

Answer 2

[Johannes Ebert]

Hello together,

I have a question: if I apply the code above with my parameter set, the result is the following:

My Input: Wc(120,200)

The result:

(7312256873090709*symsum(symsum(((exp(-(pi*n*10i)/31)*1i)/2 - (exp((pi*n*10i)/31)*1i)/2)^2/(m^2/722500 + n^2/384400)^2, n, 1, Inf)*((exp(-(pi*m*12i)/85)*1i)/2 - (exp((pi*m*12i)/85)*1i)/2)^2, m, 1, Inf))/39614081257132168796771975168

Is there a chance to get the concret value for the used parameter set?

Thanks und kind regards,

Johannes

[Walter Roberson]

syms m n real
Wc = (sym('7312256873090709')*symsum(symsum(((exp(-(pi*n*10i)/31)*1i)/2 - (exp((pi*n*10i)/31)*1i)/2)^2/(m^2/722500 + n^2/384400)^2, n, 1, Inf)*((exp(-(pi*m*12i)/85)*1i)/2 - (exp((pi*m*12i)/85)*1i)/2)^2, m, 1, Inf))/sym('39614081257132168796771975168')

Wc = 

vpa(Wc)

ans = 

S = rewrite(Wc,'sincos')

S = 

simplify(S,'steps',50)

ans = 

[Walter Roberson]

Hmmm, this does not seem to help.

syms m n integer
part1 = ((exp(-(pi*m*12i)/85)*1i)/2 - (exp((pi*m*12i)/85)*1i)/2)^2;
part2 = symsum(((exp(-(pi*n*10i)/31)*1i)/2 - (exp((pi*n*10i)/31)*1i)/2)^2/(m^2/722500 + n^2/384400)^2, n, 1, Inf)

part2 = 

part2 = simplify(rewrite(part2, 'sin'), 'steps', 50)

part2 = 

part3 = symsum( part2 .* part1, m, 1, Inf)

part3 = 

part3 = simplify(rewrite(part3, 'sin'), 'steps', 50)

part3 = 

Wc = (sym('7312256873090709') .* part3)/sym('39614081257132168796771975168')

Wc = 

[Walter Roberson]

Maple is able to find an explicit formula for part2 (that is, the inner symsum() ) in terms of a fair number of psi() calls and some trig terms. But it struggles to find anything definite further than that or to simplify the mess.