Fzero giving weird answer

I'm trying to find a maximum upper bound for an ODE by using ode15s, and fzero to find root of function that will determine the maximum of a variable expressed by the ODE. I have used a while loop to reiterate using absolute error values until it works however it runs forever. On closer inspection its because the first iteration (even before while loop) returns the wrong "solution" through fzero. Am i doing something fundamentally wrong or should I use another way to solve this?

%% at point where Q_generated = Q_removed , func = 0 @ Q_generated - Q_removed
%% initalising guess and coressponding outcomes at that time, 

guess = 46 ; %mins
CA0 = (ni_incident(1))/V_incident ;
CB0 = (ni_incident(2))/V_incident ;
T0 = 175+273.15 ;

Y0(1) = CA0 ;
Y0(2) = CB0 ;
Y0(3) = T0 ;

[X,Y] = ode15s(@(t,y)DYDT_incident_C(t,y),[0 guess],Y0) ;

CA_final = Y(end,1) ;
CB_final = Y(end,2) ;
T_final = Y(end,3) ;

func = @(Temp)((V_incident * ((-1*(k0*(exp(((-1*Ea)/(R*Temp))))))/(CA_final*CB_final)) * H_reaction)-(UA * (Temp - T_cooling_water))) ;
solution = fzero(func,T_final) ;
error = solution - T_final ;
%% while loop is created to progressively re-iterate increasing upper bound values of time for ode solver until solution satisfied

while abs(error) > (0.1)
    guess = guess + (1/60) ;
    [X,Y] = ode15s(@(t,y)DYDT_incident_C(t,y),[45 guess],Y0) ;
    CA_final = Y(end,1) ;
    CB_final = Y(end,2) ;
    T_final = Y(end,3) ;

    func = @(Temp)((V_incident * (-1*(k0*(exp(((-1*Ea)/(R*Temp))))))/(CA_final*CB_final) * H_reaction)-(UA * (Temp - T_cooling_water))) ;
    solution = fzero(func,T_final);
    error = solution - T_final ;
end

function dydt = DYDT_incident_C(t,y)

    global V_incident H_reaction sum_niCp k0 Ea R

    dydt = zeros(3,1) ;

    CA = y(1) ;
    CB = y(2) ;
    T = y(3) ;

    k = k0*(exp(-((Ea)/(R*T)))) ;

    %% dCAdt = rA
    rA = (-k)*CA*CB ;

    %% dCBdt = 2*rA = rB
    rB = 2*rA ;

    if t > 45
        Q_gen = (V_incident * rA * H_reaction) ;
        Q_rem = 0 ;
    else
        Q_gen = 0 ;
        Q_rem = 0 ;
    end
   
    %% Temp variation over time
    dT_dt = ((Q_gen-Q_rem)./(sum_niCp)) ;
    
    dydt(1) = rA ;
    dydt(2) = rB ;
    dydt(3) = dT_dt ;

end

[Walter Roberson]

Please enough for us to be able to execute to test.

[Matt J]

We can't run your code because it relies on a number of variables that haven't been provided.

Aside from that, you need to describe in what way the result of fzero is "wrong". Did you check whether func(solution)=0? Did you plot func() to see if it even has any roots?

[Torsten]

You should use the "Event" facility of the ODE integrators and stop integration at the time when y(3) is a root of your algebraic equation you now solve with fzero. This way, no iteration from your side is necessary.

[Josh]

I will take a look into bvp4c, thanks.

I have solved the problem on my graphical calculator - this is how i determined the "solution to be wrong".

I will put the code (in its runable entirety in the chat).

[Josh]

global T_cooling_water ni_normal V_normal ni_incident V_incident

%% cooling water temperature in heat exchanger
T_cooling_water = 25 + 273.15 ; % K

%% initial molar qunatities of reactants under normal operation
ni_ONCB_normal = 3.17E3 ; % mol
ni_NH3_normal = 33E3 ; % mol
ni_H2O_normal = 103.6E3 ; % mol
ni_normal = [ni_ONCB_normal ni_NH3_normal ni_H2O_normal] ;

%% reaction volume under these conditions, in litres
V_normal = 3.265E3 ; % L

%% initial molar qunatities of reactants under circumstances of incident
ni_ONCB_incident = 9.044E3 ; % mol
ni_NH3_incident = 33E3 ; % mol
ni_H2O_incident = 103.6E3 ; % mol
ni_incident = [ni_ONCB_incident ni_NH3_incident ni_H2O_incident] ;

%% reaction volume under these conditions
V_incident = 5.119E3 ; % L

global Cp UA Ea H_reaction R k0 sum_niCp

%% heat capcities of reactants 
Cp_ONCB = 167.4 ; % J mol^-1 K^-1
Cp_NH3 = 35.1 ; % J mol^-1 K^-1
Cp_H2O = 75.3 ; % J mol^-1 K^-1
Cp = [Cp_ONCB Cp_NH3 Cp_H2O] ;

%% sum of initial molar quantities and heat capcities under incident conditions
sum_niCp = sum(ni_incident.*Cp) ;

%% value of product of overall heat transfer coeffcient and heat transfer surface area 
UA = 1.51E5 ; % J min^-1 K^-1

%% activation energy, heat of reaction, gas constant and rate constant values 
Ea = 47166 ; % J mol^-1
H_reaction = -2.469E6 ; % J mol^-1
R = 8.314 ; % J mol^-1 K^-1
k_ref = 1.7E-4 ;
k0 = (k_ref)/(exp(-((Ea)/(R*(188+273.15))))) ; % L mol^-1 min^-1

%% at point where Q_generated = Q_removed , func = 0 @ Q_generated - Q_removed
%% initalising guess and coressponding outcomes at that time, 

guess = 50 ; %mins
CA0 = (ni_incident(1))/V_incident ;
CB0 = (ni_incident(2))/V_incident ;
T0 = 175+273.15 ;

Y0(1) = CA0 ;
Y0(2) = CB0 ;
Y0(3) = T0 ;

[X,Y] = ode15s(@(t,y)DYDT_incident_C(t,y),[0 guess],Y0) ;

CA_final = Y(end,1) ;
CB_final = Y(end,2) ;
T_final = Y(end,3) ;

func = @(Temp)((V_incident * ((-1*(k0*(exp((-1*Ea)/(R*Temp)))))/(CA_final*CB_final)) * H_reaction)-(UA * (Temp - T_cooling_water))) ;
solution = fzero(func,T_final) ;
error = solution - T_final ;

%% while loop is created to progressively re-iterate increasing upper bound values of time for ode solver until solution satisfied

while abs(error) > (0.1)
    guess = guess + (1/60) ;
    [X,Y] = ode15s(@(t,y)DYDT_incident_C(t,y),[45 guess],Y0) ;
    CA_final = Y(end,1) ;
    CB_final = Y(end,2) ;
    T_final = Y(end,3) ;

    func = @(Temp)((V_incident * (-1*(k0*(exp(((-1*Ea)/(R*Temp))))))/(CA_final*CB_final) * H_reaction)-(UA * (Temp - T_cooling_water))) ;
    solution = fzero(func,T_final);
    error = solution - T_final ;
end

function dydt = DYDT_incident_C(t,y)

    global V_incident H_reaction sum_niCp k0 Ea R

    dydt = zeros(3,1) ;

    CA = y(1) ;
    CB = y(2) ;
    T = y(3) ;

    k = k0*(exp(-((Ea)/(R*T)))) ;

    %% dCAdt = rA
    rA = (-k)*CA*CB ;

    %% dCBdt = 2*rA = rB
    rB = 2*rA ;

    if t > 45
        Q_gen = (V_incident * rA * H_reaction) ;
        Q_rem = 0 ;
    else
        Q_gen = 0 ;
        Q_rem = 0 ;
    end
   
    %% Temp variation over time
    dT_dt = ((Q_gen-Q_rem)./(sum_niCp)) ;
    
    dydt(1) = rA ;
    dydt(2) = rB ;
    dydt(3) = dT_dt ;

end

Answer 1

[Matt J]

As you can see from the plot, the root lies quite close to solution, so fzero succeeded. If this violates your expectations, it is likely that you haven't implemented the intended equations.

guess = 46 ; %mins
CA0 = (ni_incident(1))/V_incident ;
CB0 = (ni_incident(2))/V_incident ;
T0 = 175+273.15 ;
Y0(1) = CA0 ;
Y0(2) = CB0 ;
Y0(3) = T0 ;
[X,Y] = ode15s(@(t,y)DYDT_incident_C(t,y),[0 guess],Y0) ;
CA_final = Y(end,1) ;
CB_final = Y(end,2) ;
T_final = Y(end,3) ;
func = @(Temp)((V_incident * ((-1*(k0*(exp(((-1*Ea)/(R*Temp))))))/(CA_final*CB_final)) * H_reaction)-(UA * (Temp - T_cooling_water))) ;
solution = fzero(func,T_final)

solution = 298.1516

fplot(@(x)func(x+solution),[-.1,.1])

Warning: Function behaves unexpectedly on array inputs. To improve performance, properly vectorize your function to return an output with the same size and shape as the input arguments.

[Josh]

The equation is simply a reaarangement of terms you can see in the function file, on my GDC the solution for fzero comes @ t= 53.58 which corresponds to T = 463.54

[Matt J]

Plots don't lie. Neither does substitution. Let's plug T=463.54 into your function and see what we get:

>> func(463.54)

ans =

-2.476175486149412e+07

Clearly T=463.54 is not a root of the function you've provided.

[Josh]

I appreciae your help, really not sure why i cant get it to gel (other than the obvious provided info ie not a root) but the equation is correct im sure of that much.

[Torsten]

@Matt J

At T = 463.54, CA_final and CB_final will have changed because of the ODE.