How to find the intersection points on two functions

how to find the intersection points of dH_rem and dH_gen within the limits specified below

To=500:1:850; %outlet temp
Ti=570; %inlet temp
y_A=0.003;   %proportion of benzene in feed
V=8.5;
P=3e5; %Pa
R=8.3145; %kJ/mol.K
Cp=1.09 %kJ/kgK
Mr_air=29e-3; %kg/mol
dH_1=-1850; %kJ/mol
dH_2=-1423;
dH_3=-3273;
F_ao=0.1; %molar flow rate of benzene
k1 = 1e7.*exp(-12700./To);
k2 = 5e4.*exp(-10800./To);
k3 = 7e7.*exp(-15000./To);
t=(V.*y_A.*P)./(F_ao.*R.*To);
X=(t.*(k1+k3))./(1+t.*(k1+k3));
S=k1./((1+(k2.*t)).*(k1+k3));
Y_B = (k1.*t)./(((k2.*t)+1).*(1+(t.*(k1+k3))));

dH_gen=-((V.*y_A.*P)./(R.*To)).*((((k1.*dH_1)+(k3.*dH_3)).*(1-X))+(k2.*X.*S.*dH_2));
plot(To,dH_gen)
hold on
dH_rem=(Mr_air./y_A).*F_ao.*Cp.*(To-Ti);
plot(To,dH_rem)
hold off

Answer 1

[John D'Errico]

Standard question: how to find the intersection(s) of two functions.

1. Subtract them. Where the difference is zero, there lies an intersection.

2. Use a root finder. That could be anything from fzero, solve, vpasolve, fsolve, etc.

Note that all standard optimization based root finders will find ONE root, and only one root. It will depend on your starting values. Solve might be able to find the three points of intersection I saw on the plot.

The lazy solution to finding an approximate root is to use a tool like Doug Schwarz's intersections code on the FEX. Evaluate each function at a few hundred points, then call intersections. It will use linear interpolation to find the crossing points. But a virtue of that solution is it will report all crossings between the curves.

[sophp]

I shortened the limited between 500-650 as I know from the graph that only one solution lies between these limits. Then I put after the code given in the question

syms To
u=solve(dH_gen-dH_rem == 0, To)

The output was

Empty sym: 0-by-1

Implying there is no solution, why is this the case?

[Walter Roberson]

When you set To to numeric at the beginning of the code, you created a numeric vector for dH_gen and dH_rem: the numeric values of To are copied as needed to form the expressions. Then you assign a symbol to To and suddenly expect that you can solve() the numeric vector for a symbol.

You need to do the syms To in place of the To=500:1:850;

If you also syms Ti in place of assigning 370 to it, and then go through the solve() you will get

dH_rem =
(3161*To)/3000 - (3161*Ti)/3000

You can see at a glance that the equation will be satisfied when To == Ti .

Notice that you compute Y_B but never use it.

[sophp]

Thank you!! Okay, so I replaced

To=500:1:850

with

syms To

Matlab then tells me that the equation dH_rem-dH_gen is unsolvable symbolically and gives me a numerical approximation instead.

The solution is correct but how do I alter the code to give the lowest intersection value.

[Walter Roberson]

When I use the code you posted, except replacing that assignment with syms To, then dH_rem = (3161*To)/3000 - 60059/100 which is linear with the single solution To == 570. Perhaps you are using different code by now?