How can I solve this differential equation?

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.

[Torsten]

If "dsolve" gives a solution for c1 in implicit form, then usually "solve" cannot solve this implicit solution for c1 explicitly.

Answer 1

[Walter Roberson]

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); 

[Taeksang Yoon]

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]

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]

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]

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]

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]

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.

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]

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]

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.

https://kr.mathworks.com/matlabcentral/answers/521824-fitting-parameters-in-ode-for-a-kinetic-reaction?s_tid=srchtitle

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]

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]

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.)

[t,y] = ode45(odefun,tspan,y0,options)

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]

@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?

Answer 2

[Torsten]

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

[Taeksang Yoon]

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?

[Walter Roberson]

function in1 = main

[Torsten]

Does lsqnonlin not end properly ?

in1, the result, should be printed at the end of the optimization.

If this numerical solution does not yield a result, the above symbolic variant won't do either. So my guess.

[Taeksang Yoon]

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]

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]

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]

@Torsten

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.

[Torsten]

For each value i, sol(end,1) is the solution of your differential equation at the endpoint of tspan,

i.e. at t = 0.1.

[Taeksang Yoon]

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]

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]

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!