CODEWITHSUNDEEP
matlabmatlab-figuresymbolic-math
asys
asys

Reputation: 729

Double fourier series in MATLAB

I'd like to plot a mesh for both f and g functions below in MATLAB: f and it's double Fourier series g

I'd tried this for f and 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')

The result of mesh is:

the result of mesh

The section plotting f is running without any error. The other section plotting g show nothing in the figure. How can I plot g?

Upvotes: 3

Views: 1266

Answers (2)

horchler
horchler

Reputation: 18484

If you're going to work with symbolic math, it's a good idea to get comfortable with assumptions, especially when dealing with periodic functions and functions with discontinuities. You may also want to use fmesh (or ezmesh in older versions) to plot meshes of symbolic expressions:

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]);

This creates a plot like this: fmesh plot

Another option is to evaluate g numerically using subs and double and then use mesh to plot:

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

or use matlabFunction to create a numeric function:

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

In both these latter cases, real must be used to clip the insignificant imaginary parts due to numerical imprecision.

Upvotes: 3

EBH
EBH

Reputation: 10440

In case you want the explicit way to do it with a for loop, here it is:

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')

And this is a bit more compact way, which should be faster:

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)

and the result:

fourier

Upvotes: 0

Related Questions