Problem with mesh grid and intersection of two surfaces

Hello everyone,

I’m trying to plot the intersection points between two second-order polynomial equations. These equations depend on two variables, a and b, which I define as vectors.

To evaluate the solutions, I use meshgrid to create all combinations of a and b. However, this leads to an issue: for a single root (solution), the same value of a appears repeatedly in the result. This causes my plot to show vertical lines instead of individual points, since multiple identical a values are associated with different b values (or vice versa).

Ideally, for each unique pair (a, b), I want to obtain a unique solution. But using meshgrid evaluates all possible combinations, including repeated values for a or b.

I’d appreciate any advice on how to better approach this problem, particularly in terms of how to structure the evaluation and plotting to avoid this repetition and get cleaner, more accurate plots.

Thanks in advance!

[Matt J]

We will need you to run your code for us (here in the online environment) to illustrate the problem. Off the top there are a few things I don't understand.

(1) If a and b are the dependent variables, it's hard to see why repititions in either a or b would cause vertical lines. They should be on the xy axes.

(2) Ideally, for each unique pair (a, b), I want to obtain a unique solution. Why would that be expected? If these are 2nd order equations, there could be as many as two real solutions. Or none.

[Zac Lee]

clear all; clc; close all
pas = 0.001;
N11list = 0:pas:0.01;
N31list = 0:pas:0.01;
[N11vec, N31vec] = meshgrid(N11list, N31list);
N = [N11vec(:), N31vec(:)];
racine = [];
N11sol = [];
N31sol = [];
% Stockage des valeurs déjà rencontrées (pour chaque variable séparément)
valeurs_n11 = [];
valeurs_n31 = [];
valeurs_racine = [];
tol = 1e-8;
for i = 1:length(N)
    n11 = N(i, 1);
    n31 = N(i, 2);
    
    %% Polynôme 1
    A1 = n11^2;
    B1 = 2*n11 * n31^2;
    C1 = (5 * n11^2 + 3 * n31);
    %% Polynôme 2
    A2 = n31^2;
    B2 = n31 + n11 * n31^3;
    C2 = n31 * n11^2 + n31 * n11^3+1;
    %% Intersection
    A = A1 - A2;
    B = B1 - B2;
    C = C1 - C2;
    if A == 0 && B == 0
        continue
    end
    roots_poly = roots([A, B, C]);
    for k = 1:length(roots_poly)
        x = roots_poly(k);
        if isreal(x)
            % Vérifier que la racine n’a pas déjà été vue
            if any(abs(valeurs_racine - x) < tol)
                continue
            end
            % Tout est unique : on stocke
            racine(end+1) = x;
            N11sol(end+1) = n11;
            N31sol(end+1) = n31;
            valeurs_racine(end+1) = x;
            valeurs_n11(end+1) = n11;
            valeurs_n31(end+1) = n31;
        end
    end
end
% Affichage
scatter3(racine,N11sol,N31sol,'.b');
xlabel('ROOT');
ylabel('N11');
zlabel('N31');

[Zac Lee]

Dear @Matt J, thanks for your answer.

I would like to keep only real solutions. If there are two real solutions, they shouyld be kept. My ptoblem, is that, i'm having multiple times the same value ffor a and b for the same root.

[Catalytic]

My ptoblem, is that, i'm having multiple times the same value ffor a and b for the same root.

Are you trying to develop a mapping from (a,b) to the roots, or the other way around? If you are mapping from (a,b) to the roots then it shouldn't matter if there are multiple (a,b) that map to the same root.

If you are mapping from the roots back to (a,b) then you first have to map the roots to a set of polynomial coefficients and then map the coefficients to (a,b). Bear in mind though that polynomial coefficients are unique only up to a scale factor - so you need to decide on some way to remove that ambiguity, before you can map the coefficients to a unique (a,b).

[Zac Lee]

Dear @Catalytic, thank you verymuch for your answer. Effectively, I want to develop a mapping from (a, b) to the roots. The main problem is that I have several repeated instances of the same (a, b) pairs, which leads to the issue shown in the picture. I think I’ll need to apply some kind of filtering to ensure that each (a, b) pair appears only once. That way, I can be sure to avoid this problem — or at least reduce it significantly. I think i have to approach the problem in another way

Answer 1

[Catalytic]

Effectively, I want to develop a mapping from (a, b) to the roots

If so, then @Matt J already told you what you are doing wrong. (a,b), or I guess (N11,N13), should be on the xy axis -

pas = 0.001;
[N11, N31] = meshgrid(0:pas:0.01); 

%% Polynôme 1
A1 = N11.^2;
B1 = 2 .* N11 .* N31.^2;
C1 = (5 .* N11.^2 + 3 .* N31);
%% Polynôme 2
A2 = N31.^2;
B2 = N31 + N11 .* N31.^3;
C2 = N31 .* N11.^2 + N31 .* N11.^3 + 1;
%% Intersection
A = A1 - A2;
B = B1 - B2;  %normalize leading coefficient
C = C1 - C2;  %normalize leading coefficient


D=B.^2-4.*A.*C; %Discriminant
R1=(-B-sqrt(D))./2./A; %First root
R2=(-B+sqrt(D))./2./A; %Second root

k=D>0 & (A|B);
R1=R1(k); R2=R2(k); N11=N11(k); N31=N31(k);


scatter3(N11,N31,R1,'.b'); hold on
scatter3(N11,N31,R2,'.r'); hold off
xlabel('N11');
ylabel('N12');
zlabel('ROOT');
legend ROOT1 ROOT2

●

Answer 2

[Matt J]

Perhaps as follows:

tol = 1e-8;

pas = 0.001;
[N11, N31] = meshgrid(0:pas:0.01); N11=N11(:); N31=N31(:);

%% Polynôme 1
A1 = N11.^2;
B1 = 2 .* N11 .* N31.^2;
C1 = (5 .* N11.^2 + 3 .* N31);
%% Polynôme 2
A2 = N31.^2;
B2 = N31 + N11 .* N31.^3;
C2 = N31 .* N11.^2 + N31 .* N11.^3 + 1;
%% Intersection
A = A1 - A2;
B = (B1 - B2)./A;  %normalize leading coefficient
C = (C1 - C2)./A;  %normalize leading coefficient


D=B.^2-4*C; %Discriminant
R1=(-B-sqrt(D))/2; %First root
R2=(-B+sqrt(D))/2; %Second root

%Post-process: discard redundant or illegal solutions
T=[R1,R2,N11,N31];

 T( any(~isfinite(T),2) | D<0 ,  : )=[]; %discard non-finite or complex solutions
 
 [~,Iu]=uniquetol(T(:,1:2),tol,'ByRows',true); %discard repeat solutions (within tolerance)

 T=num2cell(T(Iu,:),1);

[R1,R2,N11,N13]=deal(T{:});

figure
scatter3(R1,R2,N11,'.b'); 
xlabel('ROOT #1');
ylabel('ROOT #2');
zlabel('N11');

figure
scatter3(R1,R2,N13,'.b'); 
xlabel('ROOT #1');
ylabel('ROOT #2');
zlabel('N13');

[Zac Lee]

Dear @Matt J and @Catalytic, it seems working npow. Thank you verymuch for your help.