Find values of a function in the same region and save them into different arrays for each region

Hello,

I have created a meshgrid divided in 13x17 regions and I have a function that travels through this region. This function that iterates in a for loop is an array and what I want to do is create different arrays that groups the values of this function if they are in the same region. Sorry if I'm not explaining very well, english is not my first language and I its a complex work for me. I attached an image that shows this function in red travelling through the grid.

I'm having a big trouble trying to separate those values from the function and grouping them in other arrays if they meet a condition (the condition is to be in the same region of the grid). Anyone who can help me?

Attached code:

%% this script evaluates the response of two second order functions and identifies the plant using a look up table
%% INIT
clear variables;
close all;
clc;

%% variable initialization
% time
time=150;
deltaT=0.01;
t = (0:deltaT:time);

%% Reference
A=2;
w=0.8;
r=A*sin(w*t);
rd=gradient(r,t);

[a,b,x1init,x2init,d1,d2,w1,w2,w12,w22,d1sigma,d1sigma2,d2tau,d2tau2,mup,x1p,x2p,fintp,x1d,x2d,fintd,x1cp,x2cp,fintcp] = initVar(t);

% initial conditions
x1p(1)=x1init; %2
x2p(1)=x2init;

x1d(1)=x1init;
x2d(1)=x2init;

Kp=3;
Kd=0.8;
%% SOLVE SYSTEM AND CREATE THE PLANE

% fd function for desired plant PD
fd=@(pd1,pd2)(a*pd1+b*pd2);

% fp function for van der pol system plant (P plant)
fp=@(fp1,fp2)(mup.*(1-(fp1.^2)).*fp2-fp1);

% define working area to plot the plane 
[X,Y]=meshgrid(d1,d2);

%values of both functions x2prima at vertexes
PDLUT= fd(X, Y)

PDLUT = 13×17

4.5000 4.3500 4.2000 4.0500 3.9000 3.7500 3.6000 3.4500 3.3000 3.1500 3.0000 2.8500 2.7000 2.5500 2.4000 2.2500 2.1000 3.9500 3.8000 3.6500 3.5000 3.3500 3.2000 3.0500 2.9000 2.7500 2.6000 2.4500 2.3000 2.1500 2.0000 1.8500 1.7000 1.5500 3.4000 3.2500 3.1000 2.9500 2.8000 2.6500 2.5000 2.3500 2.2000 2.0500 1.9000 1.7500 1.6000 1.4500 1.3000 1.1500 1.0000 2.8500 2.7000 2.5500 2.4000 2.2500 2.1000 1.9500 1.8000 1.6500 1.5000 1.3500 1.2000 1.0500 0.9000 0.7500 0.6000 0.4500 2.3000 2.1500 2.0000 1.8500 1.7000 1.5500 1.4000 1.2500 1.1000 0.9500 0.8000 0.6500 0.5000 0.3500 0.2000 0.0500 -0.1000 1.7500 1.6000 1.4500 1.3000 1.1500 1.0000 0.8500 0.7000 0.5500 0.4000 0.2500 0.1000 -0.0500 -0.2000 -0.3500 -0.5000 -0.6500 1.2000 1.0500 0.9000 0.7500 0.6000 0.4500 0.3000 0.1500 0 -0.1500 -0.3000 -0.4500 -0.6000 -0.7500 -0.9000 -1.0500 -1.2000 0.6500 0.5000 0.3500 0.2000 0.0500 -0.1000 -0.2500 -0.4000 -0.5500 -0.7000 -0.8500 -1.0000 -1.1500 -1.3000 -1.4500 -1.6000 -1.7500 0.1000 -0.0500 -0.2000 -0.3500 -0.5000 -0.6500 -0.8000 -0.9500 -1.1000 -1.2500 -1.4000 -1.5500 -1.7000 -1.8500 -2.0000 -2.1500 -2.3000 -0.4500 -0.6000 -0.7500 -0.9000 -1.0500 -1.2000 -1.3500 -1.5000 -1.6500 -1.8000 -1.9500 -2.1000 -2.2500 -2.4000 -2.5500 -2.7000 -2.8500

PLUT= fp(X, Y)

PLUT = 13×17

49.0000 37.2500 27.0000 18.2500 11.0000 5.2500 1.0000 -1.7500 -3.0000 -2.7500 -1.0000 2.2500 7.0000 13.2500 21.0000 30.2500 41.0000 41.5000 31.6250 23.0000 15.6250 9.5000 4.6250 1.0000 -1.3750 -2.5000 -2.3750 -1.0000 1.6250 5.5000 10.6250 17.0000 24.6250 33.5000 34.0000 26.0000 19.0000 13.0000 8.0000 4.0000 1.0000 -1.0000 -2.0000 -2.0000 -1.0000 1.0000 4.0000 8.0000 13.0000 19.0000 26.0000 26.5000 20.3750 15.0000 10.3750 6.5000 3.3750 1.0000 -0.6250 -1.5000 -1.6250 -1.0000 0.3750 2.5000 5.3750 9.0000 13.3750 18.5000 19.0000 14.7500 11.0000 7.7500 5.0000 2.7500 1.0000 -0.2500 -1.0000 -1.2500 -1.0000 -0.2500 1.0000 2.7500 5.0000 7.7500 11.0000 11.5000 9.1250 7.0000 5.1250 3.5000 2.1250 1.0000 0.1250 -0.5000 -0.8750 -1.0000 -0.8750 -0.5000 0.1250 1.0000 2.1250 3.5000 4.0000 3.5000 3.0000 2.5000 2.0000 1.5000 1.0000 0.5000 0 -0.5000 -1.0000 -1.5000 -2.0000 -2.5000 -3.0000 -3.5000 -4.0000 -3.5000 -2.1250 -1.0000 -0.1250 0.5000 0.8750 1.0000 0.8750 0.5000 -0.1250 -1.0000 -2.1250 -3.5000 -5.1250 -7.0000 -9.1250 -11.5000 -11.0000 -7.7500 -5.0000 -2.7500 -1.0000 0.2500 1.0000 1.2500 1.0000 0.2500 -1.0000 -2.7500 -5.0000 -7.7500 -11.0000 -14.7500 -19.0000 -18.5000 -13.3750 -9.0000 -5.3750 -2.5000 -0.3750 1.0000 1.6250 1.5000 0.6250 -1.0000 -3.3750 -6.5000 -10.3750 -15.0000 -20.3750 -26.5000

%Controller plant LUT is calculated directly
CLUT=PDLUT-PLUT;

% dimensions of the PB model
X1=X(1,:);
X2=Y(:,1);

%% LOOP
% interpolates the system using the working area X Y, LUT table(f at
% vertexese) and the pos and speed at instant (x1p, x2p, x1d, x2d) (emulates Ode45)
nIterations=length(t);

x1dError=zeros(1,length(t));
x2dError=zeros(1,length(t));
x1pError=zeros(1,length(t));
x2pError=zeros(1,length(t));
u=zeros(1,length(t));
ud=zeros(1,length(t));
up=zeros(1,length(t));

F=zeros(length(t),1);
Xmat=zeros(length(t),4);
phi=zeros(4,1);

tic
for ind2=2:nIterations
    %% P PLANT + CP PLANT
    % introduce reference signal
    r1=r(1,ind2-1);
    rd1=rd(1,ind2-1);

    %  if the function is between conditions, then the controller C starts
    %  working.
     % using the C controller plant values and P plant observations (x1p, x2p), interpolates the value of U
    u(1,ind2)=interp2(X,Y,CLUT,x1p(1,ind2-1),x2p(1,ind2-1));

    % calculate response of van der pol system:
    % 1.- first interpolates fp funtion value using the last observations
    % of the P plant (x1pant, x2pant)(initial values in the first loop x1,x2=2)
    % 2.- using u and the new f, calculates the next observations of the P
    % plant.
    [x1p(1,ind2),x2p(1,ind2),fintp(1,ind2)]=ffunc(X,Y,PLUT,x1p(1,ind2-1),x2p(1,ind2-1),deltaT,u(1,ind2));

    % Position and velocity error calculated 
    x1pError(1,ind2)=r1-x1p(1,ind2-1);
    x2pError(1,ind2)=rd1-x2p(1,ind2-1);

    % PID output "U"
    ufb(1,ind2)= Kp*x1pError(1,ind2)+Kd*x2pError(1,ind2);

    [x1p(1,ind2),x2p(1,ind2),fintp(1,ind2)]=fCP(X,Y,PLUT,x1p(1,ind2-1),x2p(1,ind2-1),deltaT,u(1,ind2),ufb(1,ind2-1));
    
    %% PD PLANT
    % Position and velocity error calculated 
    x1dError(1,ind2)=r1-x1d(1,ind2-1);
    x2dError(1,ind2)=rd1-x2d(1,ind2-1);

    % PID output "U"
    ud(1,ind2)= Kp*x1dError(1,ind2)+Kd*x2dError(1,ind2);
    % calculate response of desired plant:
    % 1.- first interpolates fd funtion value using the last observations
    % of the PD plant (x1dant, x2dant)(initial values in the first loop x1,x2=2)
    % 2.- using u and the new f, calculates the next observations of the DP
    % plant. 
    [x1d(1,ind2),x2d(1,ind2),fintd(1,ind2)]=ffunc(X,Y,PDLUT,x1d(1,ind2-1),x2d(1,ind2-1),deltaT,ud(1,ind2));

    %% CP PLANT
    
    x1cp=x1p;
    x2cp=x2p;
    %calc x1cp and x2cp

    % get the f function of the CP close loop plant 
    fintcp(1,ind2)= fintp(1,ind2)+ u(1,ind2);

    %% Calculation of weights alpha and beta
    %find in which part of the table corresponds x1
    X1find=find(X1<x1cp(1,ind2));
    X1ind=X1find(end);
    d1sigma(1,ind2)=X1(X1ind);
    d1sigma2(1,ind2)=X1(X1ind+1);
    
    %find in which part of the table corresponds x2
    X2find=find(X2<x2cp(1,ind2));
    X2ind=X2find(end);
    d2tau(1,ind2)=X2(X2ind);
    d2tau2(1,ind2)=X2(X2ind+1);

    %aisle values of fintcp from a certain region (HERE THE ISSUE)
    % I want to save all those values of fintcp that are in the same region
    % into different arrays for each region

    %calculate alpha and beta 
    w1(1,ind2)=(x1cp(1,ind2)-d1sigma(1,ind2))/(d1sigma2(1,ind2)-d1sigma(1,ind2));
    w2(1,ind2)=(x2cp(1,ind2)-d2tau(1,ind2))/(d2tau2(1,ind2)-d2tau(1,ind2));
    w12(1,ind2)=1-w1(1,ind2);
    w22(1,ind2)=1-w2(1,ind2);

end 
toc

Elapsed time is 7.469420 seconds.

    
phi=inv(Xmat'*Xmat)*Xmat'*F;

Warning: Matrix is singular to working precision.

   

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%FUNCTION

function [x1act,x2act,f]= ffunc(X,Y,LUT,x1ant,x2ant,deltaT,u)

% calculate f(x1,x2) with the vertexes (LUT, X and Y) and instant
% pos/speed(x1,x2)
f=interp2(X,Y,LUT,x1ant,x2ant);

x1act= deltaT*x2ant+x1ant;
x2act=(f+u)*deltaT+x2ant;

end

%%%%%%%%%%%

function [x1act,x2act,f]= fCP(X,Y,LUT,x1ant,x2ant,deltaT,u,r)

% calculate f(x1,x2) with the vertexes (LUT, X and Y) and instant
% pos/speed(x1,x2)
f=interp2(X,Y,LUT,x1ant,x2ant);

x1act= deltaT*x2ant+x1ant;
x2act=(f+u+r)*deltaT+x2ant;

end



%%%%%%%%%%
%INITIALIZATION 
function [a,b,x1init,x2init,d1,d2,w1,w2,w12,w22,d1sigma,d1sigma2,d2tau,d2tau2,mup,x1p,x2p,fintp,x1d,x2d,fintd,x1cp,x2cp,fintcp] = initVar(t)
% a refers to realtion of electromotive force and the rotor inertia a=K/J
% b refers to relation of viscuous friction parameter and the rotor b=-b/J
% inertia
a=-0.3;
b=-1.1;

%initial conditions
x1init= 2;
x2init= 2;

% vertex
d1=(-4:0.5:4);
d2=(-3:0.5:3);

w1=zeros(1,length(t));
w2=zeros(1,length(t));
w12=zeros(1,length(t));
w22=zeros(1,length(t));
d1sigma=zeros(1,length(t));
d1sigma2=zeros(1,length(t));
d2tau=zeros(1,length(t));
d2tau2=zeros(1,length(t));

%van der pol inestability and damping force parameter
mup=1;

%van der pol position and speed vector init
x1p=nan(1,length(t));
x2p=nan(1,length(t));

x1d=nan(1,length(t));
x2d=nan(1,length(t));

x1cp=nan(1,length(t));
x2cp=nan(1,length(t));

fintd=nan(1,length(t));
fintp=nan(1,length(t));
fintcp=nan(1,length(t));
end 

%%%%%%%%

[Matt Butts]

I am not sure I completely understand your question, but it seems like histcounts2 might be able to help. This could tell you which bin each point belongs to in your grid.

[~,~,~,binX,binY] = histcounts2(xFuncVals,yFuncVals,xEdges,yEdges);
% Points in 1,1 of mesh
xFuncVals(binX == 1 & binY == 1)
yFuncVals(binX == 1 & binY == 1)