How can I solve this differential equation?

4 views (last 30 days)
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
Torsten
Torsten on 28 Dec 2021
Edited: Torsten on 28 Dec 2021
If "dsolve" gives a solution for c1 in implicit form, then usually "solve" cannot solve this implicit solution for c1 explicitly.

Sign in to comment.

Accepted Answer

Walter Roberson
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
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.)
[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
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?

Sign in to comment.

More Answers (1)

Torsten
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
Torsten
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
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!

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!