How to solve coupled odes with two time dependent variables with ode45?

30 次查看（过去 30 天）
Ari Dillep Sai
Ari Dillep Sai 2022-5-31
编辑： Torsten 2022-5-31
I am having two coupled odes in which there are two dependent variables are present. I tried solving them with the below but I am getting error. Can someone help me in solving the below equations?
%% equations which I need to solve are:
%% dc/dt = -((1-eps)*rhop/eps)*(dq/dt);
%% dq/dt = Kl*((qm*Keq*c*R*T/(1+(Keq*c*R*T)^n)^(1/n)) - q);
eps = 0.43;
rhop = 1228.5;
qm = 5.09;
R = 8.314;
T = 301;
Kl = 0.226;
n = 0.429;
K0 = 4.31e-9; % in pascal-1
delh = -29380; % heat of adsorption in j/mol
Keq = K0*exp(-delh/(R*T));
t = 0:1:100;
y0= zeros(2,1);
[tsol,ysol] = ode45(@(t,y) odfun(t,y), t, y0);
plot(tsol,ysol(:,1))
function dy = odfun(t,y,Kl,qm,Keq,R,T,n,rhop,eps)
c = y(1);
q = y(2);
dy(1) = Kl*((qm*Keq*c*R*T/(1+(Keq*c*R*T)^n)^(1/n))-q);
dy(2) = -((1-eps)*rhop/eps)*dy(1);
end
The error which I am getting when I am running this code is as follows:
Not enough input arguments.
Error in odcase>odfun (line 21)
dy(1) = Kl*((qm*Keq*c*R*T/(1+(Keq*c*R*T)^n)^(1/n))-q);
Error in odcase>@(t,y)odfun(t,y) (line 14)
[tsol,ysol] = ode45(@(t,y) odfun(t,y), t, y0);
Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.
Error in ode45 (line 115)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
Error in odcase (line 14)
[tsol,ysol] = ode45(@(t,y) odfun(t,y), t, y0);
Someone please let me know what are the mistakes in my code.

请先登录，再回答此问题。

采纳的回答

Torsten
Torsten 2022-5-31
编辑：Torsten 2022-5-31
function dy = odfun(t,y,Kl,qm,Keq,R,T,n,rhop,eps)
q = y(1);
c = y(2);
dy(1) = Kl*((qm*Keq*c*R*T/(1+(Keq*c*R*T)^n)^(1/n))-q);
dy(2) = -((1-eps)*rhop/eps)*dy(1);
end
instead of
function dy = odfun(t,y,Kl,qm,Keq,R,T,n,rhop,eps)
c = y(1);
q = y(2);
dy(1) = Kl*((qm*Keq*c*R*T/(1+(Keq*c*R*T)^n)^(1/n))-q);
dy(2) = -((1-eps)*rhop/eps)*dy(1);
end

更多回答（1 个）

Sam Chak
Sam Chak 2022-5-31
Hi @Ari Dillep Sai
Not sure what went wrong and why it is unstable. If you are absolutely sure that the absorption dynamics is stable (converging to a steady-state value), then the ODEs must be incorrect. Please check all parameters and the signs. Sometimes, a single change of sign can make a huge difference. For example, as a test, I simply added a minus sign '–' in front of K1 on the right-hand side of dydt(1), and the system becomes stable.
Please countercheck the ODEs against various textbooks and journal papers. If you only rely on a single reference and there is a misprint, then you know what happens...
function dydt = odefcn(t, y)
dydt = zeros(2,1);
ep = 0.43;
rhop = 1228.5;
qm = 5.09;
R = 8.314;
T = 301;
Kl = 0.226;
n = 0.429;
K0 = 4.31e-9; % in pascal-1
delh = -29380; % heat of adsorption in j/mol
Keq = K0*exp(-delh/(R*T));
c = y(1);
q = y(2);
dydt(1) = -Kl*( ( qm*Keq*c*R*T/(1 + (Keq*c*R*T)^n)^(1/n) ) - q );
dydt(2) = -((1 - ep)*rhop/ep)*dydt(1);
end
tspan = [0 0.1];
init = [1; 0];
[t, y] = ode45(@odefcn, tspan, init);
plot(t, y(:, 1), 'linewidth', 1.5)
  3 个评论
显示 1更早的评论
Ari Dillep Sai
Ari Dillep Sai 2022-5-31
Yeah, those are same equations. But I am not getting solution for those equations using method of lines. So first I thought of solving them by neglecting the spatial variations. So if I assume that there is no need of boundary conditions as they are varying only with z-axis.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by