How can I solve this differential equation?
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.
1 Comment
Accepted Answer
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)
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)
c1guess = 1;
C1 = @(varargin) @(c1) fsolve(C1inner(c1, varargin{:}), c1guess);
11 Comments
Taeksang Yoon
on 29 Dec 2021
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
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)
c1guess = 1;
C1 = @(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)
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)
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)
c1guess = 1;
C1 = @(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)
options = optimset('Display','iter','TolX', 1e-15, 'TolFun', 1e-15, 'MaxFunEvals', 4000, 'MaxIter', 4000,'Algorithm','Levenberg-Marquardt')
[in1,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(myfunct,in,T,c,[],[],options)
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.
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
[C1] = dsolve(Eq1,c1(0) == c10)
C1 = matlabFunction([C1], 'Vars',{[a1,a2,a3],t,c10,rho_cat,V, R, T})
%% 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)
options = optimset('Display','iter','TolX', 1e-15, 'TolFun', 1e-15, 'MaxFunEvals', 4000, 'MaxIter', 4000,'Algorithm','Levenberg-Marquardt')
[in1,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(myfunct,in,T,c,[],[],options);
in1 % show a1, a2, a3
%% 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
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)
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
Taeksang Yoon
on 6 Jan 2022
More Answers (1)
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
Walter Roberson
on 3 Jan 2022
function in1 = main
Taeksang Yoon
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
Taeksang Yoon
on 14 Jan 2022
Taeksang Yoon
on 17 Jan 2022
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
Categories
Find more on Conversion Between Symbolic and Numeric in Help Center and File Exchange
on 28 Dec 2021
on 20 Jan 2022
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
