Fzero giving weird answer

2023 2 2
1 个回答
更新时间：2023 2 3
1 次查看（30 天）

0 个投票

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更早的评论 隐藏 3更早的评论

Walter Roberson
Walter Roberson 2023-2-2
Please enough for us to be able to execute to test.
Matt J
Matt J 2023-2-2
编辑：Matt J 2023-2-2
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
Torsten 2023-2-2
编辑：Torsten 2023-2-2
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
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:

0 个投票

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更早的评论 隐藏 2更早的评论

Josh
Josh 2023-2-3
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
Matt J 2023-2-3
编辑：Matt J 2023-2-3
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
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.

请先登录，再进行评论。

类别

在 帮助中心 和 File Exchange 中查找有关 Programming 的更多信息

产品

标签

提问:

Josh
Josh
2023-2-2

评论：

Torsten
Torsten
2023-2-3

Community Treasure Hunt

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

Start Hunting!

Translated by