MATLAB中的双重傅里叶级数

Question

我可以为下面的 f 和 g 函数在 MATLAB 中绘制网格图： 我已经尝试了以下代码来绘制 f 和 g 的图像：

%% plot f
[x,y] = meshgrid(linspace(-pi,pi,50));
f = x.*y ;
subplot(1,2,1)
mesh(f)
title('f')

%% plot g
syms  m n
A = 4*(-1)^(m+n)*(sin(m*x)*sin(n*y))/(m*n);
g = symsum(symsum(A,n,1,inf),m,1,inf);
subplot(1,2,2)
mesh(g)
title('g')

< p > mesh 的结果如下：

绘制 f 的部分没有任何错误。另一个绘制 g 的部分在图中没有显示。我该如何绘制 g？

- asys

---

这不是代码错误。抱歉，我已经修改了我的问题。 - asys

为什么你在 g 中使用符号数学？不要这样做，因为你无法绘制符号。 - Adriaan

那么如果没有symsum，我该如何定义求和？ - asys

你试图去无限远，那将需要无限长的时间。编写一个 for 循环，并在特定的 n 和 m 顺序下退出，看看结果有多接近。 - Adriaan

有人能帮我写for循环吗？我已经很累了，对这段代码完全感到困惑，抱歉。 - asys

不，因为Stack Overflow并不是一个代码编写服务。只需要循环遍历m和n两个维度，并在最后对两个维度的和求和即可。 - Adriaan

Answer 1

如果你要处理符号数学，尤其是涉及周期函数和具有不连续性的函数时，最好先熟悉假设。您可能还想使用 fmesh （或旧版本中的 ezmesh ）来绘制符号表达式的网格：

syms m n x y
assume(in(m,'integer') & m>=1);
assume(in(n,'integer') & n>=1);
assume(x>-pi & x<pi);
assume(y>-pi & y<pi);
A = 4*(-1)^(m+n)*(sin(m*x)*sin(n*y))/(m*n);

g = symsum(symsum(A,n,1,Inf),m,1,Inf);
fmesh(g,5*[-pi pi -pi pi],'MeshDensity',1e2); % or ezmesh(g,5*[-pi pi -pi pi]);

这将创建一个如下所示的图表: 另一个选项是使用subs和double对g进行数值计算，然后使用mesh进行绘图:

[X,Y] = meshgrid(linspace(-5*pi,5*pi,100));
g2 = real(double(subs(g,{x,y},{X,Y})));
mesh(g2);

或者使用matlabFunction创建一个数值函数：

g2 = matlabFunction(g);
[X,Y] = meshgrid(linspace(-5*pi,5*pi,100));
mesh(real(g2(X,Y)));

在这两种情况下，必须使用real来裁剪由于数值不精确而产生的无关紧要的虚部。

- horchler

Answer 2

如果你想使用显式的 for 循环来实现，可以按照以下方式：

XY = 2*linspace(-pi,pi,50);
N = 100;
M = 100;
G = zeros(numel(XY));
tmp = zeros(M*N,1);
cx = 1;
for x = XY
    cy = 1;
    for y = XY
        c = 1;
        for m = 1:M
            for n = 1:N
                tmp(c) = 4*(-1)^(m+n)*(sin(m*x)*sin(n*y))/(m*n);
                c = c+1;
            end
        end
        G(cx,cy) = sum(tmp);
        cy = cy+1;
    end
    cx = cx+1;
end
mesh(G)
title('g')

这是一种更紧凑的方式，应该会更快一些：

XY = 2*linspace(-pi,pi,50);
N = 100;
M = 100;
G = zeros(numel(XY));
cx = 1;
for x = XY
    cy = 1;
    for y = XY
        gfun = @(m,n) 4.*(-1).^(m+n).*(sin(m.*x).*sin(n.*y))./(m.*n);
        tmp = bsxfun(gfun,(1:M).',1:N);
        G(cx,cy) = sum(tmp(:));
        cy = cy+1;
    end
    cx = cx+1;
end
mesh(G)

并且结果是：

- EBH