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.
3 个评论
Lucas Resende
2022-4-4
Torsten
2022-4-4
Just out of curiosity: Which PDE does the above sum solve ?
Lucas Resende
2022-4-5
采纳的回答
Walter Roberson
2022-4-5
编辑：Walter Roberson
2022-4-5
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)
4 个评论
Lucas Resende
2022-4-5
Torsten
2022-4-5
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
2022-4-5
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
2022-4-5
Thank you for the explanation.
更多回答（1 个）
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
类别
在 帮助中心 和 File Exchange 中查找有关 Calculus 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
