How can i plot surfaces for a recursive function?

I have written a recursive function from a given recurrence formula, using the form "function z=F(n,x,y)" which evaluates the function at (x,y).

Now i want to plot the surfaces of Fn(x,y) for a range of n (0:5) and -1≤ x,y≤ 1. To get them all in the same figure i know i need to use subplot, but i am struggling with how to use surf for a recursive function?

Any help much appreciated :-)

[Tommy]

Could you provide your code, or a simple example?

Perhaps to generate the plot for a given n, you could use

ax = subplot(..)
fsurf(ax, @(x,y) F(n,x,y), [-1 1])

[Natalie Jalili]

function z=F(n,x,y)
[X,Y]=meshgrid(x,y);
if n==0
    z=ones(size(X));
elseif n==1
    z=(0.5*X)-Y.^2;
else
    z=2*n*X*Y*F((n-1),x,y)-((2*n)+1)*F((n-2),x,y);
    z=z./(2*n.^2);
end
end

This is my code for the recursive function

[John D'Errico]

Actually, No. That is not the code for a recursive function when written properly.

You did not use the .* operator where needed. I've re-written the one line of importance here:

function z=F(n,x,y)
[X,Y]=meshgrid(x,y);
if n==0
    z=ones(size(X));
elseif n==1
    z=(0.5*X)-Y.^2;
else
    z=2*n*X.*Y.*F((n-1),x,y)-((2*n)+1)*F((n-2),x,y);
    z=z./(2*n.^2);
end
end

Note the difference in this line:

    z=2*n*X.*Y.*F((n-1),x,y)-((2*n)+1)*F((n-2),x,y);

Use of the * there instead makes it a MATRIX multiplication, which while it would not create an error, is not valid for what you are doing.

Regardless, you can use this function in the answer I wrote.

Now, is your implementation of this function a good one? NO! Why not? There are several reasons why not.

1. The recursion recomputes the meshgrid EVERY time. That is wasteful, since X and Y never change.
2. This is a recursive function, much like the Fibonacci sequence. Like that sequence, it is a bad idea to use direct recursion. That is because you keep on recomputing the same values. That is, consider how many times F(x,y,0) will be called? In fact, the number of recursive calls grows exponentially with n.

Instead, you want to compute a 3 term recurrence very differently. A simple loop is often best. Again, the Fibonacci sequence is a splendid example of just that, where well coded evaluation of the sequence requires O(n) calls to evaluate F(n).

I'll show in my answer how to EFFICIENTLY compute this recursive function.

Answer 1

[Tommy]

Try this:

for n = 0:5
    ax = subplot(2,3,n+1);
    fsurf(ax, @(x,y) F(n,x,y), [-1 1])
end

[Natalie Jalili]

Ah yes this is just what i was after!! Thank you. However an orange warning message appears in the command window: function behaves unexpectedly on array inputs

[Tommy]

fsurf does the work of creating the grid for x and y, so there is no need to call meshgrid inside your function. All your function has to do is take in two inputs x and y and return z, whose size should equal the size of x (and y), where the values contained within z tell fsurf the values of your function at each (x,y) pair contained within x and y:

for n = 0:5
    ax = subplot(2,3,n+1);
    fsurf(ax, @(x,y) F2(n,x,y), [-1 1])
end

function z=F2(n,x,y)
    if n==0
        z=ones(size(x));
    elseif n==1
        z=(0.5*x)-y.^2;
    else
        z=2*n*x.*y.*F2((n-1),x,y)-((2*n)+1)*F2((n-2),x,y);
        z=z./(2*n^2);
    end
end

Note that the function you posted above isn't necessarily wrong, it just doesn't match what fsurf expects. If I call

F(1,[1 2],[1 2])

for your function, the result is a 2x2 matrix, but fsurf expects a 1x2 matrix because both inputs are 1x2 matrices. You can still accomplish what you're trying to do with your original function, using surf instead of fsurf:

x = linspace(-1,1,30);
for n = 0:5
    ax = subplot(2,3,n+1);
    surf(ax, x, x, F(n, x, x));
end

function z=F(n,x,y)
    [X,Y]=meshgrid(x,y);
    if n==0
        z=ones(size(X));
    elseif n==1
        z=(0.5*X)-Y.^2;
    else
        z=2*n*X.*Y.*F((n-1),x,y)-((2*n)+1)*F((n-2),x,y);
        %      ^  ^ note, added some dots because X and Y and F(...) are not scalars
        z=z./(2*n.^2);
    end
end

(edit) yes admittedly it is redundant to call F(n,x,y) without using the results of F(n-1,x,y) and F(n-2,x,y) which you've already computed.

[Natalie Jalili]

Thank you for your help, i have learnt a lot and i really appreciate it! Especially adding the dots- i had overlooked that completely. Its all in the detail! :-)

Answer 2

[John D'Errico]

Being recursive is not relevant to anything material. As long as you can evaluate it, do so. For example, you might create an array of gridded surfaces. I'm not feeling very creative, and you have not supplied any function.

f(x,y,n) = f(x,y,n-1) + f(x,y,n-2)

So it will act very Fibonacci-like, but the initial conditions will vary as a function of x and y, and of course, why make things integers? ;-)

f(x,y,0) = 1

f(x,y,1) = x + y

The nice thing is, I have no idea how this will all pan out as n increases, but we can just have fun doing it. Because of the recursive nature of this beast, it will be best if we build the arrays up with n. Now for some MATLAB code.

[x,y] = meshgrid(linspace(-1,1,11));
Nmax = 5;
F = zeros([size(x),Nmax+1]);
F(:,:,1) = 2; % remember that MATLAB has an index origin of 1, NOT 0.
F(:,:,2) = x + y;
for i = 3:Nmax+1
  F(:,:,i) =  F(:,:,i-1) + F(:,:,i-2);
end

% now plot the surfaces
for iplot = 1:Nmax + 1
  subplot(2,3,iplot)
  surf(x,y,F(:,:,iplot))
  xlabel x
  ylabel y
  title(['n = ',num2str(iplot-1)])
end

As I said, the recursive nature of the beast is not that relevant. I evaluated the function, stored the result in an array, then plotted it. Could I have done things differently? Of course. But lacking any information about your function, I chose a scheme that works nicely for the function in question.

Edit:

Now that we know the actual recursive function in question, here is how I would write it. Since the goal is still to plot all surfaces for n= 0:5, I'll compute them all in one loop. There is no need to even write the function separately.

Nmax = 5;
x = linspace(-1,1,21);
y = linspace(-1,1,21);
[X,Y]=meshgrid(x,y);
z = zeros([size(X),Nmax+1]);
z(:,:,1) = 1; % remember that MATLAB has an index origin of 1, NOT 0.
z(:,:,2) = (0.5*X)-Y.^2;
for n = 3:Nmax+1
  z(:,:,n) = (2*n*X.*Y.*z(:,:,n-1)-((2*n)+1)*z(:,:,n-2))./(2*n.^2);
end

% now plot the surfaces
for iplot = 1:Nmax + 1
  subplot(2,3,iplot)
  surf(X,Y,z(:,:,iplot))
  xlabel x
  ylabel y
  title(['n = ',num2str(iplot-1)])
end

Could I have written it differently? Again, yes.

[Natalie Jalili]

Thank you so much for this John! I am very new to Matlab and have much to learn, but your explanantions helped me enormously. I had to write the function and plots seperately as part of a homework question, but it does make more sense to compute in one loop. Thanks again!! :-)

Answer 3

[Bobby Fischer]

Hi, I have written this code. Is this what you meant?

function []=bobby_98(n)
lon=20;
x=linspace(-1,1,lon);
y=x;
z=zeros(lon);
close(figure(1))
figure(1)
hold on
for k=1:n
    x=x+1.3;
    for i=1:length(x)
        for j=1:length(y)
            z(i,j)=-5*sin((x(i)-1.3*k)^2+y(j)^2);
        end
    end
    surf(x,y,z)
end
axis equal
end

function []=bobby_99(n)
lon=20;
x=linspace(-1,1,lon);
y=x;
z=zeros(lon);
close(figure(1))
figure(1)
for k=1:n
    subplot(1,n,k)
    for i=1:length(x)
        for j=1:length(y)
            z(i,j)=-5*sin((x(i))^2+y(j)^2);
        end
    end
    surf(x,y,z)
    title(num2str(k))
    axis equal
end
end

Best regards.