You are now following this question
- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.
How can I solve this differential equation?
4 views (last 30 days)
Show older comments
Hello, I'm trying to solve the differential equation however, I have some problem during solving my problem.
clear, clc
format long
syms c1(t) rho_cat V a1 a2 a3 a4 c10 R T t
Eq1 = diff(c1) == (-(rho_cat / V) * a1* exp(-a2 / (R * (T+273.15))) * (c1^a3)) / (1+a4*c1);
[C1] = dsolve(Eq1,c1(0) == c10)
C1 = matlabFunction([C1], 'Vars',{[a1,a2,a3,a4],t,c10,rho_cat,V, R, T})
When I try to get c1 using dsolve, and define c1 as a matlab function, the matlab gets error with this message.
Warning: Function '_minus' not verified to be a valid MATLAB function.
Warning: Finite sets ('DOM_SET') not supported. Using element '0' instead.
I think the reason of warning is because C1 has a form like below when I use dsolve to get C1
solve([c1^(1 - a3)*(a3 - a4*c1 + a3*a4*c1 - 2) == ((c10*(a3 - a4*c10 + a3*a4*c10 - 2))/(c10^a3*(a3^2 - 3*a3 + 2)) + (a1*rho_cat*t*exp(-(20*a2)/(5463*R + 20*R*T)))/V)*(a3^2 - 3*a3 + 2), c1 ~= 0], c1) minus {0}
Do I have to use other method to solve this problem like ode45 or ode23 instead of dsolve?
I hope I can get a good reply.
Thanks.
Accepted Answer
Walter Roberson
on 28 Dec 2021
Something like this:
syms c1(t) rho_cat V a1 a2 a3 a4 c10 R T t
Eq1 = diff(c1) == (-(rho_cat / V) * a1* exp(-a2 / (R * (T+273.15))) * (c1^a3)) / (1+a4*c1);
[C1] = dsolve(Eq1,c1(0) == c10)
C1 =
C11 = children(children(C1,1),1);
C1f = lhs(C11)-rhs(C11)
C1f =
V1 = [a1,a2,a3,a4]
V1 =
V2 = {t,c10,rho_cat,V, R, T}
V2 = 1×6 cell array
{[t]} {[c10]} {[rho_cat]} {[V]} {[R]} {[T]}
V3 = setdiff(setdiff(symvar(C1f), V1), [V2{:}])
V3 =
vars = {V3, V1, V2{:}}
vars = 1×8 cell array
{[c1]} {[a1 a2 a3 a4]} {[t]} {[c10]} {[rho_cat]} {[V]} {[R]} {[T]}
C1inner = matlabFunction(C1f, 'Vars', vars)
C1inner = function_handle with value:
@(c1,in2,t,c10,rho_cat,V,R,T)[-((c10.*c10.^(-in2(:,3)).*(in2(:,3)-in2(:,4).*c10+in2(:,3).*in2(:,4).*c10-2.0))./(in2(:,3).*-3.0+in2(:,3).^2+2.0)+(in2(:,1).*rho_cat.*t.*exp((in2(:,2).*-2.0e+1)./(R.*5.463e+3+R.*T.*2.0e+1)))./V).*(in2(:,3).*-3.0+in2(:,3).^2+2.0)+c1.^(-in2(:,3)+1.0).*(in2(:,3)-in2(:,4).*c1+in2(:,3).*in2(:,4).*c1-2.0),c1]
c1guess = 1;
C1 = @(varargin) @(c1) fsolve(C1inner(c1, varargin{:}), c1guess);
11 Comments
Taeksang Yoon
on 29 Dec 2021
Thank you for the answer, Walter!
If you don't mind, may I request more explanation about the code from defining vars?
I'm not sure about whether I understand your code.
If my explanation is wrong, please comment me.
V1 = [a1,a2,a3,a4]
I think this code is for defining the unknown parameter sets in equation.
V2 = {t,c10,rho_cat,V, R, T}
This code is for defining parameter sets in equation I set in the question.
V3 = setdiff(setdiff(symvar(C1f), V1), [V2{:}])
This code is for defining rest of parameters not included in V1 nor V2.
vars = {V3, V1, V2{:}}
Since we get all parameters in equation through defining V1 ~ V3, we make it as a one set. I think the {:} for V2 is to classify each parameter in V2 individually.
C1inner = matlabFunction(C1f, 'Vars', vars)
c1guess = 1;
C1 = @(varargin) @(c1) fsolve(C1inner(c1, varargin{:}), c1guess);
from here, I don't understand since 'C1inner' already has the form of 'C1' in my question
C1 = matlabFunction([C1], 'Vars',{[a1,a2,a3,a4],t,c10,rho_cat,V, R, T})
What is the difference between C1inner and C1?
What is the purpose for c1guess and C1 in your code?
Walter Roberson
on 29 Dec 2021
Your matlabFunction call defines two kinds of variables. The part you coded as [a1, a2, a3, a4] indicates that you want those four variables to be bundled into a single vector for the purpose of the call. The other ones are to be individually passed. This arrangement would be common for the case where the vector holds coefficients that need to be found and the others are parameters to the system.
The code I created needs to keep that structure for compatibility with your code.
Now, for any given set of values of the vector and parameters, the solution to the dsolve is the constant "c1" such that c1 is the root(s) of an equation, excluding 0. This effectively introduces a new variable into the system, the constant(s) to be found.
But we cannot assume that we know the name of the variable that dsolve will introduce. So we have to take the list of variables from the equation and remove from that the known variable names, and what is left is the name of the introduced variable.
We need the two kinds of known variables to be distinct so that we can construct the right calling sequence, and the calling sequence needs us to use a cell array to define it for matlabFunction purposes. But we need the variable names as lists to filter out the known names. So we end up needing to either define the variables twice, once in list form and once in cell form, or else we need to convert between list and cell as needed, under the restrictions that some of the functions that would normally be used to convert between forms are not permitted on symbolic vectors.
Walter Roberson
on 29 Dec 2021
Edited: Walter Roberson
on 29 Dec 2021
The C1 you had is in the form of a call to solve() with an additional set difference to remove 0. But matlabFunction cannot handle set differences, so we first have to strip off the set difference level.
Once we have stripped off the set difference we can hope that matlabFunction can process the solve that remains. But it cannot: it generates a function that calls solve() and which uses the undefined variable c1. solve() cannot work on numeric parameters, and matlabFunction cannot introduce a symbolic variable to solve over. matlabFunction has bugs with respect to "bound" variables that are to be solved for or which are variables of integration or limit() variables.
So we have to pull apart the solve() to get the expression. And then we need to generate a numeric root finder function.
Your C1 was intended to take a particular set of parameters and my C1 takes the same parameter list. For those particular parameters, the solution has to be found numerically as the c1 such that the expression becomes 0. C1inner is the function that has to be called by fsolve to try to find the root, with the fsolve solution being what should be returned by C1.
The code is not perfect. For any given set of parameters then C1 should be a set of values that satisfy the expression, with 0 removed from the set, whereas this code only tries to find a single solution and does not filter out 0 even.
Taeksang Yoon
on 29 Dec 2021
Edited: Taeksang Yoon
on 29 Dec 2021
Thanks for the explanation!
So, the reason for defining C1 twice is to define the function as two kinds of form (once in list form and once in cell form) to use it properly.
As you commented, the code is not perfect since my question is only part of my full code. What I'm trying is estimate the parameters {a1,a2,a3,a4} by using 'lsqcurvefit' and experiment data and finding the solution of differential equation is previous step before the data fitting.
syms c1(t) rho_cat V a1 a2 a3 a4 c10 R T t
Eq1 = diff(c1) == (-(rho_cat / V) * a1* exp(-a2 / (R * (T+273.15))) * (c1^a3)) / (1+a4*c1);
[C1] = dsolve(Eq1,c1(0) == c10);
C11 = children(children(C1,1),1);
C1f = lhs(C11)-rhs(C11);
V1 = [a1,a2,a3,a4];
V2 = {t,c10,rho_cat,V, R, T};
V3 = setdiff(setdiff(symvar(C1f), V1), [V2{:}]);
vars = {V3, V1, V2{:}};
C1inner = matlabFunction(C1f, 'Vars', vars)
C1inner = function_handle with value:
@(c1,in2,t,c10,rho_cat,V,R,T)[-((c10.*c10.^(-in2(:,3)).*(in2(:,3)-in2(:,4).*c10+in2(:,3).*in2(:,4).*c10-2.0))./(in2(:,3).*-3.0+in2(:,3).^2+2.0)+(in2(:,1).*rho_cat.*t.*exp((in2(:,2).*-2.0e+1)./(R.*5.463e+3+R.*T.*2.0e+1)))./V).*(in2(:,3).*-3.0+in2(:,3).^2+2.0)+c1.^(-in2(:,3)+1.0).*(in2(:,3)-in2(:,4).*c1+in2(:,3).*in2(:,4).*c1-2.0),c1]
c1guess = 1;
C1 = @(varargin) @(c1) fsolve(C1inner(c1, varargin{:}), c1guess)
C1 = function_handle with value:
@(varargin)@(c1)fsolve(C1inner(c1,varargin{:}),c1guess)
fsolve(C1inner(c1, varargin{:}), c1guess)
Unable to use curly braces to index into varargin.
When I separated 'fsolve' part and run it, it returns error message. Is it because it needs more expression as you commented?
Walter Roberson
on 29 Dec 2021
Using completely made-up values:
syms c1(t) rho_cat V a1 a2 a3 a4 c10 R T t
Eq1 = diff(c1) == (-(rho_cat / V) * a1* exp(-a2 / (R * (T+273.15))) * (c1^a3)) / (1+a4*c1);
[C1] = dsolve(Eq1,c1(0) == c10);
C11 = children(children(C1,1),1);
C1f = lhs(C11)-rhs(C11);
V1 = [a1,a2,a3,a4];
V2 = {t,c10,rho_cat,V, R, T};
V3 = setdiff(setdiff(symvar(C1f), V1), [V2{:}]);
vars = {V3, V1, V2{:}};
C1inner = matlabFunction(C1f, 'Vars', vars)
C1inner = function_handle with value:
@(c1,in2,t,c10,rho_cat,V,R,T)[-((c10.*c10.^(-in2(:,3)).*(in2(:,3)-in2(:,4).*c10+in2(:,3).*in2(:,4).*c10-2.0))./(in2(:,3).*-3.0+in2(:,3).^2+2.0)+(in2(:,1).*rho_cat.*t.*exp((in2(:,2).*-2.0e+1)./(R.*5.463e+3+R.*T.*2.0e+1)))./V).*(in2(:,3).*-3.0+in2(:,3).^2+2.0)+c1.^(-in2(:,3)+1.0).*(in2(:,3)-in2(:,4).*c1+in2(:,3).*in2(:,4).*c1-2.0),c1]
c1guess = 1;
C1 = @(varargin) fsolve(@(c1) C1inner(c1, varargin{:}), c1guess)
C1 = function_handle with value:
@(varargin)fsolve(@(c1)C1inner(c1,varargin{:}),c1guess)
V1_init = [.1 0 .3 .5];
t_init = 2;
c10_init = sqrt(3);
rho_cat_init = 1/37;
V_init = 7;
R_init = 8.314;
T_init = 285;
C1(V1_init, t_init, c10_init, rho_cat_init, V_init, R_init, T_init)
Warning: Trust-region-dogleg algorithm of FSOLVE cannot handle non-square systems; using Levenberg-Marquardt algorithm instead.
No solution found.
fsolve stopped because the last step was ineffective. However, the vector of function
values is not near zero, as measured by the value of the function tolerance.
ans = 1.3350
Taeksang Yoon
on 30 Dec 2021
Thanks for the additional comment!
So, it looks like I need the value of all other parameters to solve the equation.
If so, I have to find other way to solve my situation since my final purpose is to estimate parameter a1 ~ a4 classified as V1 in the code.
If it does not bother you, may I ask another question? If you want it, I will post as a new question.
I'm doing parameter estimation (a1 to a4) of my reaction equation.
My code is as follows. I fixed my code as you commented to me.
clear, clc
format long
syms c1(t) rho_cat V a1 a2 a3 a4 c10 R T t
Eq1 = diff(c1) == (-(rho_cat / V) * a1* exp(-a2 / (R * (T+273.15))) * (c1^a3)) / (1+a4*c1);
[C1] = dsolve(Eq1,c1(0) == c10)
C1 =
C11 = children(children(C1,1),1);
C1f = lhs(C11)-rhs(C11);
V1 = [a1,a2,a3,a4];
V2 = {t,c10,rho_cat,V, R, T};
V3 = setdiff(setdiff(symvar(C1f), V1), [V2{:}]);
vars = {V3, V1, V2{:}};
C1inner = matlabFunction(C1f, 'Vars', vars)
C1inner = function_handle with value:
@(c1,in2,t,c10,rho_cat,V,R,T)[-((c10.*c10.^(-in2(:,3)).*(in2(:,3)-in2(:,4).*c10+in2(:,3).*in2(:,4).*c10-2.0))./(in2(:,3).*-3.0+in2(:,3).^2+2.0)+(in2(:,1).*rho_cat.*t.*exp((in2(:,2).*-2.0e+1)./(R.*5.463e+3+R.*T.*2.0e+1)))./V).*(in2(:,3).*-3.0+in2(:,3).^2+2.0)+c1.^(-in2(:,3)+1.0).*(in2(:,3)-in2(:,4).*c1+in2(:,3).*in2(:,4).*c1-2.0),c1]
c1guess = 1;
C1 = @(varargin) @(c1) fsolve(C1inner(c1, varargin{:}), c1guess)
C1 = function_handle with value:
@(varargin)@(c1)fsolve(C1inner(c1,varargin{:}),c1guess)
%% defining parameter
P = [80;120;100;80;120;80;120;100;80;120;80;80;100;120;100]; % Pressure, bar
T = [80;80;80;80;80;100;100;100;100;100;120;120;120;120;120]; % Temperature, C
t=0.1; %reactor length 0.1m
R = 8.314 ; % gas constant J/mol*K
S = 0.001808; % cross-sectional area of reactor, m2
Q_h = [9.5290;7.9509;5.7226;12.7060;10.6018;9.7972;8.1296;5.8656;13.0635;10.8401;10.0653;6.7114;6.0086;5.5401;12.0155]; % Volumetric flowrate, L/h
Q = Q_h/1000; % Conversion of unit, L/h to m3/h
rho_cat = 500; % catalyst density, kg/m3
c10 = 100.*P./(R.*(T+273.15)); % first concentration of H2
X = [6.6;10.2;12.7;5.6;9.7;14.2;16.5;18;13.5;15.1;22;30.2;33.1;44.4;29.6]; % Conversion
c = c10.*(1-X./100); % last concentration of H2
V=Q/S; % linear velocity of fluid
%% paramter estimation a1 ~ a4
in=rand(1,4); %random value for initial value of a1 ~ a4
myfunct=@(in1,T) C1(in1,t,c10,rho_cat,V,R,T)
myfunct = function_handle with value:
@(in1,T)C1(in1,t,c10,rho_cat,V,R,T)
options = optimset('Display','iter','TolX', 1e-15, 'TolFun', 1e-15, 'MaxFunEvals', 4000, 'MaxIter', 4000,'Algorithm','Levenberg-Marquardt')
options = struct with fields:
Display: 'iter'
MaxFunEvals: 4000
MaxIter: 4000
TolFun: 1.000000000000000e-15
TolX: 1.000000000000000e-15
FunValCheck: []
OutputFcn: []
PlotFcns: []
ActiveConstrTol: []
Algorithm: 'levenberg-marquardt'
AlwaysHonorConstraints: []
DerivativeCheck: []
Diagnostics: []
DiffMaxChange: []
DiffMinChange: []
FinDiffRelStep: []
FinDiffType: []
GoalsExactAchieve: []
GradConstr: []
GradObj: []
HessFcn: []
Hessian: []
HessMult: []
HessPattern: []
HessUpdate: []
InitBarrierParam: []
InitTrustRegionRadius: []
Jacobian: []
JacobMult: []
JacobPattern: []
LargeScale: []
MaxNodes: []
MaxPCGIter: []
MaxProjCGIter: []
MaxSQPIter: []
MaxTime: []
MeritFunction: []
MinAbsMax: []
NoStopIfFlatInfeas: []
ObjectiveLimit: []
PhaseOneTotalScaling: []
Preconditioner: []
PrecondBandWidth: []
RelLineSrchBnd: []
RelLineSrchBndDuration: []
ScaleProblem: []
SubproblemAlgorithm: []
TolCon: []
TolConSQP: []
TolGradCon: []
TolPCG: []
TolProjCG: []
TolProjCGAbs: []
TypicalX: []
UseParallel: []
[in1,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(myfunct,in,T,c,[],[],options)
Error using lsqcurvefit (line 286)
Function value and YDATA sizes are not equal.
Function value and YDATA sizes are not equal.
in1
When I run lsqcurvefit to estimate a1 ~ a4, the matlab sends the error message as above.
I think the reason of difference size between function and YDATA is the presence of c1.
Is there a way to check the size of function and YDATA? I tried to check it by fixing lsqcurvefit function, but I'm stucked at here.
Walter Roberson
on 30 Dec 2021
What Maple says is
c1(t) = RootOf(Z^(-a3+2)*exp(20*a2/R/(20*T+5463))*V*a3*a4-c10^(-a3+2)*exp(20*
a2/R/(20*T+5463))*V*a3*a4-Z^(-a3+2)*exp(20*a2/R/(20*T+5463))*V*a4+c10^(-a3+2)*
exp(20*a2/R/(20*T+5463))*V*a4-a1*a3^2*rho_cat*t+Z^(1-a3)*exp(20*a2/R/(20*T+
5463))*V*a3-c10^(1-a3)*exp(20*a2/R/(20*T+5463))*V*a3+3*a1*a3*rho_cat*t-2*V*Z^(
1-a3)*exp(20*a2/R/(20*T+5463))+2*V*c10^(1-a3)*exp(20*a2/R/(20*T+5463))-2*t*a1*
rho_cat, Z)
Here, RootOf(expression in Z, Z) means the set of Z such that the expression becomes 0.
Notice that is an expression in t: for every different t you have to solve it again.
It is difficult to tell with your code what t corresponds to.
The most natural possibility for a curve-fitting situation is if t were your independent variable, and your ydata was the measured response. If that is what you were doing, then t would be T, the variable you are passing as XData. But your T does not contain unique values; to use lsqcurvefit() your XData must contain unique values.
Another possibility is that t is a final time after the reaction. That would fit in with your t=0.1 line. But then what is your independent variable? It cannot be P or T as those both have duplicates. It cannot be (P,T) pairs, as the second and 5th pair are the same.
As I look at what you have, it does not look to me as if you have a Y vs X situation: it looks to me as if you have scattered multidimensional data. In such a case you should not be using lsqcurvefit() as that is not designed for more than 3 dimensions.
What should you do instead of lsqcurvefit()? Probably you should construct a sum-of-squared-residues and fmincon() or fminsearch() the parameters. Each residue would be a predicted concentration minus the corresponding measured concentration. The predicted concentration might require fsolve(). Hmmm, maybe the whole thing could be written as fsolve() of a series of equations in which each row used corresponding parameters of the known variables, along with the parameters to be estimated, and subtract off the known concentration... fsolve() would then be trying to find that parameter values that matched concentrations. But since the system would be overdetermined, instead fsolve() would end up trying to minimize the error, which is fine for your purposes.
Taeksang Yoon
on 30 Dec 2021
What a complex root !
t is the reactor length fixed with 0.1m. However, as the velocity is different for each experiment case, the reaction time varies also.
As you told me my situation is with scattered multidimensional data. However, even it has multidimensional data, the lsqcurvefit could estimate the parameter well in the case of different c1 equation(though it sometimes returns imaginary number or wrong answer). I modified the code based on this link. My code is as follows, you may check the code and comment to me.
syms c1(t) rho_cat V a1 a2 a3 c10 R T t
Eq1 = diff(c1) == -(rho_cat / V) * a1* exp(-a2 / (R * (T+273.15))) * c1^a3
Eq1(t) =
[C1] = dsolve(Eq1,c1(0) == c10)
C1 =
C1 = matlabFunction([C1], 'Vars',{[a1,a2,a3],t,c10,rho_cat,V, R, T})
C1 = function_handle with value:
@(in1,t,c10,rho_cat,V,R,T)(-(in1(:,3)-1.0).*(c10./(c10.^in1(:,3)-in1(:,3).*c10.^in1(:,3))-(in1(:,1).*rho_cat.*t.*exp((in1(:,2).*-2.0e+1)./(R.*5.463e+3+R.*T.*2.0e+1)))./V)).^(-1.0./(in1(:,3)-1.0))
%% defining parameter
P = [80;120;100;80;120;80;120;100;80;120;80;80;100;120;100]; % pressure, bar
T = [80;80;80;80;80;100;100;100;100;100;120;120;120;120;120]; % temperature, C
t=0.1; %reactor length 0.1m
R = 8.314 ; % gas constant J/mol*K
S = 0.001808; % cross-sectional area, m2
Q_h = [9.5290;7.9509;5.7226;12.7060;10.6018;9.7972;8.1296;5.8656;13.0635;10.8401;10.0653;6.7114;6.0086;5.5401;12.0155]; % volumetric flowrate, L/h
Q = Q_h/1000; % unit conversion, L/h to m3/h
rho_cat = 500; % catalyst density, kg/m3
c10 = 100.*P./(R.*(T+273.15)); % First concentration of H2
X = [6.6;10.2;12.7;5.6;9.7;14.2;16.5;18;13.5;15.1;22;30.2;33.1;44.4;29.6]; % Conversion X, %
c = c10.*(1-X./100); % last concentration of H2
V=Q/S; % linear velocity m/h
%% parameter estimation
in=rand(1,3); % random initial value of a1 to a3
myfunct=@(in1,T) C1(in1,t,c10,rho_cat,V,R,T)
myfunct = function_handle with value:
@(in1,T)C1(in1,t,c10,rho_cat,V,R,T)
options = optimset('Display','iter','TolX', 1e-15, 'TolFun', 1e-15, 'MaxFunEvals', 4000, 'MaxIter', 4000,'Algorithm','Levenberg-Marquardt')
options = struct with fields:
Display: 'iter'
MaxFunEvals: 4000
MaxIter: 4000
TolFun: 1.0000e-15
TolX: 1.0000e-15
FunValCheck: []
OutputFcn: []
PlotFcns: []
ActiveConstrTol: []
Algorithm: 'levenberg-marquardt'
AlwaysHonorConstraints: []
DerivativeCheck: []
Diagnostics: []
DiffMaxChange: []
DiffMinChange: []
FinDiffRelStep: []
FinDiffType: []
GoalsExactAchieve: []
GradConstr: []
GradObj: []
HessFcn: []
Hessian: []
HessMult: []
HessPattern: []
HessUpdate: []
InitBarrierParam: []
InitTrustRegionRadius: []
Jacobian: []
JacobMult: []
JacobPattern: []
LargeScale: []
MaxNodes: []
MaxPCGIter: []
MaxProjCGIter: []
MaxSQPIter: []
MaxTime: []
MeritFunction: []
MinAbsMax: []
NoStopIfFlatInfeas: []
ObjectiveLimit: []
PhaseOneTotalScaling: []
Preconditioner: []
PrecondBandWidth: []
RelLineSrchBnd: []
RelLineSrchBndDuration: []
ScaleProblem: []
SubproblemAlgorithm: []
TolCon: []
TolConSQP: []
TolGradCon: []
TolPCG: []
TolProjCG: []
TolProjCGAbs: []
TypicalX: []
UseParallel: []
[in1,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(myfunct,in,T,c,[],[],options);
First-Order Norm of
Iteration Func-count Residual optimality Lambda step
0 4 56.3783 222 0.01
1 11 38.4363 559 10 0.323216
2 15 1.74014 47.8 1 0.0663336
3 19 1.30998 1.27 0.1 0.0145596
4 23 1.30804 0.475 0.01 0.0691817
5 27 1.30774 0.134 0.001 0.0381017
6 31 1.30774 0.000952 0.0001 0.0328154
7 35 1.30771 0.000292 1e-05 0.326987
8 39 1.3075 0.00267 1e-06 3.27088
9 43 1.30536 0.0323 1e-07 32.6869
10 47 1.28449 0.785 1e-08 323.766
11 52 1.26401 0.827 1e-07 318.798
12 56 1.24387 0.739 1e-08 315.437
13 61 1.22413 0.657 1e-07 312.369
14 65 1.20477 0.583 1e-08 309.325
15 70 1.18579 0.517 1e-07 306.294
16 74 1.16719 0.46 1e-08 303.274
17 79 1.14896 0.409 1e-07 300.267
18 83 1.13109 0.364 1e-08 297.274
19 88 1.11359 0.325 1e-07 294.297
20 92 1.09644 0.29 1e-08 291.337
21 97 1.07963 0.259 1e-07 288.395
22 101 1.06317 0.232 1e-08 285.473
23 106 1.04705 0.207 1e-07 282.572
24 110 1.03125 0.186 1e-08 279.693
25 115 1.01578 0.167 1e-07 276.836
26 119 1.00063 0.15 1e-08 274.003
27 124 0.985788 0.134 1e-07 271.195
28 128 0.971253 0.121 1e-08 268.412
29 133 0.957016 0.109 1e-07 265.655
30 137 0.943073 0.0981 1e-08 262.925
31 142 0.929418 0.0885 1e-07 260.221
32 146 0.916044 0.0799 1e-08 257.545
33 151 0.902945 0.0722 1e-07 254.896
34 155 0.890117 0.0653 1e-08 252.275
35 160 0.877552 0.0591 1e-07 249.683
36 164 0.865246 0.0536 1e-08 247.12
37 169 0.853193 0.0486 1e-07 244.585
38 173 0.841387 0.0441 1e-08 242.079
39 178 0.829823 0.0401 1e-07 239.602
40 182 0.818496 0.0364 1e-08 237.154
41 187 0.807399 0.0332 1e-07 234.735
42 191 0.796529 0.0302 1e-08 232.345
43 196 0.78588 0.0276 1e-07 229.984
44 200 0.775447 0.0251 1e-08 227.652
45 204 0.772985 0.666 1e-09 2096.68
46 209 0.709465 0.677 1e-08 1739.8
47 213 0.642718 0.647 1e-09 1595.32
48 218 0.572688 0.567 1e-08 1481.44
49 222 0.51487 0.501 1e-09 1388.06
50 227 0.46501 0.434 1e-08 1300.02
51 231 0.423416 0.376 1e-09 1219.03
52 236 0.388217 0.325 1e-08 1143.39
53 240 0.358425 0.282 1e-09 1073.33
54 245 0.333013 0.244 1e-08 1008.43
55 249 0.311206 0.211 1e-09 948.461
56 254 0.292366 0.184 1e-08 893.082
57 258 0.275977 0.16 1e-09 842.505
58 263 0.261694 0.141 1e-08 794.683
59 267 0.249087 0.123 1e-09 751.072
60 272 0.23797 0.108 1e-08 710.74
61 276 0.228102 0.0953 1e-09 673.38
62 281 0.21931 0.0843 1e-08 638.72
63 285 0.211446 0.0747 1e-09 606.507
64 290 0.204392 0.0663 1e-08 576.51
65 294 0.198043 0.059 1e-09 548.523
66 299 0.192317 0.0526 1e-08 522.35
67 303 0.187139 0.0469 1e-09 497.824
68 308 0.182449 0.0419 1e-08 474.786
69 312 0.178193 0.0374 1e-09 453.1
70 317 0.174326 0.0334 1e-08 432.643
71 321 0.170807 0.0298 1e-09 413.308
72 326 0.167602 0.0266 1e-08 395.002
73 330 0.16468 0.0238 1e-09 377.64
74 335 0.162014 0.0212 1e-08 361.158
75 339 0.159579 0.0189 1e-09 345.495
76 344 0.157355 0.0168 1e-08 330.604
77 348 0.15532 0.0149 1e-09 316.442
78 353 0.153459 0.0133 1e-08 302.976
79 357 0.151754 0.0118 1e-09 290.178
80 362 0.150192 0.0105 1e-08 278.022
81 366 0.148758 0.00934 1e-09 266.484
82 371 0.147442 0.00831 1e-08 255.542
83 375 0.146231 0.00739 1e-09 245.177
84 380 0.145117 0.00659 1e-08 235.365
85 384 0.144089 0.00587 1e-09 226.086
86 389 0.14314 0.00525 1e-08 217.318
87 393 0.142263 0.00469 1e-09 209.036
88 398 0.141451 0.0042 1e-08 201.219
89 402 0.140697 0.00377 1e-09 193.843
90 407 0.139997 0.00339 1e-08 186.884
91 411 0.139346 0.00306 1e-09 180.321
92 416 0.138739 0.00276 1e-08 174.13
93 420 0.138172 0.0025 1e-09 168.287
94 425 0.137641 0.00227 1e-08 162.771
95 429 0.137145 0.00206 1e-09 157.564
96 433 0.137017 0.12 1e-10 1339.13
97 438 0.132375 0.0648 1e-09 1047.22
98 442 0.130166 0.0355 1e-10 861.963
99 447 0.128894 0.0211 1e-09 732.681
100 451 0.128015 0.0137 1e-10 638.685
101 456 0.127347 0.00952 1e-09 567.44
102 460 0.126813 0.00692 1e-10 511.574
103 465 0.126373 0.00522 1e-09 466.541
104 469 0.126001 0.00405 1e-10 429.423
105 474 0.125682 0.00322 1e-09 398.264
106 478 0.125404 0.00262 1e-10 371.696
107 483 0.12516 0.00216 1e-09 348.751
108 487 0.124943 0.00181 1e-10 328.712
109 492 0.124749 0.00153 1e-09 311.041
110 496 0.124574 0.00131 1e-10 295.33
111 501 0.124416 0.00113 1e-09 281.257
112 505 0.124271 0.000989 1e-10 268.569
113 509 0.124061 0.053 1e-11 2129.69
114 514 0.122942 0.0232 1e-10 1589.8
115 518 0.122503 0.0124 1e-11 1277.84
116 523 0.122243 0.00746 1e-10 1068.99
117 527 0.122061 0.0049 1e-11 918.513
118 532 0.121924 0.00341 1e-10 804.392
119 536 0.121818 0.00248 1e-11 714.566
120 541 0.121733 0.00186 1e-10 641.835
121 545 0.121663 0.00143 1e-11 581.625
122 550 0.121605 0.00113 1e-10 530.889
123 554 0.121557 0.000906 1e-11 487.493
124 559 0.121515 0.000737 1e-10 449.947
125 563 0.12148 0.000607 1e-11 417.106
126 567 0.121436 0.0214 1e-12 2595.39
127 572 0.121259 0.00768 1e-11 1681.1
128 576 0.121214 0.00323 1e-12 1158.45
129 581 0.121195 0.00152 1e-11 829.767
130 585 0.121186 0.000773 1e-12 609.739
131 589 0.121179 0.0038 1e-13 1394.91
132 593 0.121174 0.000529 1e-14 527.224
133 597 0.121174 7.18e-06 1e-15 47.019
134 601 0.121174 1.46e-08 1e-16 2.14934
135 605 0.121174 1.44e-08 1e-17 0.00878556
136 614 0.121174 1.49e-09 1e-12 0.205549
137 618 0.121174 5.47e-09 1e-13 0.088525
138 627 0.121174 1.7e-08 1e-08 6.85128e-05
Local minimum possible.
lsqcurvefit stopped because the final change in the sum of squares relative to
its initial value is less than the value of the function tolerance.
in1 % show a1, a2, a3
in1 =
1.0e+04 *
3.1568 - 0.0000i 4.7740 - 0.0000i 0.0002 - 0.0000i
%% validation data
P_V = [100;80;100;100;80;100;100]; % Pressure bar
T_V = [80;80;80;100;100;100;120]; %temperature C
c10_V = 100.*P_V./(R.*(T_V+273.15)); % first concentration of H2,
Q_hV = [8.5822;6.3539;11.4435;8.7967;6.5326;11.7295;9.0112]; % Volumetric flowrate, L/h
Q_V = Q_hV/1000; % unit conversion L/h to m3/h
X_V = [7.8;9.7;7.3;15;16.3;14.1;31.2]; % Conversion
c_V = c10_V.*(1-X_V./100); % Final concentration of H2
V_V = Q_V/S; %linear velocity
c_Calc = 1./(-(real(in1(3)) - 1).*(c10_V./(c10_V.^real(in1(3)) - real(in1(3)).*c10_V.^real(in1(3))) - (real(in1(1)).*rho_cat.*t.*exp(-(20.*real(in1(2)))./(5463.*R + 20.*R.*T_V)))./V_V)).^(1./(real(in1(3)) - 1)); % calculated final concentration of H2 data
X_Calc = 100-100.*c_Calc./c10_V; % Calculated conversion
Rsq = 1 - sum((X_V - X_Calc).^2) / sum((X_V - mean(X_V)).^2) % Rsquare calculation
Rsq = 0.9576
As you can see above I could estimate parameter with lsqcurvefit, so I tried it with different shape of equation c1. However, as the equation becomes complex, the situation becomes much more complex to solve than I thought. when I try it with below equation, the matlab returns error as follows. That's why I asked how to check the size of function and YDATA.
Error using lsqcurvefit (line 286)
Function value and YDATA sizes are not equal.
As you commented in the last paragraph, maybe I have to re-write the whole code with different function like fsolve for finding {a1,a2,a3,a4}. Maybe the problem cannot be solved with this equation. I'll try my best.
syms c1(t) rho_cat V a1 a2 a3 a4 c10 R T t
Eq1 = diff(c1) == (-(rho_cat / V) * a1* exp(-a2 / (R * (T+273.15))) * (c1^a3)) / (1+a4*c1)
Eq1(t) =
Torsten
on 30 Dec 2021
I see that you already simplified the ODE for c1. My guess for the error message is that your function "myfunct" only returns one value for c1. Use "arrayfun" to return a data vector of the length of c calculated with the actual parameters "in1" .
My advice would be to use a numerical integrator (e.g. ODE45) to determine c1 depending on P, T, V and c10. This will make things much easier.
Taeksang Yoon
on 31 Dec 2021
Thanks for comment, Torsten! I will try "arrayfun".
I tried to use numerical integrator like ode45 or ode23 also. However, when I searched about them in MATLAB help, the code seems like it can only solve problem with one variable, t. If so, I can use only small part of my data since they are multi-dimensional (P,T,V etc.)
Maybe I misunderstood ode45. If it can solve it with multivariables, can you suggest the example of it? It would be grateful if you help me.
Taeksang Yoon
on 6 Jan 2022
@Torsten I tried to use "arrayfun" as you recommended however, I have no idea of combining "arrayfun" and "lsqcurvefit". Do you have any advice for it?
More Answers (1)
Torsten
on 31 Dec 2021
Edited: Torsten
on 31 Dec 2021
function main
%% defining parameter
P = [80;120;100;80;120;80;120;100;80;120;80;80;100;120;100]; % Pressure, bar
T = [80;80;80;80;80;100;100;100;100;100;120;120;120;120;120]; % Temperature, C
t = 0.1; %reactor length 0.1m
R = 8.314 ; % gas constant J/mol*K
S = 0.001808; % cross-sectional area of reactor, m2
Q_h = [9.5290;7.9509;5.7226;12.7060;10.6018;9.7972;8.1296;5.8656;13.0635;10.8401;10.0653;6.7114;6.0086;5.5401;12.0155]; % Volumetric flowrate, L/h
Q = Q_h/1000; % Conversion of unit, L/h to m3/h
rho_cat = 500; % catalyst density, kg/m3
c10 = 100.*P./(R.*(T+273.15)); % first concentration of H2
X = [6.6;10.2;12.7;5.6;9.7;14.2;16.5;18;13.5;15.1;22;30.2;33.1;44.4;29.6]; % Conversion
c = c10.*(1-X./100); % last concentration of H2
V = Q/S; % linear velocity of fluid
%% paramter estimation a1 ~ a4
in = rand(1,4); %random value for initial value of a1 ~ a4
options = optimset('Display','iter','TolX', 1e-15, 'TolFun', 1e-15, 'MaxFunEvals', 4000, 'MaxIter', 4000,'Algorithm','Levenberg-Marquardt')
[in1,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqnonlin(@(in1)myfunct(in1,c,t,c10,rho_cat,P,V,R,T),in,[],[],options)
in1
end
function res = myfunct(in1,c,t,c10,rho_cat,P,V,R,T)
tspan = [0 t];
for i = 1:numel(P)
func = @(tt,y)(-(rho_cat / V(i)) * in1(1)* exp(-in1(2) / (R * (T(i)+273.15))) * ...
(y(1)^in1(3))) / (1+in1(4)*y(1));
[T,sol] = ode15s(func,tspan,c10(i));
res(i) = c(i)-sol(end,1);
end
end
11 Comments
Taeksang Yoon
on 3 Jan 2022
Thanks for the answer!
Happy new year. It's been a year to see you again.
I have additional question about your comment. I ran the code in the MATLAB but it does not return the result. Does your code need additional code to run properly?
Taeksang Yoon
on 10 Jan 2022
When I run the code, it usually returns the error message with 'Index exceeds the number of array elements'.
If this numerical solution does not yield a result, the above symbolic variant won't do either. So my guess.
->Okay, I got it.
Torsten
on 10 Jan 2022
Edited: Torsten
on 10 Jan 2022
Replace
[T,sol] = ode15s(func,tspan,c10(i));
by
[~,sol] = ode15s(func,tspan,c10(i));
to avoid conflicts with the temperature array T.
Of course, the arrays for P, T, V, c and c10 must be of the same size for the code to run properly.
If the solver does not converge is a different problem.
Taeksang Yoon
on 11 Jan 2022
Replace
[T,sol] = ode15s(func,tspan,c10(i));
[~,sol] = ode15s(func,tspan,c10(i));
This make the code run without error, Thank you!
Taeksang Yoon
on 14 Jan 2022
Hello, If it doesn't bother you, may I ask additional question?
I got some another question on your code.
function res = myfunct(in1,c,t,c10,rho_cat,P,V,R,T)
tspan = [0 t];
for i = 1:numel(P)
func = @(tt,y)(-(rho_cat / V(i)) * in1(1)* exp(-in1(2) / (R * (T(i)+273.15))) * ...
(y(1)^in1(3))) / (1+in1(4)*y(1));
[T,sol] = ode15s(func,tspan,c10(i));
res(i) = c(i)-sol(end,1);
end
In this code, it looks like the sol(end,1) is indicating a specific value in matrix which does not change for every iteration. Did I understand correctly?
If so, Why the sol(end,1) is selected among the matrix sol?
The reason I'm asking this is because I'm on the verification of the parameter I estimated with other data. I made the code for validation based on the function you made using ode15s to get the calculated final value of reactants. The file in1.m is the code that you made and the other code is for validating parameter estimated in in1.m. When I run the code, I got 11x7 double matrix for sol.
I wonder which is the right value of c, the final concentration of H2.
Always thank you.
Taeksang Yoon
on 17 Jan 2022
Thanks for the reply.
If sol(end,1) is the solution of the differential equation, do I have to write the code for validation as a for expression? I want to check the estimation result also fits to other experiment data.
Torsten
on 17 Jan 2022
Yes, set P_new, T_new, V_new, c_new and c10_new for the other data and call ode15s with the parameters "in1" obtained from the optimization:
in = rand(1,4); %random value for initial value of a1 ~ a4
options = optimset('Display','iter','TolX', 1e-15, 'TolFun', 1e-15, 'MaxFunEvals', 4000, 'MaxIter', 4000,'Algorithm','Levenberg-Marquardt')
[in1,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqnonlin(@(in1)myfunct(in1,c,t,c10,rho_cat,P,V,R,T),in,[],[],options)
P_new = ...;
T_new = ...;
V_new = ...;
c_new = ...;
c10_new = ...;
res = myfunct(in1,c_new,t,c10_new,rho_cat,P_new,V_new,R,T_new)
Taeksang Yoon
on 20 Jan 2022
Thank you!
I made a code with for expression and the code can show the validation result properly.
Your comments were really helpful to solve the problem!
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)