Solve a system of coupled differential equations with ODE23s

19 次查看（过去 30 天）
Thomas Stone-Wigg
Thomas Stone-Wigg 2024-1-11
评论： William Rose 2024-1-12
Hi,
Im trying to solve the system of coupled differential equations (TMB, CMB and EB at the bottom of the code) with ODE23s but I'm unsure of the correct method to set up the equations to pass them to the solver.
This snippet is just to show how they are coupled:
dPdt(1:n) = ...* dTdt(1:n);
dy1dt(1:n) = ...*dPdt(1:n) ... *dTdt(1:n);
dTdt(1:n) = ...*dPdt(1:n);
I've attached the three equations and their non-dimensional and discretized forms for reference
At the moment im not too worried about the efficiency of the code but any extra advice would be much appreciated
Many thanks
Tom
clear,clc
%% MAIN
% Physical Parameters
L = 1; % Column length m
r = 0.25; % Bed radius m
t0 = 0; % Initial Time
tf = 500; % Final time
dt = 0.1; % Time step
t = t0:dt:tf; % Time vector
dz = 0.1; % Mesh size
z = 0:dz:L; % Mesh generation
n = numel(z);
%%%%%%% blended gas parameters %%%%%%%%%%%
%%% calc u input
TPF = [1 300 6e5 1.75e5 1e5]; % Feed: Velocity (m/s), tempurature (K) and pressure (Pa) [Vfeed Tfeed PH PI PL]
yF = 0.7; % Mole fraction for methane
%%%%%%% Preliminary calculations %%%%%%%%%%%
A = pi*r^2; % Bed cross-sectional area m2
V = A*L; % Bed volume m^3
mBed = A*L*426.7; % Mass of adsorbent (in bed) 426.7 = bed density
%%%%%%%%% loop around here to work out p and n and mole fraction for each
%%%%%%%%% loop
% y is an array size n*4 of y1 = 1:n, q1 = n+1:2*n, y2 = 2*n+1:3*n,
% q2 =3*n+1:4*n T = 4*n+1:5*n
sol = adsorptionSolver(t,z,yF,TPF);
% sol.x is time steps
% sol.y is
format long
purityh = 100*sol.y(3*n,end) /(sol.y(n,end) + sol.y(3*n,end));
fprintf('Purity of hydrogen is %f%% after %4.2fs\n', purityh, tf)
% purityh = sol.y(3*n,:) ./(sol.y(n,:) + sol.y(3*n,:));
%purityh = (sol.y(n,:) + sol.y(3*n,:)); % mole fraction
% figure
% plot(sol.x, purityh)
% xlabel('time')
% ylabel('purity of hydrogen')
figure(1)
subplot(1,2,1)
mesh(sol.x,z,sol.y(1:n,1:end))
xlabel('t')
ylabel('bed length')
zlabel('mole fraction y1')
title('mole fraction of methane')
subplot(1,2,2)
mesh(sol.x,z,sol.y(2*n+1:3*n,1:end))
xlabel('t')
ylabel('bed length')
zlabel('mole fraction y2')
title('mole fraction of hydrogen')
figure(2)
subplot(1,2,1)
mesh(sol.x,z,sol.y(n+1:2*n,1:end))
xlabel('t')
ylabel('bed length')
zlabel('loading q1 mol/kg')
title('loading of methane')
subplot(1,2,2)
mesh(sol.x,z,sol.y(3*n+1:4*n,1:end))
xlabel('t')
ylabel('bed length')
zlabel('loading q2 mol/kg')
title('loading of hydrogen')
figure(3)
subplot(1,2,1)
mesh(sol.x,z,sol.y(4*n+1:5*n,1:end))
xlabel('t')
ylabel('bed length')
zlabel('Temp')
title('temp of system')
subplot(1,2,2)
mesh(sol.x,z,sol.y(5*n+1:6*n,1:end))
xlabel('time')
ylabel('bed length')
zlabel('Pressure')
title('Pressure of system')
%% Solver function
function out = adsorptionSolver(t,z,yFeed,TPFeed)
n = numel(z); % Size of mesh grid
h = diff(z);
h = h(1);
% start up condidtions
Tw = 300; % Ambient Tempurature K
Pw = 5.95e5; % Ambient Pressure Pa
y1w = 0.7;
% Initial Conditions
y1_0 = zeros(n,1) + y1w;
q1_0 = zeros(n,1);
q2_0 = zeros(n,1);
T_0 = zeros(n,1) + Tw;
P_0 = zeros(n,1) + Pw;
% y is an array size n*4 of c1 = 1:n, q1 = n+1:2*n, c2 = 2*n+1:3*n, q2 = 3*n+1:4*n
y0 = [y1_0; q1_0; q2_0; T_0; P_0];
% Solving
out = ode23s(@(t,y) adsorptionModel(t,y,h,n,yFeed,TPFeed),t,y0);
end
%% Adsorption model
function dydt = adsorptionModel(t,y,h,n,yiFeed,TPFeed)
% Variables allocation
% y is an array size n*5 of y1 = 1:n, q1 = n+1:2*n,
% q2 = 2*n+1:3*n, T = 3*n+1:4*n, P = 4*n+1:5*n
y1 = y(1:n);
y2 = 1 - y1;
q1 = y(n+1:2*n);
q2 = y(2*n+1:3*n);
T = y(3*n+1:4*n);
P = y(4*n+1:5*n);
% Z half points
yhalf = zeros(n+1,1);
Thalf = zeros(n+1,1);
Phalf = zeros(n+1,1);
vhalf = zeros(n+1,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Constant physical properties and basic simulation parameters
R = 8.3145; % ideal gas constant J/mol/K
D = 1.3e-5; % Axial Dispersion coefficient m2/s
epl = 0.4043; % Void fraction of bed
eplp = 0.546; % Void fraction of particle
eplt = epl + (1-epl)*eplp; % Total void fraction
k = [0.136;0.259]; % (lumped) Mass Transfer Coefficient 1/s
rhop = 716.3; % Particle density kg/m3
rhob = 426.7; % Bed density kg/m3 = (1-epl)*rhop
Cps = 1046.7; % heat capacity of solid J/kg/K
Kl = 1.2e-6; % Thermal diffusity J/m/s/K
Kz = 0.09; % Effective gas thermal conductivity
deltaH = [24124; 8420]; % heat of adsorption J/mol
Vfeed = TPFeed(1);
Tfeed = TPFeed(2);
PH = TPFeed(3);
PI = TPFeed(4);
PL = TPFeed(5);
lambda = 0.5; % rate of pressure change (1/s)
% Ergun equation parameters
visc = 3.73e-8; % gas viscosity kg/m/s
Rp = 5.41e-3; % particle radius m
% Constants for CH4
c1a = -0.703029; c1b = 108.4773; c1c = -42.52157; c1d = 5.862788; c1e = 0.678565;
% Constants for H2
c2a = 33.066178; c2b = -11.363417; c2c = 11.432816; c2d = -2.772874; c2e = -0.158558;
% Cpg = heat capacity (J/mol*K)
Cpg1 = c1a + c1b*(T/1e3) + c1c*(T/1e3).^2 + c1d*(T/1e3).^3 + c1e./((T/1e3).^2);
Cpg2 = c2a + c2b*(T/1e3) + c2c*(T/1e3).^2 + c2d*(T/1e3).^3 + c2e./((T/1e3).^2);
Cpg = (y1.*Cpg1 + y2.*Cpg2);
%%%%%%%% Boundary conditions: j = 0.5 and n + 0.5 %%%%%%%%%%%%%%%%%%%%%%%%
%%%% Pressurization
Phalf(1) = PH - (PH - PL)*exp(-lambda*t);
Phalf(n+1) = P(n);
vhalf(1) = -(2/visc/h)*(4/150*(epl/(1-epl))^2)*(Rp^2)*(P(1)-Phalf(1));
vhalf(n+1) = -(2/visc/h)*(4/150*(epl/(1-epl))^2)*(Rp^2)*(Phalf(n+1)-P(n));
yhalf(1) = (y1(1)+yiFeed*vhalf(1)*h) / (2*D+vhalf(1)*h);
yhalf(n+1) = y1(n);
% Inlet gas condtions
rhogInlet = (Phalf(1) .* (yhalf(1)*16.04e-3 + (1-yhalf(1))*2.02e-3)) ./ (R*Tfeed);
Cpg1in = c1a + c1b*(Tfeed/1e3) + c1c*(Tfeed/1e3).^2 + c1d*(Tfeed/1e3).^3 + c1e./((Tfeed/1e3).^2);
Cpg2in = c2a + c2b*(Tfeed/1e3) + c2c*(Tfeed/1e3).^2 + c2d*(Tfeed/1e3).^3 + c2e./((Tfeed/1e3).^2);
CpgInlet = (yhalf(1).*Cpg1in + (1-yhalf(1)).*Cpg2in);
%
Thalf(1) = (T(1) + Tfeed*vhalf(1)*epl*rhogInlet*CpgInlet*h)/(2*Kz + vhalf(1)*epl*rhogInlet*CpgInlet*h);
Thalf(n+1) = T(n);
%%%%%%%% Velocity wall values: j 1.5 to n - 0.5 %%%%%%%%%%%%%%%%%%%%%%%%
vhalf(2:n) = -(2/visc/h)*(4/150*(epl/(1-epl))^2)*(Rp^2)*(P(2:n) - P(1:n-1));
%%%%%%%% Wall values: j = 1.5 %%%%%%%%%%%%%%%%%%%%%%%%
alpha0P = (2/3) / (P(2)-P(1) + 1e-10)^4;
alpha1P = (1/3) / (2*(P(1)-Phalf(1)) + 1e-10)^4;
Phalf(2) = alpha0P/(alpha0P+alpha1P)*(0.5*(P(1)+P(2))) + alpha1P/(alpha0P+alpha1P)*(2*P(1)-Phalf(1));
alpha0y = (2/3) / (y1(2)-y1(1) + 1e-10)^4;
alpha1y = (1/3) / (2*(y1(1)-yhalf(1)) + 1e-10)^4;
yhalf(2) = alpha0y/(alpha0y+alpha1y)*(0.5*(y1(1)+y1(2))) + alpha1y/(alpha0y+alpha1y)*(2*y1(1)-yhalf(1));
alpha0T = (2/3) / (T(2)-T(1) + 1e-10)^4;
alpha1T = (1/3) / (2*(T(1)-Thalf(1)) + 1e-10)^4;
Thalf(2) = alpha0T/(alpha0T+alpha1T)*(0.5*(T(1)+T(2))) + alpha1T/(alpha0T+alpha1T)*(2*T(1)-Thalf(1));
%%%%%%%% Wall values: j = 2.5 to n - 0.5 %%%%%%%%%%%%%%%%%%%%%%%%
alpha0P = (2/3) ./ (P(3:n)-P(2:n-1) + 1e-10).^4;
alpha1P = (1/3) ./ (P(2:n-1)-P(1:n-2) + 1e-10).^4;
Phalf(3:n) = alpha0P./(alpha0P+alpha1P).*(0.5.*(P(2:n-1)+P(3:n))) + alpha1P./(alpha0P+alpha1P).*(1.5.*P(2:n-1)- 0.5* P(1:n-2));
alpha0y = (2/3) ./ (y1(3:n)-y1(2:n-1) + 1e-10).^4;
alpha1y = (1/3) ./ (y1(2:n-1)-y(1:n-2) + 1e-10).^4;
yhalf(3:n) = alpha0y/(alpha0y+alpha1y)*(0.5*(y1(2:n-1)+y1(3:n))) + alpha1y./(alpha0y+alpha1y).*(1.5.*y1(2:n-1)- 0.5* y1(1:n-2));
alpha0T = (2/3) ./ (T(3:n)-T(2:n-1) + 1e-10).^4;
alpha1T = (1/3) ./ (T(2:n-1)-T(1:n-2)+ 1e-10).^4;
Thalf(3:n) = alpha0T/(alpha0T+alpha1T)*(0.5*(T(2:n-1)+T(3:n))) + alpha1T./(alpha0T+alpha1T).*(1.5.*T(2:n-1)- 0.5* T(1:n-2));
rhog = (P .* (y1*16.04e-3 + y2*2.02e-3)) ./ (R*T); % ρ = PM/RT, M: H2 = 2.02 g/mol CH4 = 16.04 g/mol Air= 28.97 g/mol 50/50 H2 and CH4 = 9.03
%%%%%%%%%%%%%%% Langmuir Equation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % CH4
a11 = -9.5592; a21 = 4638.5; b11 = 3.725e-4/100000; b21 = 1.443e4;
% % H2
a12 = -23.131; a22 = 8069.1; b12 = 2.248/100000; b22 = -1.435e4;
qs1 = a11 + (a21./T);
qs2 = a12 + (a22./T);
B1 = b11.*exp(b21./(8.3145*T));
B2 = b12.*exp(b22./(8.3145*T));
qstar1 = (qs1.*B1.*P.*y1) ./ (1+((P.*y1.*B1)+(P.*(1-y1).*B2))); %mols/kg
qstar2 = (qs2.*B2.*P.*(1-y1)) ./ (1+((P.*y1.*B1)+(P.*(1-y1).*B2)));
LT1 = qstar1-q1;
LT2 = qstar2-q2;
%%%%%%%%%%%%%%% Eqautions for loading of component i %%%%%%%%%%%%%%
dq1dt = k(1)*LT1;
dq2dt = k(2)*LT2;
%%%%%%%%%%%%%% Balance equations
dPdt = zeros(n,1);
dy1dt = zeros(n,1);
dTdt = zeros(n,1);
% TMB
dPdt(1:n) = (T(1:n)./h).*(vhalf(2:n+1).*Phalf(2:n+1)./Thalf(2:n+1)-vhalf(1:n).*Phalf(1:n)./Thalf(1:n)) - ((R*T(1:n)*(1-epl))/epl).*(dq1dt(1:n)+dq2dt(1:n)) + P(1:n)/T(1:n) * dTdt(1:n);
% CMB
dy1dt(1) = (D.*T(1)./(h.*P(1))) * ( (Phalf(2)./Thalf(2)).*((y(2)-y(1))/h)-(Phalf(1)./Thalf(1)).*(2/h*(y(1)-yhalf(1)))) - (T(1)./P(1)/h).*(vhalf(2).*Phalf(2).*yhalf(2)./Thalf(2) - vhalf(1).*Phalf(1).*yhalf(1)./Thalf(1)) - ((R*T(1)*(1-epl))/epl./P(1)).*dq1dt(1) - (y1(1)/P(1))*dPdt(1) + (y1(1)/T(1))*dTdt(1);
dy1dt(2:n-1) = (D.*T(2:n-1)./(h.*P(2:n-1))) .* ( (Phalf(3:n)./Thalf(3:n)).*((y(3:n)-y(2:n-1))/h)-(Phalf(2:n-1)./Thalf(2:n-1)).*((y(2:n-1)-y(1:n-2))/h)) - (T(2:n-1)./P(2:n-1)/h).*(vhalf(3:n).*Phalf(3:n).*yhalf(3:n)./Thalf(3:n) - vhalf(2:n-1).*Phalf(2:n-1).*yhalf(2:n-1)./Thalf(2:n-1)) - ((R.*T(2:n-1).*(1-epl))./epl./P(2:n-1)).*dq1dt(2:n-1) - (y1(2:n-1)./P(2:n-1)).*dPdt(2:n-1) + (y1(2:n-1)./T(2:n-1)).*dTdt(2:n-1);
dy1dt(n) = (D.*T(n)./(h.*P(n))) * ( (Phalf(n+1)./Thalf(n+1)).*(2/h*(yhalf(n+1)-y(n)))-(Phalf(n)./Thalf(n)).*((y(n)-y(n-1))/h)) - (T(n)./P(n)/h).*(vhalf(n+1).*Phalf(n+1).*yhalf(n+1)./Thalf(n+1) - vhalf(n).*Phalf(n).*yhalf(n)./Thalf(n)) - ((R*T(n)*(1-epl))/epl./P(n)).*dq1dt(n) - (y1(n)/P(n))*dPdt(n) + (y1(n)/T(n))*dTdt(n);
% EB
Cpa = (q1.*Cpg1)./(q1+q2+1e-10) + (q2.*Cpg2)./(q1+q2+1e-10); % 1e-10 is there to stop nan
ebcoeff = 1./( ((1-epl)/epl).*(rhop.*Cps + Cpa.*(dq1dt+dq2dt)) );
dTdt(1) = ebcoeff(1)*(Kl/epl/h*((T(2)-T(1))/h-2*(T(1)-Thalf(1))/h)-Cpg(1)./(R*h).*(Phalf(2).*vhalf(2)-Phalf(1).*vhalf(1))-((1-epl)/epl).*Cpa(1).*T(1).*(dq1dt(1)+dq2dt(1)) + ((1-epl)/epl)*(-deltaH(1)*dq1dt(1)-deltaH(2)*dq2dt(1)) - Cpg(1)/R.*dPdt(1));
dTdt(2:n-1) = ebcoeff(2:n-1).*(Kl/epl/h*((T(3:n)-2.*T(2:n-1)+T(1:n-2))./h)-Cpg(2:n-1)./(R*h).*(Phalf(3:n).*vhalf(3:n)-Phalf(2:n-1).*vhalf(2:n-1))-((1-epl)/epl).*Cpa(2:n-1).*T(2:n-1).*(dq1dt(2:n-1)+dq2dt(2:n-1)) + ((1-epl)/epl).*(-deltaH(1).*dq1dt(2:n-1)-deltaH(2).*dq2dt(2:n-1)) - Cpg(2:n-1)./R.*dPdt(2:n-1));
dTdt(n) = ebcoeff(n)*(Kl/epl/h*(2*(Thalf(n+1)-T(n))/h-(T(n)+T(n-1))/h)-Cpg(n)./(R*h).*(Phalf(n+1).*vhalf(n+1)-Phalf(n).*vhalf(n))-((1-epl)/epl).*Cpa(n).*T(n).*(dq1dt(n)+dq2dt(n)) + ((1-epl)/epl)*(-deltaH(1)*dq1dt(n)-deltaH(2)*dq2dt(n)) - Cpg(n)/R.*dPdt(n));
%%%%%%%%%%%%%%%% Concatenate vectors of time derivatives
dydt = [dy1dt;dq1dt;dq2dt;dTdt;dPdt];
end

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

采纳的回答

William Rose
William Rose 2024-1-12
编辑：William Rose 2024-1-12
@Thomas Stone-Wigg,
Does your code run? Do you want help making it run at all, or help making it run better? If better, better how? You said speed is not an issue. If it doesn't run, what is the error message?
I notice that you use t in multiple roles on the right hand side below. That may not turn out well.
t = t0:dt:tf; % Time vector
% Solving
out = ode23s(@(t,y) adsorptionModel(t,y,h,n,yFeed,TPFeed),t,y0);
I would do the following instead
tspan = t0:dt:tf; % Time vector
% Solving
out = ode23s(@(t,y) adsorptionModel(t,y,h,n,yFeed,TPFeed),tspan,y0);
Good luck.
  2 个评论
Thomas Stone-Wigg
Thomas Stone-Wigg 2024-1-12
Hi William, thank you for your reply.
My code was not functioning at the time of posting, as I initially believed the issue was solely with the solver. I have since fixed this by resolving the coupled differential equations as simultaneous equations. This allowed me to pass the individual differential equations to the solver:
for i=1:n
aX = (T(i)./h).*(vhalf(i+1).*Phalf(i+1)./Thalf(i+1)-vhalf(i).*Phalf(i)./Thalf(i)) - ((R*T(i)*(1-epl))/epl).*(dq1dt(i)+dq2dt(i));
bX = P(i)/T(i) ;
Cpa = (q1.*Cpg1)./(q1+q2+1e-10) + (q2.*Cpg2)./(q1+q2+1e-10); % 1e-10 is there to stop nan
ebcoeff = 1./( ((1-epl)/epl).*(rhop.*Cps + Cpa.*(dq1dt+dq2dt)) );
% CMB and EB
if i==1
cX = (D.*T(1)./(h.*P(1))) * ( (Phalf(2)./Thalf(2)).*((y(2)-y(1))/h)-(Phalf(1)./Thalf(1)).*(2/h*(y(1)-yhalf(1)))) - (T(1)./P(1)/h).*(vhalf(2).*Phalf(2).*yhalf(2)./Thalf(2) - vhalf(1).*Phalf(1).*yhalf(1)./Thalf(1)) - ((R*T(1)*(1-epl))/epl./P(1)).*dq1dt(1) ;
dX = (D.*T(1)./(h.*P(1))) * (y1(1)/P(1));
eX = (D.*T(1)./(h.*P(1))) * (y1(1)/T(1));
fX = ebcoeff(1)*(Kl/epl/h*((T(2)-T(1))/h-2*(T(1)-Thalf(1))/h)-Cpg(1)./(R*h).*(Phalf(2).*vhalf(2)-Phalf(1).*vhalf(1))-((1-epl)/epl).*Cpa(1).*T(1).*(dq1dt(1)+dq2dt(1)) + ((1-epl)/epl)*(-deltaH(1)*dq1dt(1)-deltaH(2)*dq2dt(1)) );
gX = Cpg(1)/R;
elseif i==n
cX = (D.*T(i)./(h.*P(i))) * ( (Phalf(i+1)./Thalf(i+1)).*(2/h*(yhalf(i+1)-y(i)))-(Phalf(i)./Thalf(i)).*((y(i)-y(i-1))/h)) - (T(i)./P(i)/h).*(vhalf(i+1).*Phalf(i+1).*yhalf(i+1)./Thalf(i+1) - vhalf(i).*Phalf(i).*yhalf(i)./Thalf(i)) - ((R*T(i)*(1-epl))/epl./P(i)).*dq1dt(i) ;
dX = (D.*T(i)./(h.*P(i))) * (y1(i)/P(i));
eX = (D.*T(i)./(h.*P(i))) * (y1(i)/T(i));
fX = ebcoeff(i)*(Kl/epl/h*(2*(Thalf(i+1)-T(i))/h-(T(i)+T(i-1))/h)-Cpg(i)./(R*h).*(Phalf(i+1).*vhalf(i+1)-Phalf(i).*vhalf(i))-((1-epl)/epl).*Cpa(i).*T(i).*(dq1dt(i)+dq2dt(i)) + ((1-epl)/epl)*(-deltaH(1)*dq1dt(i)-deltaH(2)*dq2dt(i)) );
gX = Cpg(i)/R;
elseif (i > 1) && (i < n)
cX = (D.*T(i)./(h.*P(i))) .* ( (Phalf(i+1)./Thalf(i+1)).*((y(i+1)-y(i))/h)-(Phalf(i)./Thalf(i)).*((y(i)-y(i-1))/h)) - (T(i)./P(i)/h).*(vhalf(i+1).*Phalf(i+1).*yhalf(i+1)./Thalf(i+1) - vhalf(i).*Phalf(i).*yhalf(i)./Thalf(i)) - ((R.*T(i).*(1-epl))./epl./P(i)).*dq1dt(i) ;
dX = (D.*T(i)./(h.*P(i))) .* (y1(i)./P(i));
eX = (D.*T(i)./(h.*P(i))) .* (y1(i)./T(i));
fX = ebcoeff(i).*(Kl/epl/h*((T(i+1)-2.*T(i)+T(i-1))./h)-Cpg(i)./(R*h).*(Phalf(i+1).*vhalf(i+1)-Phalf(i).*vhalf(i))-((1-epl)/epl).*Cpa(i).*T(i).*(dq1dt(i)+dq2dt(i)) + ((1-epl)/epl).*(-deltaH(1).*dq1dt(i)-deltaH(2).*dq2dt(i)) );
gX = Cpg(i)./R;
else
end
% Numerical computation of derivatives
DPdt = aX / (1 - bX * gX); % Revised from symbolic solve
Dy1dt = (cX - dX * DPdt) / (1 - eX * DPdt); % Revised
DTdt = (fX - gX * DPdt); % Revised
dPdt(i) = DPdt;
dy1dt(i) = Dy1dt;
dTdt(i) = DTdt;
end
Regarding speed, i have a lot of repeating equations and im sure there is a more efficient way to code them.
I have managed to get the code to run but mole fraction increases rapidly when i run it with ode15s and a mass matrix. The run stops with an error message:
Error using daeic12 (line
168)
Need a better guess y0 for
consistent initial
conditions.
Error in ode15s (line 304)
[y,yp,f0,dfdy,nFE,nPD,Jfac]
=
daeic12(odeFcn,odeArgs,t,ICtype,Mt,y,yp0,f0,...
Error in
Rogue2>adsorptionSolver
(line 132)
out = ode15s(@(t,y)
adsorptionModel(t,y,h,n,yFeed,TPFeed),t,y0,opts);
Error in Rogue2 (line 33)
sol =
adsorptionSolver(t,z,yF,TPF);
% % Creating mass matrix
M = eye(5*n); % Mass matrix for temporal domain integrator
M(1,1) = 0; % Algebraic equation for yi1 at z = 0;
M(n,n) = 0; % Algebraic equation for yi1 at z = L
M(2*n+1,2*n+1) = 0; % Algebraic equation for yi2 at z = 0;
M(3*n,3*n) = 0; % Algebraic equation for yi2 at z = L
M(4*n+1,4*n+1) = 0;
opts = odeset('Mass',M,'MassSingular','yes'); % creates mass matrix for the function ode15s
% Solving
out = ode15s(@(t,y) adsorptionModel(t,y,h,n,yFeed,TPFeed),t,y0,opts);
The code runs indefinitly with the ode23s in the original post.
Theres a very high possibility that i have made a mistake in the balance equations/boundary equations/wall value equations also, ill check the method and them over again.
Thanks again

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Particle & Nuclear Physics 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by