ODE45 problem with concentration of CO

Smog begins to build up again immediately after a Santa Ana wind passes through the basin. The volumetric flow rate through the basin has dropped to 1.67 x 10^12 ft3/hr (1.5 mph). Plot the concentration of carbon monoxide in the basin as a function of time for up to 72 hours and just after a Santa Ana wind. The initial concentration of CO is 0.08 ppm (2 x 10^-10 lbmol/ft3).

Given F_CO,A + F_CO,S - Q*C_CO = V(dC_CO/dt)

F_CO,A = a + b sin(π*t/6)

Where a = 35,000 lb mol/h and b = 30,000 lb mol/h

V=4x10^13 ft3

Q =1.67x10^13 ft3/hr for t <72, Q= 1.67x10^12 ft3/hr for t ≥72

F_CO,S= C_CO,S*Q

C_CO,S = 2x10^-10 lb mol/ft3

Initial condition: C_CO(t=0)= 2x10^-10 lb mol/ft3

Hello, I want to give a little context for this problem so that you can understand what is going on. FCO,S is the flow rate of carbon monoxide entering the basin with a concentration of carbon monoxide (CO) CCO,s. FCO,A is the flow rate of carbon monoxide created by the cars driving inside the basin. Q*C_CO is the flow rate of carbon monoxide exiting the basin, and C_CO is changing with time because this is an unsteady state process. My problem is essentially solving the differential equation given in the problem

F_CO,A + F_CO,S - Q*C_CO = V(dC_CO/dt)

There are a couple of problems for me. First, there is the fact that Q changes depending on if t < 72 or if t >= 72, and I'm not sure how to remedy this in the ode45 built in function. Anyways, here is what I wrote so far. Keep in mind my a and b are different from the problem statement, but it should have no effect on where my issue lies here.

tSpan = [0 100]; % hr
if t < 72
    Q = 1.67e13; % ft^3/hr
else
    Q = 1.67e12; % ft^3/hr
end
C_COs = 2e-10; % lb mol/ft^3
F_COs = C_COs*Q; % lb mol/hr

a = 20254815; % lb mol/hr
b = 0.5*a; % lb mol/hr

F_COa = a + b*sin(pi*t/6); % lb mol/hr
Volume = 4e13; % ft^3

fHan = @(t,C_CO) (F_COa + F_COs - C_CO*Q)/Volume;
[T, C_CO] = ode45(fHan,tSpan,C_COs);
plot(T,C_CO)
xlabel('time (hours)');
ylabel('CO Concentration (lb mol/ft^3)');
title('CO Concentration vs. time');

and I get an error Undefined function or variable 't'. Well, I remember with ode45 you need to have your function handle be @(t,x). Also, how do I make it so that it knows with my if/else conditions how to handle the changing volumetric flow rate, Q?

Answer 1

[Star Strider]

What you want to do is a bit challenging with an anonymous function for your ODE function, but it can be done.

You need to create three functions of t:

Q = @(t) [1.67e13*(t<72) + 1.67e12*(t>=72)];
F_COs = @(t) C_COs*Q(t); % lb mol/hr
F_COa = @(t) a + b*sin(pi*t/6); % lb mol/hr

The Q function uses a combination of numerical and logical statements to force it to take on the appropriate values at various values of t. If (t < 72), that logical value equates to 1, so Q takes equates to 1.67e13. Since (t >= 72) is by definition false at the same value of t, that term evaluates to zero. The reverse occurs when (t >= 72). It’s a one-line if-then-else with only two conditions.

The other two lines simply create your previous statements as functions of t so your ODE function can call them with the values of t it gets from ode45.

Your ODE then becomes:

fHan = @(t,C_CO) (F_COa(t) + F_COs(t) - C_CO*Q(t))/Volume;

Your resulting code (with some statements commented-out) is:

tSpan = [0 100]; % hr
% if t < 72
%     Q = 1.67e13; % ft^3/hr
% else
%     Q = 1.67e12; % ft^3/hr
% end
Q = @(t) [1.67e13*(t<72) + 1.67e12*(t>=72)];
C_COs = 2e-10; % lb mol/ft^3
F_COs = @(t) C_COs*Q(t); % lb mol/hr

a = 20254815; % lb mol/hr
b = 0.5*a; % lb mol/hr

F_COa = @(t) a + b*sin(pi*t/6); % lb mol/hr
Volume = 4e13; % ft^3

% fHan = @(t,C_CO) (F_COa + F_COs - C_CO*Q)/Volume;
fHan = @(t,C_CO) (F_COa(t) + F_COs(t) - C_CO*Q(t))/Volume;
[T, C_CO] = ode45(fHan,tSpan,C_COs);
plot(T,C_CO)
xlabel('time (hours)');
ylabel('CO Concentration (lb mol/ft  ^3)');
title('CO Concentration vs. time');

It runs, and seems to give the correct result, with the appropriate jump starting at (t = 72). Only you know if it does.

[Rick]

Thanks, out of curiosity, how might I do this if I were just to write a function instead of using an anonymous function?

[Rick]

Here is my function

function dC_COdt = fHan (t, C_CO)

     if t < 72
         Q = 1.67e13; % ft^3/hr
     else
         Q = 1.67e12; % ft^3/hr
     end

     C_COs = 2e-10; % lb mol/ft^3
     F_COs = C_COs*Q; % lb mol/hr

     a = 20254815; % lb mol/hr
     b = 0.5*a; % lb mol/hr

     F_COa = a + b*sin(pi*t/6); % lb mol/hr
     Volume = 4e13; % ft^3

     dC_COdt = (F_COa + F_COs - C_CO*Q)/Volume;

end

then in a separate script I put in

tSpan = [0 100]; % hr

C_COs = 2e-10; % lb mol/ft^3

[T, C_CO] = ode45(fHan(0,C_COs),tSpan,C_COs);
plot(T,C_CO)
xlabel('time (hours)');
ylabel('CO Concentration (lb mol/ft^3)');
title('CO Concentration vs. time');

but I get this error

Error using exist
The first input to exist must be a string.

Error in odearguments (line 60)
  if (exist(ode)==2)

Error in ode45 (line 114)
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs,
odeFcn, ...

[Star Strider]

I didn’t run your newest code, but with function files, you need to refer to them as function handles in the code that calls them. If they aren’t already declared as functions (with function handles) in your code, you have to declare them as such with a leading ‘@’ sign.

This should work:

[T, C_CO] = ode45(@fHan,tSpan,C_COs);

It’s not necessary to include the argument list since you aren’t passing other arguments to your ODE function ‘fHan’, and in the example you posted, t would always be 0 anyway. That’s not what you want.

If you were going to pass other arguments (in this illustration a vector of parameters called ‘parms’) to your ODE function using an anonymous function construction, you would do it this way:

[t,y] = ode45(@(t,y) odefun(t,y,parms),tspan,y0);