Fzero giving weird answer

1 次查看（过去 30 天）
Josh
Josh 2023-2-2
评论： Torsten 2023-2-3
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
  5 个评论
显示 3更早的评论
Josh
Josh 2023-2-2
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
Josh 2023-2-2
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

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Matt J
Matt J 2023-2-3
编辑：Matt J 2023-2-3
Ran in:
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.
  4 个评论
显示 2更早的评论
Josh
Josh 2023-2-3
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
Torsten 2023-2-3
@Matt J
At T = 463.54, CA_final and CB_final will have changed because of the ODE.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

产品

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by