Trying to find an x value from a certain y value on my graph

I have a graph of reactor temperature against time and I am trying to find the time the reactor temperature is at 455K.

To make this graph I used ode15 function [t,y] = ode15s(@(t,y)fun(t,y(1),y(2),y(3),y(4),y(5),y(6)), tspan, initial_conditions);

then plotted the graph using plot: figure(3) plot(t,T_2)

Does anyone know how to write code that would be able to find the time value?

I would appreciate any help, thanks.

[KSSV]

How the plot is? How plot(t,T_2) looks?

[Tom Goodland]

Answer 1

[KSSV]

idx = knnsearch(T_2,455) ;
t(idx)

If there is one-to-one mapping between (t,T_2)

iwant = interp1(T_2,t,455)

[Tom Goodland]

Thanks so much it worked.

[Tom Goodland]

Now I have to use this time value I just calculated using your code (t=0.9378) and use the t value to find ca,cb,cs and cd value on a graph i have so I know all initial conditions. Do you know how I could do this?Here's the graph for ca,cb,cs and cd against time:

[KSSV]

Try interp1 again.

[Tom Goodland]

I used it again for the concentration plot, it worked the first time I did it with Ca, however when I tried to use it cb, cs and cd I got some errors, could this be due to the one-on-one mapping you mentioned?

plot(t,[ca,cb,cs,cd])

idx = knnsearch(0.9378,ca) ;

t(idx)

iwant = interp1(ca,t,0.9378);% ca = 1.1587

idx = knnsearch(0.9378,cb) ;

t(idx)

iwant = interp1(cb,t,0.9378);

idx = knnsearch(0.9378,cs) ;

t(idx)

iwant = interp1(cs,t,0.9378);

idx = knnsearch(0.9378,cd) ;

t(idx)

iwant = interp1(cd,t,0.9378);

The errors I got were:

Error using matlab.internal.math.interp1

Sample points must be unique.

Error in interp1 (line 188)

VqLite = matlab.internal.math.interp1(X,V,method,method,Xqcol);

Error in Twob (line 21)

iwant = interp1(cb,t,0.9378);

Answer 2

[Star Strider]

I will re-post my original Answer because I believe it addresses the problem (although without the data, it is difficult to determine that with certainty) —

t = linspace(0, 10);
y = sin(2*pi*t/3) + t/10;

yval = 1;
idx = find(diff(sign(y-yval)));                                         % Approximate Index Values For 'yval'
for k = 1:numel(idx)
        idxrng = [max([1 idx(k)-2]) : min([numel(t) idx(k)+2])];        % Index Range For Interpolation
        tval(k) = interp1(y(idxrng), t(idxrng), yval);                  % Interpolate
end
tval

tval = 1×7

0.5893 0.9546 3.3487 4.2029 6.1877 7.3731 9.0459

figure
plot(t,y)
hold on
plot(tval, ones(size(tval))*yval, '+r', 'MarkerSize',7.5)
hold off
grid

.

[Tom Goodland]

Thanks for the repost, the most recent graph I posted is of multiple different component concentrations agaisnt time. Every line doesn't have any repeating y values however, I just need to know the y value (conc.) for a certain x value (time) for each line. Do you know how to do this I would really appreciate any help.

[Star Strider]

Please provide the data vectors (or the code that creates them) and I will see what I can do.

If data, a text or Excel file is preferable.

If code, I would prefer the code as text pasted into a Comment (or original post) rather than an attached .m file that I then have to download and copy over to here to run. I only need the relevant parts of the code necessary to do my best to solve the problem — relevant variables, functions, the ode15s call and subsequent plot calls.

.

[Tom Goodland]

I didn't want to miss anything out so this is the whole script that makes all the graphs:

function main

tspan = [0 200]; %just as an example

initial_conditions = [4.3 5.1 3 0 4.4 422]; %its the values of the variables you are solving for at t=0

%there must be as many inputs to the initial conditions as variables solved

%Solve your ODE function (I picked ode45 in this case,could be another):

[t,y] = ode15s(@(t,y)fun(t,y(1),y(2),y(3),y(4),y(5),y(6)), tspan, initial_conditions);

%output t is your time range, and y an array with all your variables where:

ca = y(:,1);

cb = y(:,2);

cs = y(:,3);

cd = y(:,4);

P = y(:,5);

T_2 = y(:,6);

figure(1)

plot(t,[ca,cb,cs,cd])

figure(2)

plot(t,P)

figure(3)

plot(t,T_2)

end

function out = fun(tspan,ca,cb,cs,cd,P,T_2)

%Get your inputs and assign them to variables:

%(keep the same order as in your initial conditions!)

%Define your relationships:

k_0_1 = 4E14;

k_0_2 = 1E84;

Ea1 = 128000;

Ea2 = 800000;

R_gas = 8.314;

V0 = 4000;% L

VH = 5000;

Cv1 = 3360; % mol / h * atm

Cv2 = 53600; % mol / h * atm

Hr1 = -45400; % J / mol

Hr2 = -3.2E5;

Ta = 373;

UA = 0;

NCp = 1.26E7; % J / K

%choose your values of k1 & k2:

k1 = k_0_1*exp(-Ea1/(R_gas*T_2));

k2 = k_0_2*exp(-Ea2/(R_gas*T_2));

%Reaction rates:

r1 = k1*ca*cb;

r2 = k2*cs;

%Fd

Fd = (-0.5*r1-3*r2)*V0;

%

if Fd < 11.4 % mol / h

Fvent = Fd;

else

if P < 28.2 % atm

Fvent = (P - 1)*Cv1;

else

Fvent = (P - 1)*(Cv1 + Cv2);

end

end

%Define your ODEs:

dca = -r1;

dcb = -r1;

dcs = -r2;

dcd = r1/2 + 3*r2;

dp = (Fd - Fvent)*((R_gas*T_2)/VH);

dT = (V0*(r1*-Hr1 + r2*-Hr2) - UA*(T_2 - Ta))/NCp;

%Assign your ourput:

out = [dca; dcb; dcs; dcd;dp;dT];

end

Thanks for the help

[Torsten]

Define an event that interrupts integration when T reaches 455 K.

Then you have all variables available at that time.

[Star Strider]

What do you want to get, and from what plot or plots? (I added the legend, axis labels, and figure number titles to make them easier to refer to.)

main

Warning: Failure at t=3.426260e+00. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (7.105427e-15) at time t.

T_2 = 455K at 0.92355

At t = 0.92355:
	ca	 = 2.01035
	cb	 = 2.81035
	cs	 = 3.00000
	cd	 = 1.14482

function main
tspan = [0 200]; %just as an example
initial_conditions = [4.3 5.1 3 0 4.4 422]; %its the values of the variables you are solving for at t=0
%there must be as many inputs to the initial conditions as variables solved
%Solve your ODE function (I picked ode45 in this case,could be another):
[t,y] = ode15s(@(t,y)fun(t,y(1),y(2),y(3),y(4),y(5),y(6)), tspan, initial_conditions);
%output t is your time range, and y an array with all your variables where:
ca = y(:,1);
cb = y(:,2);
cs = y(:,3);
cd = y(:,4);
P = y(:,5);
T_2 = y(:,6);
figure(1)
plot(t,[ca,cb,cs,cd])
xlabel('t')
ylabel('Concentrations (?)')
title('Figure 1')
legend('ca','cb','cs','cd', 'Location','best')
figure(2)
plot(t,P)
xlabel('t')
ylabel('P')
title('Figure 2')
figure(3)
plot(t,T_2)
xlabel('t')
ylabel('T_2')
title('Figure 3')

T455K_t = interp1(T_2, t, 455);
ca_cb_cs_cd = interp1(t, [ca,cb,cs,cd], T455K_t);
fprintf('\n\nT_2 = 455K at %.5f\n\nAt t = %.5f:\n\tca\t = %.5f\n\tcb\t = %.5f\n\tcs\t = %.5f\n\tcd\t = %.5f\n',T455K_t,T455K_t,ca_cb_cs_cd)

end

function out = fun(tspan,ca,cb,cs,cd,P,T_2)
%Get your inputs and assign them to variables:
%(keep the same order as in your initial conditions!)
%Define your relationships:      
k_0_1 = 4E14;
k_0_2 = 1E84;
Ea1 = 128000;
Ea2 = 800000;
R_gas = 8.314;
V0 = 4000;% L
VH = 5000;
Cv1 = 3360; % mol / h * atm
Cv2 = 53600; % mol / h * atm
Hr1 = -45400; % J / mol
Hr2 = -3.2E5;
Ta = 373;
UA = 0;
NCp = 1.26E7; % J / K
%choose your values of k1 & k2:
k1 = k_0_1*exp(-Ea1/(R_gas*T_2));
k2 = k_0_2*exp(-Ea2/(R_gas*T_2));
%Reaction rates:
r1 = k1*ca*cb;
r2 = k2*cs;
%Fd
Fd = (-0.5*r1-3*r2)*V0;
%
if Fd < 11.4 % mol / h
    Fvent = Fd;
else 
    if P < 28.2 % atm
        Fvent = (P - 1)*Cv1;
    else
        Fvent = (P - 1)*(Cv1 + Cv2);
    end
end
%Define your ODEs:
dca = -r1;
dcb = -r1;
dcs = -r2;
dcd = r1/2 + 3*r2;
dp = (Fd - Fvent)*((R_gas*T_2)/VH);
dT = (V0*(r1*-Hr1 + r2*-Hr2) - UA*(T_2 - Ta))/NCp;
%Assign your ourput:
out = [dca; dcb; dcs; dcd;dp;dT];
end

EDIT — (17 Jan 2022 at 14:54)

Added these lines at the end of ‘main’ —

T455K_t = interp1(T_2, t, 455);

ca_cb_cs_cd = interp1(t, [ca,cb,cs,cd], T455K_t);

fprintf('\n\nT_2 = 455K at %.5f\n\nAt t = %.5f:\n\tca\t = %.5f\n\tcb\t = %.5f\n\tcs\t = %.5f\n\tcd\t = %.5f\n',T455K_t,T455K_t,ca_cb_cs_cd)

.

[Tom Goodland]

Thanks for that star strider, I want to find out what ca, cb, cs and cd are when T_2 = 455 K, I've already used someone's advice previously to find out the time T_2 = 455 K which is t = 0.9378 hours.

[Tom Goodland]

Thorsten I was looking at the ode event page on matlab last night and got quite confused, could you help define this event when T = 455 K please?

[Torsten]

https://de.mathworks.com/matlabcentral/answers/1628960-help-with-odes-and-two-different-reactions?s_tid=srchtitle