How to plot an Implicit function with certain conditions

im trying to plot in MATLAB the red implicit function given these 3 conditions, in DESMOS it's super easy but unfortunately i have no idea how to pull it off in MATLAB, i've searched the internet but my luck fell short, would appreciate the help.

Answer 1

[Pratyush Swain]

Hi Shai,

To plot this implicit functions you need to leverage the fimplicit function in MATLAB.Here is an demonstrated example on how we can use this function for your usecase:

warning('off','all');
% Plotting the function using fimplicit
% NOTE: Here I have arbitrarily chosen input values for x and y
% x range --> 1<=x<=2.5, y range --> 0<=y<=1.5 
fimplicit(@fun, [1, 2.5, 0, 1.5]);
title('Implicit Function with Conditions');
xlabel('x');
ylabel('y');

% Callback for Impicit to calculate the values, Here x and y are retreived as vector of input values 
function values = fun(x, y)
    % Initialize output array to NaN %
    values = NaN(size(x));
    
    % Define the condition for y to avoid division by zero %
    valid = y ~= 0;
    
    % Calculate the expression (2x/y)*(1-1/2y) under valid condition %
    % NOTE: Here the ./ operater ensures division occurs on all valid
    % elements of vector x
    expression = (2*x(valid)./y(valid)).*(1-1./(2*y(valid)));
    
    % Apply all conditions and obtain a final logical vector
    condition = x > 1.166 & expression > 0 & expression < pi/2 & y < 1;
    
    % Apply the implicit function where conditions are met
    values(condition) = (2*x(condition)./y(condition)) - tan(expression(condition));
end

I have referred to this thread to form a solution for this usecase: https://www.mathworks.com/matlabcentral/answers/1756760-how-to-plot-implicit-function-with-conditions

For more information on fimplicit function, please refer to https://www.mathworks.com/help/matlab/ref/fimplicit.html

[Shai Zipori]

unfortunately, the plot that's being displayed isn't appropriate (when increasing the x range beyond 2.5, it displays crushed line) could you help me with getting the points (coordinates) of given x as input (lets say 1.166:0.01:5.1) and y as output that fullfills these constraints?

0<((2*x)/(y))*(1-((1)/(2*y)))<(pi/2)

y<1

and suits this (same) equation?

((2*x)/(y))=tan(((2*x)/(y))*(1-((1)/(2*y))))

unfortunately, i just can't seem to grasp how to use the logical vectors.

Thank you so much for your help!

[Mathieu NOE]

you can improve the result by increasing the 'MeshDensity' factor :

h = fimplicit(@fun, [1, 5, 0.5, 1.1], 'MeshDensity',300);

then you get rid of the waves and NaNs

% Plotting the function using fimplicit
% NOTE: Here I have arbitrarily chosen input values for x and y
% x range --> 1<=x<=5, y range --> 0.5<=y<=1.1
h = fimplicit(@fun, [1, 5, 0.5, 1.1], 'MeshDensity',300);

Warning: Function behaves unexpectedly on array inputs. To improve performance, properly vectorize your function to return an output with the same size and shape as the input arguments.

title('Implicit Function with Conditions');
xlabel('x');
ylabel('y');

% if you need to access to the x & y data
x = h.XData;
y = h.YData;

figure
plot(x,y);
title('Implicit Function with Conditions');
xlabel('x');
ylabel('y');

% Callback for Impicit to calculate the values, Here x and y are retreived as vector of input values 
function values = fun(x, y)
    % Initialize output array to NaN %
    values = NaN(size(x));
    
    % Define the condition for y to avoid division by zero %
    valid = y ~= 0;
    
    % Calculate the expression (2x/y)*(1-1/2y) under valid condition %
    % NOTE: Here the ./ operater ensures division occurs on all valid
    % elements of vector x
    expression = (2*x(valid)./y(valid)).*(1-1./(2*y(valid)));
    
    % Apply all conditions and obtain a final logical vector
    condition = x > 1.166 & expression > 0 & expression < pi/2 & y < 1;
    
    % Apply the implicit function where conditions are met
    values(condition) = (2*x(condition)./y(condition)) - tan(expression(condition));
end

Answer 2

[Mathieu NOE]

hello @Shai Zipori

I have to say I don't do much with implicit function problems , tried using fsolve and fminbnd but was lacking how to implement the conditions

now, a poor's man solution is to create a x, y grid , evaluate your function, and find the x,y points where you function is minimal

% solving implicit function 
% (2*x./y) - tan((2*x./y).*(1-1./(2*y))) = 0;

x = linspace(1.167,5,1000); % create x array with condition x>1.166
y = linspace(0.5,1,1000);

% create a X Y meshgrid and evaluate function 
[X,Y] = meshgrid(x,y);
C = (2*X./Y).*(1-1./(2*Y));
fun = (2*X./Y) - tan(C); % this is your implicit function 

% apply condition (C>0 & C<pi/2)
ind = (C>0 & C<pi/2);
Z = NaN(size(fun));
Z(ind) = fun(ind); % keep only valid fun values according to condition (C>0 & C<pi/2)

% plot function Z=f(x,y) to show minimum line (is what we are looking for)
figure
h = imagesc(x,y,log(abs(Z)));
colorbar('vert')
set(gca,'YDir','normal')
set(h, 'AlphaData', 1-isnan(Z))

% find x,y of minimum line 
for ci = 1:numel(x)
    zz = Z(:,ci);
    % find y coordinate to get minimum z value 
    [val,ind] = min(abs(zz));
    if ~isempty(ind)
        xc(ci) = x(ci);
        yc(ci) = y(ind);
    end
end

figure
plot(xc,yc);

[Mathieu NOE]

the method above will give you a correct plot if you have a refined grid resolution . if you reduce the grid size , the final curve will exibit some staircase behaviour

of course , one solution is to smooth the coarse curve

% solving implicit function 
% (2*x./y) - tan((2*x./y).*(1-1./(2*y))) = 0;

x = linspace(1.167,5,300); % create x array with condition x>1.166
y = linspace(0.5,1,300);

% create a X Y meshgrid and evaluate function 
[X,Y] = meshgrid(x,y);
C = (2*X./Y).*(1-1./(2*Y));
fun = (2*X./Y) - tan(C); % this is your implicit function 

% apply condition (C>0 & C<pi/2)
ind = (C>0 & C<pi/2);
Z = NaN(size(fun));
Z(ind) = fun(ind); % keep only valid fun values according to condition (C>0 & C<pi/2)

% plot function Z=f(x,y) to show minimum line (is what we are looking for)
figure
h = imagesc(x,y,log(abs(Z)));
colorbar('vert')
set(gca,'YDir','normal')
set(h, 'AlphaData', 1-isnan(Z))

% find x,y of minimum line 
for ci = 1:numel(x)
    % find y coordinate to get minimum z value 
    [val,ind] = min(abs(Z(:,ci)));
    if ~isempty(ind)
        xc(ci) = x(ci);
        yc(ci) = y(ind);
    end
end


yf = smoothdata(yc,'loess',25);

figure
plot(xc,yc,'b',xc,yf,'r');
legend('coarse','smoothed');

[Mathieu NOE]

or you could do this :

% solving implicit function 
% (2*x./y) - tan((2*x./y).*(1-1./(2*y))) = 0;

x = linspace(1.167,5,300); % create x array with condition x>1.166
y = linspace(0.5,1,300);

% create a X Y meshgrid and evaluate function 
[X,Y] = meshgrid(x,y);
C = (2*X./Y).*(1-1./(2*Y));
fun = (2*X./Y) - tan(C); % this is your implicit function 

% apply condition (C>0 & C<pi/2)
ind = (C>0 & C<pi/2);
Z = NaN(size(fun));
Z(ind) = fun(ind); % keep only valid fun values according to condition (C>0 & C<pi/2)

% plot function Z=f(x,y) to show minimum line (is what we are looking for)
figure
h = imagesc(x,y,log(abs(Z)));
colorbar('vert')
set(gca,'YDir','normal')
set(h, 'AlphaData', 1-isnan(Z))

% find x,y of minimum line 
k1 = 0;
k2 = 0;

for ci = 1:numel(x)
    
    % find x, y coordinates of minimum Z value 
    % 1/ coarse method (relies on grid resolution)
    [val,ind] = min(abs(Z(:,ci)));
    if ~isempty(ind)
        k1 = k1 + 1;
        xc(k1) = x(ci);
        yc(k1) = y(ind);
    end
    
   % 2/ refined method with linear interpolation
    [ZxP,ZxN] = find_zc(y,Z(:,ci),0);
    if ~isempty(ZxN)
        k2 = k2 + 1;
        xf(k2) = x(ci);
        yf(k2) = ZxN;
    end
end

figure
plot(xc,yc,'b',xf,yf,'r');
legend('coarse','refined');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [ZxP,ZxN] = find_zc(x,y,threshold)

% put data in rows
x = x(:);
y = y(:);

% positive slope "zero" crossing detection, using linear interpolation
y = y - threshold;
zci = @(data) find(diff(sign(data))>0);  %define function: returns indices of +ZCs
ix=zci(y);                      %find indices of + zero crossings of x
ZeroX = @(x0,y0,x1,y1) x0 - (y0.*(x0 - x1))./(y0 - y1); % Interpolated x value for Zero-Crossing 
ZxP = ZeroX(x(ix),y(ix),x(ix+1),y(ix+1));

% negative slope "zero" crossing detection, using linear interpolation
zci = @(data) find(diff(sign(data))<0);  %define function: returns indices of +ZCs
ix=zci(y);                      %find indices of + zero crossings of x
ZeroX = @(x0,y0,x1,y1) x0 - (y0.*(x0 - x1))./(y0 - y1); % Interpolated x value for Zero-Crossing 
ZxN = ZeroX(x(ix),y(ix),x(ix+1),y(ix+1));

end

[Shai Zipori]

Thank you Mathieu NOE, worked like a charm.

[Mathieu NOE]

my pleasure !