How to find the intersection of two curves (2 points of intersection)

I have a piece of code that I eventually will turn into a function to find the minimum and maximum velocities an aircraft can fly at, in relation to the power available from the engine, and the power required to fly at those conditions. In the figure from this code, there is a flat line that represents the power available, and a curved line that represents the power required.

Is there a simple way to output the intersections of these lines? Ideally so further down the line this can all be used as a function?

Appreciate all the help I can get!

close all
clear
clc
%Code to determine the maximum and minimum speeds at a given altitude for
%the Cessna 172P

H=0;        %Eventually this will become an input in a function where the 
            % user just call the function and recieve the miniumum and maximum velocities for a given altitude 


%Constants of the aircraft
S = 16.16499443;                %Wing area [m^2]
g = 9.81;                       %acceleration due to gravity [ms^-2]
W = 887.0748299*g;              %Weight of the aircraft [kgf]
C_D0 = 0.026;                   %Zero lift drag coefficient 
K = 0.121;                      %Drag constant

%Range of velocities we wish to investigate [Mach]
M = linspace(0,0.8,100);


[~,rho,~,a]=atmosProperties(H,0,0);     %Density and speed of sound at given altitude [kgm^-3] [ms^-1]
[~,maxPA]=getPAmax(H);                  %Maximum power available at given altitude [W]
V = M.*a;                               %Velocity in [ms^-1]
C_L = (2./rho).*(W./S).*(1./(V.^2));    %Lift coefficient at given altitude
C_D = C_D0 + K .* C_L.^2;               %Drag coefficient at given altitude
P_req = 0.5 .* rho .* V.^3 .* S .* C_D;    %Power required curve

b = [maxPA, maxPA];                     %B and c are used to plot a line that corresponds to the maximum power available
c = [0,272.235190420871];

%Goal of the figure is to find the intersections of the power available curve (flat line) and the power required curve (curved line)
%Legends will be added later on
plot (V, P_req);
hold on
plot (c, b);
xlabel("TAS [ms^-1]");
ylabel("Power reequired (W)")

Answer 1

[Star Strider]

I can’t run your code, since I don’t have the functions.

As a general rule however, it is straightforward to find the intersections between two curves, even if the functions to calculate them do not exist.

Try this with your vectors:

x = linspace(0, 10);                                            % Independent Variable
y1 = 10 - (x-5).^2;                                             % Dependent Variable #1
y2 = 5*sin(2*pi*x/5);                                           % Dependent Variable #1

curvdif = y1 - y2;                                              % Subtract Curves
intidx = find(diff(sign(curvdif)));                             % Approximate Indices Of Zero-Crossings

for k = 1:numel(intidx)
    idxrng = [max(1,intidx(k)-1):min(numel(x),intidx(k)+1)];    % Index Range For Interpolation
    xv(k) = interp1(curvdif(idxrng), x(idxrng), 0);             % Interpolate To Estimate ‘x’ At Intersection
    yv(k) = interp1(x(idxrng),y1(idxrng), xv(k));               % Interpolate To Estimate ‘y’ At Intersection
end

figure
plot(x, y1,    x, y2)
hold on
plot(xv, yv, 'pg', 'MarkerSize',15)
hold off
grid
legend('y_1', 'y_2', 'Intersections', 'Location','S')

I have many happy memories of flying C172s! Have fun!

[Rhiannon Elliott]

Hi! I'm not sure I really follow what is happening here. I'm relatively new to Matlab, unfortunately I've had little guidance this year due to Covid! What exactly is happening in this process to find the intersections?

I'm not sure of this helps give more context to my code but I have included the functions that were missing, and what the figure looks like:

function [P,rho,T,a] = atmosProperties(H, tempDiff, lapseDiff)

g = 9.8065;
R = 287.05287;
gamma = 1.4;
lNaught = -0.0065 + lapseDiff;
tNaught = 288.15 + tempDiff;
tEleven = 216.65 + tempDiff;
pNaught = 101325;
T= tNaught+(lNaught*H);
pEleven = pNaught*(1 + ((lNaught*11000)/tNaught))^(-g/(R*lNaught));

pressureTrop= pNaught*(1 + ((lNaught*H)/tNaught))^(-g/(R*lNaught));
pressureLStrat= pEleven*exp(1)^((-g/(R*tEleven))*(H-11000));
densTrop = pressureTrop/(R*T);
densLStrat = pressureLStrat/(R*tEleven);

if (H<11000) %troposphere
    rho = densTrop;
    P = pressureTrop;
    T=T;
    
elseif H >= 11000 %stratosphere
    rho = densLStrat;
    P = pressureLStrat;
    T = tEleven;
    
end

a = sqrt(gamma*R*T);

end


function [max_power_hp,max_power_W] = getPAmax (H)
hp_sealevel = 180;                                  %Maximum power at sealevel [hp] taken from Flight Gear source code

[~,rho_0,~,~] = atmosProperties (0,0,0);      %Sea level desnity [kgm^-3]
[~,rho_H,~,~] = atmosProperties (H,0,0);      %Density at given altitude [kgm^-3]
prop_eff = 0.85;                                    %Propeller efficiency

max_power_hp = prop_eff * hp_sealevel * (1.132*(rho_H/rho_0) - 0.132);
max_power_W = max_power_hp * 746;
end

[Star Strider]

My code requires two vectors that are functions of the same independent variable (or can be calculated to be). After that, simply subtract them and use my code to find the intersections.

I adapted my posted code to your problem, and added or changed this part of the code:

bclr = [c; ones(size(c))].' \ b(:);     % Linear Regression Coefficients
bvec = [V; ones(size(V))].' * bclr;     % Linear Regression Fit

curvdif = bvec - P_req(:);                                              % Subtract Curves
curvdif(isnan(curvdif)) = 1E4;
intidx = find(diff(sign(curvdif)));                                     % Approximate Indices Of Zero-Crossings

for k = 1:numel(intidx)
    idxrng = [max(1,intidx(k)-1):min(numel(V),intidx(k)+1)]             % Index Range For Interpolation
    xv(k) = interp1(curvdif(idxrng), V(idxrng), 0,'pchip','extrap');    % Interpolate To Estimate ‘V’ At Intersection
    yv(k) = interp1(V(idxrng),P_req(idxrng), xv(k));                    % Interpolate To Estimate ‘P_req’ At Intersection
end

fprintf(1,'\tIntersections:\n')
fprintf(1,'\t\tV = %9.4f, P_req = %9.1f\n', [xv(2:end); yv(2:end)])

%Goal of the figure is to find the intersections of the power available curve (flat line) and the power required curve (curved line)
%Legends will be added later on
plot (V, P_req);
hold on
plot (c, b);
plot(xv, yv, 'rx', 'MarkerSize',15)
xlabel("TAS [ms^{-1}]");
ylabel("Power reequired (W)")
legend('P_{req}', 'b', 'Intersections', 'Location','NW')
ylim([0 5E5])                                                   % Zoom To Show Intersections (Delete)

so the complete code is now:

%Code to determine the maximum and minimum speeds at a given altitude for
%the Cessna 172P

H=0;        %Eventually this will become an input in a function where the 
            % user just call the function and recieve the miniumum and maximum velocities for a given altitude 


%Constants of the aircraft
S = 16.16499443;                %Wing area [m^2]
g = 9.81;                       %acceleration due to gravity [ms^-2]
W = 887.0748299*g;              %Weight of the aircraft [kgf]
C_D0 = 0.026;                   %Zero lift drag coefficient 
K = 0.121;                      %Drag constant

%Range of velocities we wish to investigate [Mach]
M = linspace(0,0.8,100);


[~,rho,~,a]=atmosProperties(H,0,0);     %Density and speed of sound at given altitude [kgm^-3] [ms^-1]
[~,maxPA]=getPAmax(H);                  %Maximum power available at given altitude [W]
V = M.*a;                               %Velocity in [ms^-1]
C_L = (2./rho).*(W./S).*(1./(V.^2));    %Lift coefficient at given altitude
C_D = C_D0 + K .* C_L.^2;               %Drag coefficient at given altitude
P_req = 0.5 .* rho .* V.^3 .* S .* C_D;    %Power required curve

b = [maxPA, maxPA];                     %B and c are used to plot a line that corresponds to the maximum power available
c = [0,272.235190420871];

bclr = [c; ones(size(c))].' \ b(:);     % Linear Regression Coefficients
bvec = [V; ones(size(V))].' * bclr;     % Linear Regression Fit

curvdif = bvec - P_req(:);                                              % Subtract Curves
curvdif(isnan(curvdif)) = 1E4;
intidx = find(diff(sign(curvdif)));                                     % Approximate Indices Of Zero-Crossings

for k = 1:numel(intidx)
    idxrng = [max(1,intidx(k)-1):min(numel(V),intidx(k)+1)]             % Index Range For Interpolation
    xv(k) = interp1(curvdif(idxrng), V(idxrng), 0);                     % Interpolate To Estimate ‘V’ At Intersection
    yv(k) = interp1(V(idxrng),P_req(idxrng), xv(k));                    % Interpolate To Estimate ‘P_req’ At Intersection
end

fprintf(1,'\tIntersections:\n')
fprintf(1,'\t\tV = %9.4f, P_req = %9.1f\n', [xv(2:end); yv(2:end)])

%Goal of the figure is to find the intersections of the power available curve (flat line) and the power required curve (curved line)
%Legends will be added later on
plot (V, P_req);
hold on
plot (c, b);
plot(xv, yv, 'rx', 'MarkerSize',15)
xlabel("TAS [ms^{-1}]");
ylabel("Power reequired (W)")
legend('P_{req}', 'b', 'Intersections', 'Location','NW')
ylim([0 2E5])                                                           % Zoom To Show Intersections (Delete)

producing the plot with the intersections, and:

	Intersections:
		V =    8.1607, P_req =  114138.0
		V =   73.3073, P_req =  114138.0

Make appropriate changes to get the result you want.

.

[Rhiannon Elliott]

It makes more sense seeing it in context of my code, I can't thank you enough for your help!

[Star Strider]

As always, my pleasure!

Have fun flying your C172!