Error in my ode45 equation.

7 次查看（过去 30 天）
Zara Freeman
Zara Freeman 2020-12-29
评论： Star Strider 2020-12-29
I am unsure why my code is wrong (pasted below) as using a very similar example this code worked. When running the code below, these errors are displayed:
Error in ode45 (line 115)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
Error in reactorassignmentpart1plot (line 5)
[V,y] = ode45(@fun1, Vspan, yo);
This is the code used:
clc
Vspan = [0 1];
yo = [4 0 0 1 500.15];
[V,y] = ode45(@fun1, Vspan, yo);
subplot(2,1,1)
plot(V,y(:,1),V,y(:,2),V,y(:,3));
legend('Fa','Fb','Fc');
ylabel('Molar flowrate, mol/min');
xlabel('Volume,dm3');
subplot(2,1,2)
plot(V,y(:,4));
legend('Temperature (K)');
ylabel('Temperature (K)');
xlabel('Volume,dm3');
% Compose the function
function f = fun1(V,Y)
% Define the differential equations that need to be solved
% Y is the concentration and V is the PFR volume
Fa = Y(1);
Fb = Y(2);
Fc = Y(3);
T = Y(4);
% Define initial conditions
deltaH1 = -25000; %kJ/molA
deltaH2 = 35000; %kJ/molB
deltaH2T = 35000-(80*(T-298)); %kJ/molB
CTo = 0.3996; % mol/L
To = 500.15; % K
Fio = 1; % mol/min
FTo = 4; % mol/min
Cpa = 20; % J/molK
Cpb = 80; % J/molK
Cpc = 100; % J/molK
Cpi = 20; % J/molK
Ua = 150; %J/dm3.min.K
FT = Fa + Fb + Fc + Fio;
Ca = CTo*((Fa/FT)*(To/T));
Cb = CTo*((Fb/FT)*(To/T));
Cc = CTo*((Fc/FT)*(To/T));
Ci = CTo*((Fio/FT)*(To/T));
K1(T) = 50*exp((8000/8.314)*((1/315)-(1/T))); % dm3/mol.min
Kc(T) = 10*exp((-25/8.314)*((1/315)-(1/T))); % dm3/mol.min
K2(T) = 400*exp((4000/8.314)*((1/310)-(1/T))); % dm3/mol.min-1
Ta = 523.15; % K
ra = (K1*((Ca^2)-((1/Kc)*Cb))) + (K2*Ca*(Cb^2));
rb = 1/2*((K1*((Ca^2)-((1/Kc)*Cb)))+(2*K2*(Cb^2)*Ca));
rc = K2*Ca*(Cb^2);
% Differential equations
dFadV = -ra;
dFbdV = -rb;
dFcdV = rc;
dTdV = (Ua*(Ta-T)-((ra-rc)*(deltaH1))-(2*rc*(deltaH2)))/(Cpa*Fa+Cpb*Fb+Cpc*Fc+Cpi*Fio);
f = [dFadV; dFbdV; dFcdV; dTdV];
end
Thank you !

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

回答（3 个）

Walter Roberson
Walter Roberson 2020-12-29
Ran in:
clc
Vspan = [0 1];
yo = [4 0 0 1 500.15];
[V,y] = ode45(@fun1, Vspan, yo);
Error using odearguments (line 95)
FUN1 returns a vector of length 4, but the length of initial conditions vector is 5. The vector returned by FUN1 and the initial conditions vector must have the same number of elements.

Error in ode45 (line 115)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
Sure enough, you are passing in a yo of length 5,
subplot(2,1,1)
plot(V,y(:,1),V,y(:,2),V,y(:,3));
legend('Fa','Fb','Fc');
ylabel('Molar flowrate, mol/min');
xlabel('Volume,dm3');
subplot(2,1,2)
plot(V,y(:,4));
legend('Temperature (K)');
ylabel('Temperature (K)');
xlabel('Volume,dm3');
% Compose the function
function f = fun1(V,Y)
% Define the differential equations that need to be solved
% Y is the concentration and V is the PFR volume
Fa = Y(1);
Fb = Y(2);
Fc = Y(3);
T = Y(4);
% Define initial conditions
deltaH1 = -25000; %kJ/molA
deltaH2 = 35000; %kJ/molB
deltaH2T = 35000-(80*(T-298)); %kJ/molB
CTo = 0.3996; % mol/L
To = 500.15; % K
Fio = 1; % mol/min
FTo = 4; % mol/min
Cpa = 20; % J/molK
Cpb = 80; % J/molK
Cpc = 100; % J/molK
Cpi = 20; % J/molK
Ua = 150; %J/dm3.min.K
FT = Fa + Fb + Fc + Fio;
Ca = CTo*((Fa/FT)*(To/T));
Cb = CTo*((Fb/FT)*(To/T));
Cc = CTo*((Fc/FT)*(To/T));
Ci = CTo*((Fio/FT)*(To/T));
K1(T) = 50*exp((8000/8.314)*((1/315)-(1/T))); % dm3/mol.min
Kc(T) = 10*exp((-25/8.314)*((1/315)-(1/T))); % dm3/mol.min
K2(T) = 400*exp((4000/8.314)*((1/310)-(1/T))); % dm3/mol.min-1
Ta = 523.15; % K
ra = (K1*((Ca^2)-((1/Kc)*Cb))) + (K2*Ca*(Cb^2));
rb = 1/2*((K1*((Ca^2)-((1/Kc)*Cb)))+(2*K2*(Cb^2)*Ca));
rc = K2*Ca*(Cb^2);
% Differential equations
dFadV = -ra;
dFbdV = -rb;
dFcdV = rc;
dTdV = (Ua*(Ta-T)-((ra-rc)*(deltaH1))-(2*rc*(deltaH2)))/(Cpa*Fa+Cpb*Fb+Cpc*Fc+Cpi*Fio);
f = [dFadV; dFbdV; dFcdV; dTdV];
but your f is only putting together 4 values.
end
Every yo input value is a boundary condition. You have to return the derivative for each of the entries.
You do not seem to use Y(5) in your fun1, so perhaps you should only be passing in 4 initial values instead of 5.
Mischa Kim
Mischa Kim 2020-12-29
编辑：Mischa Kim 2020-12-29
Hi Zara,
your vector of initial conditions has 5 components
yo = [4 0 0 1 500.15];
however, you only have 4 differential equations:
f = [dFadV; dFbdV; dFcdV; dTdV];
The numbers need to match, so, you are either missing a differential equation or your yo vector is too long.
Also, in your ode file you probably want to change K1(t), Kc(t), and K2(t) to K1, Kc, and K2.
  3 个评论
显示 1更早的评论
Mischa Kim
Mischa Kim 2020-12-29
Difficult to say what the issue is without digging further into the differential equation and your code. You can try to play with different solvers. And I also suggest double-checking your code. Sometimes it is one little plus vs minus sign (or vice versa) that makes all the difference. Feel free to share a screen shot of the equations.
Zara Freeman
Zara Freeman 2020-12-29
I got it to work! Thank you !

请先登录，再进行评论。

Star Strider
Star Strider 2020-12-29
The most significant problem is that you need to define the functions as anonymous functions:
K1 = @(T) 50*exp((8000/8.314)*((1/315)-(1/T))); % dm3/mol.min
Kc = @(T) 10*exp((-25/8.314)*((1/315)-(1/T))); % dm3/mol.min
K2 = @(T) 400*exp((4000/8.314)*((1/310)-(1/T))); % dm3/mol.min-1
and then call them with the appropriate arguments:
ra = (K1(T)*((Ca^2)-((1/Kc(T))*Cb))) + (K2(T)*Ca*(Cb^2));
rb = 1/2*((K1(T)*((Ca^2)-((1/Kc(T))*Cb)))+(2*K2(T)*(Cb^2)*Ca));
rc = K2(T)*Ca*(Cb^2);
See the documentation section on Anonymous Functions for details.
The other problem is that you have 4 differential equations:
f = [dFadV; dFbdV; dFcdV; dTdV];
and 5 initial conditions:
yo = [4 0 0 1 500.15];
Correct these and your code runs, however it gives a Warning:
Warning: Failure at t=7.629620e-01. Unable to meet integration tolerances without
reducing the step size below the smallest value allowed (1.776357e-15) at time t.
> In ode45 (line 360)
meaning that it has found a singularity and cannot go further. (I may not have edited ‘yo’ correctly to make the intiial conditions equal to the number of differential equations, so check that.)
The corrected ‘fun1’ code is:
function f = fun1(V,Y)
% Define the differential equations that need to be solved
% Y is the concentration and V is the PFR volume
Fa = Y(1);
Fb = Y(2);
Fc = Y(3);
T = Y(4);
% Define initial conditions
deltaH1 = -25000; %kJ/molA
deltaH2 = 35000; %kJ/molB
deltaH2T = 35000-(80*(T-298)); %kJ/molB
CTo = 0.3996; % mol/L
To = 500.15; % K
Fio = 1; % mol/min
FTo = 4; % mol/min
Cpa = 20; % J/molK
Cpb = 80; % J/molK
Cpc = 100; % J/molK
Cpi = 20; % J/molK
Ua = 150; %J/dm3.min.K
FT = Fa + Fb + Fc + Fio;
Ca = CTo*((Fa/FT)*(To/T));
Cb = CTo*((Fb/FT)*(To/T));
Cc = CTo*((Fc/FT)*(To/T));
Ci = CTo*((Fio/FT)*(To/T));
K1 = @(T) 50*exp((8000/8.314)*((1/315)-(1/T))); % dm3/mol.min
Kc = @(T) 10*exp((-25/8.314)*((1/315)-(1/T))); % dm3/mol.min
K2 = @(T) 400*exp((4000/8.314)*((1/310)-(1/T))); % dm3/mol.min-1
Ta = 523.15; % K
ra = (K1(T)*((Ca^2)-((1/Kc(T))*Cb))) + (K2(T)*Ca*(Cb^2));
rb = 1/2*((K1(T)*((Ca^2)-((1/Kc(T))*Cb)))+(2*K2(T)*(Cb^2)*Ca));
rc = K2(T)*Ca*(Cb^2);
% Differential equations
dFadV = -ra;
dFbdV = -rb;
dFcdV = rc;
dTdV = (Ua*(Ta-T)-((ra-rc)*(deltaH1))-(2*rc*(deltaH2)))/(Cpa*Fa+Cpb*Fb+Cpc*Fc+Cpi*Fio);
f = [dFadV; dFbdV; dFcdV; dTdV];
end
.
  4 个评论
显示 2更早的评论
Zara Freeman
Zara Freeman 2020-12-29
Thank you ! I managed to get it to work

请先登录，再进行评论。

类别

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

标签

Community Treasure Hunt

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

Start Hunting!

Translated by