MATLAB: Solving System of ODE that includes Algebraic Eqn’s

dsolveMATLAB C/C++ Math Libraryode45

[Deleted]

Best Answer

function main
  y0 = [2; 1; 0; 400];
  tspan = [0 10];
  Fao = 3.3;
  Cao = 0.01983;
  K = 1.6;
  [t,y] = ode15s(@(t,y)fun(t,y,Fao,Cao,K),tspan,y0);
  plot(t,y(:,1),t,y(:,2),t,y(:,3))
end
function dy = fun(t,y,Fao,Cao,K)
  Fa = y(1);
  Fb = y(2);
  Fc = y(3);
  T =  y(4);
  Xa = (Fao - Fa)/Fao;
  Ke = 0.0885*exp(4000/T);
  Ca = (Cao*(1-Xa))/(1+Xa/3);
  Cb = (Cao*(2*Xa))/(1+Xa/3);
  ra = K*(Ca-Cb^2/Ke);
  rb = 2*ra;
  rc = -ra;
  
  dy = [ra;rb;rc;(ra*200)/(Fa*4+Fb*2+Fc*4)];
end

Related Solutions

MATLAB: Extract values from ODE45

‘How Do I extract the values of Ra and Rb and plot them ?’

You don’t ‘extract’ them. Just recalculate them from the integrated values.

Make a few changes to your ‘Example2’ function:

function [] = Example2()
% clc 
% clear 
% close all
y0 = [100.; 0; 0]; % Initial values for the dependent variables 
tspan= linspace(0,20,50);
options = optimset('display','off');
[t,x]= ode45(@ODEfun,tspan,y0,options);
figure
plot(t,x)
legend('Fa','Fb','Fc')
Fa = x(:,1); 
Fb = x(:,2); 
Fc = x(:,3);
Kc = 0.01; 
Ft = Fa + Fb + Fc; 
Co = 1; 
K = 10; 
Kb = 40; 
ra = 0 - (K .* Co / Ft .* (Fa - (Co .^ 2 * Fb .* Fc .^ 2 ./ (Kc * Ft .^ 2)))); 
Ka = 1; 
Ra = Ka * Co * Fa / Ft;
Rb = Kb * Co * Fb / Ft;
figure
plot(t, [Ra Rb])
legend('R_a', 'R_b')
end

This will produce the plot you want. You could plot them in the same figure as the first plot (using the hold function), however you will likely not be able to see them well. I would use separate plots.

MATLAB: Undefined function or variable using ODE45

I believe this is what you actually intend:

function dy = fun(T,Y)
FT0 = 42;
FA0=10; %mol/min
W=100; %KG
a=.008; %kg^-1
k1C=2; %L/kg cat/min


K1C=.2; %mol/L
k2D=.4; %L/kg cat/min
k3E=5; %L/mol^2/kg cat/min
kc=1; %L/kg cat/min
v=16.67; %L/min
Fa=Y(1);
Fb=Y(2);
Fc=Y(3);
Fd=Y(4);
Fe=Y(5);
y=Y(6);
CT0=FA0;
FT=Fa+Fb+Fc+Fd+Fe;
Ca=CT0*(Fa/FT)*y;
Cb=CT0*(Fb/FT)*y;
Cc=CT0*(Fc/FT)*y;
Cd=CT0*(Fd/FT)*y;
Ce=CT0*(Fe/FT)*y;
r1C=k1C*(Ca-(Cb*Cc/K1C));
r2D=k2D*Cd;
r3E=k3E*Cd*(Ce^2);
ra=-r1C-r2D;
rb=r1C;
rc=r1C-r3E;
rd=r2D-.5*r3E;
re=r3E;
dFadw=ra;
dFbdw=rb;
dFcdw=rc;
dFddw=rd;
dFedw=re;
dydw=-a*FT/(2*FT0*y);
dy = [dFadw;dFbdw;dFcdw;dFddw;dFedw;dydw];
end
wspan=0:1:100;
x0=[10;0;0;0;0;1];
[T,X] = ode45(@fun,wspan,x0);
figure
plot(T,X)
grid

The ‘FT0’ value and assignment are missing, so I created one to test the code. (You may want to substitute a different value.) This runs without error.

Best Answer

Related Solutions

MATLAB: Extract values from ODE45

MATLAB: Undefined function or variable using ODE45

Related Question