Error in my ode45 equation.

조회 수: 3 (최근 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일
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
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
Zara Freeman
Zara Freeman 2020년 12월 29일
Thank you ! I managed to get it to work
Star Strider
Star Strider 2020년 12월 29일
My pleasure!

댓글을 달려면 로그인하십시오.

카테고리

Help CenterFile 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