How to create a surface plot of a function that uses colon operators?

I have a following function:

function p =  dist(n, theta, M, Ns, Nb, G, kappa)
    Ntheta =  55;
    logterm1 = sum( log(1: n + M - 1) ) - sum( log (1: M - 1)) - sum( log (1:n));
    logterm2 =  n.*(log (Ntheta) - log(1 + Ntheta));
    logterm3 = M*log(1 + Ntheta);
    total = logterm1 + logterm2 - logterm3;
    p = exp(total);
end


function PE = Error(x, y)

    p0 = x;
    p1 = 1-x;
    Ns = 1e-2;
    Nb =10;
    kappa = 1e-2;
    G = 1.01;
    M = 1e3;
    Ntheta0 =  150;
    Ntheta_pi =  33;
    sigma0 = sqrt( Ntheta0.*( Ntheta0 + 1));
    sigma_pi = sqrt( Ntheta_pi.*( Ntheta_pi + 1));
    
    narray = 1:floor(y);

    PX0LessthanTh = 0;
    for na = narray
        pp = dist(na, 0, M, Ns, Nb, G, kappa);
        PX0LessthanTh = PX0LessthanTh + pp;
    end

    PXpiLessthanTh = 0;
    for na = narray
        pp = dist(na, pi, M, Ns, Nb, G, kappa);
        PXpiLessthanTh = PXpiLessthanTh + pp;
    end
    PE = p0.*PX0LessthanTh + p1.*(1 - PXpiLessthanTh);

end

I am looking to create a surface plot of Error with respect to x and y. However, I because I am using a statement narray = 1:floor(y); I am getting an error:

 The following error was reported evaluating the function in FunctionLine
 update: Colon operands must be real scalars.

How can I modify my Error function above so that I can use 3D functions like surf.

Answer 1

[Harshavardhan]

Hi @Rahul Bhadani

The current “Error” function looks to be designed for scalar inputs. When used with “meshgrid” for surface plots, MATLAB passes arrays, leading to issues with expressions like “narray = 1:floor(y);”.

To resolve this, use “arrayfun” to evaluate the function for each (x,y) pair separately as shown below.

%Example x,y arrays
x = linspace(0, 1, 30);
y = linspace(1, 20, 30);
[X, Y] = meshgrid(x, y); %Make the X,Y grid

Z = arrayfun(@(xi, yi) Error(xi, yi), X, Y);% Using “arrayfun” to calculate Z(error)

figure
surf(X, Y, Z);
xlabel('x'); ylabel('y'); zlabel('Error');
title('Surface plot of Error(x, y)');
colorbar;

For more information on “arrayfun” and “meshgrid”, type the following commands in a MATLAB command window

doc meshgrid
doc arrayfun

Hope this helps.

Answer 2

[VBBV]

@Rahul Bhadani It seems you have not defined the y variable. However, If its an external input value to the function Error then you need to ensure that its a scalar. In your case its doesnt seem to be scalar but rather a vector or a matrix. Since you have a for loop inside the function, the colon operator does not parse the vector inputs to a for loop index The solution would be to reshape the matrix and use numel instead of floor as below to get rid of the warning/error message

x = linspace(0, 1, 30);
y = linspace(1, 20, 30);
[X, Y] = meshgrid(x, y); %Make the X,Y grid

PE = Error(X,Y);

narray = 1×900

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

y = 1×900

1.0000 1.6552 2.3103 2.9655 3.6207 4.2759 4.9310 5.5862 6.2414 6.8966 7.5517 8.2069 8.8621 9.5172 10.1724 10.8276 11.4828 12.1379 12.7931 13.4483 14.1034 14.7586 15.4138 16.0690 16.7241 17.3793 18.0345 18.6897 19.3448 20.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

figure
surf(X,Y,PE)

function PE = Error(x, y)

    p0 = x;
    p1 = 1-x;
    Ns = 1e-2;
    Nb =10;
    kappa = 1e-2;
    G = 1.01;
    M = 1e3;
    Ntheta0 =  150;
    Ntheta_pi =  33;
    sigma0 = sqrt( Ntheta0.*( Ntheta0 + 1));
    sigma_pi = sqrt( Ntheta_pi.*( Ntheta_pi + 1));
    
    narray = 1:numel(y)
    y = reshape(y,1,[])  % 
    PX0LessthanTh = 0;
    for na = narray
        pp = dist(y(na), 0, M, Ns, Nb, G, kappa);
                % ^^ 
        PX0LessthanTh = PX0LessthanTh + pp;
    end

    PXpiLessthanTh = 0;
    for na = narray
        pp = dist(y(na), pi, M, Ns, Nb, G, kappa);
                % ^^
        PXpiLessthanTh = PXpiLessthanTh + pp;
    end
    PE = p0.*PX0LessthanTh + p1.*(1 - PXpiLessthanTh);
    
   

end

function p =  dist(n, theta, M, Ns, Nb, G, kappa)
    Ntheta =  55;
    logterm1 = sum( log(1: n + M - 1) ) - sum( log (1: M - 1)) - sum( log (1:n));
    logterm2 =  n.*(log (Ntheta) - log(1 + Ntheta));
    logterm3 = M*log(1 + Ntheta);
    total = logterm1 + logterm2 - logterm3;
    p = exp(total);
end