Concentration of CO in Los Angeles Basin (matlab)

[gfd43tg]

Homework Statement

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*10^{12} \frac{ft^3}{hr}$. 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 $2*10^{-10} \frac{lb mol}{ft^3}$.

Given Information:

F_{CO,A} + F_{CO,S} - Q*C_{CO} = V\frac{dC_{CO}}{dt}
F_{CO,A} = a + b sin(\pi*t/6)
Where $a = 20254815 \frac {lb mol}{hr}$ and $b = 0.5*a \frac {lb mol}{hr}$
V=4*10^{13} ft^3
$Q=1.67*10^{13} \frac{ft^3}{hr}$ for $t \lt 72$ and $Q= 1.67*10^{12} \frac{ft^3}{hr}$ for $t \ge 72$
F_{CO,S}=C_{CO,S}*Q
C_{CO,S} = 2*10^{-10} \frac{lb mol}{ft^3}
Initial condition: C_{CO,t=0} = 2*10^{-10} \frac{lb mol}{ft^3}

Homework Equations

The crux of my problem is using ode45 to solve this differential equation
$F_{CO,A} + F_{CO,S} - Q*C_{CO} = V\frac{dC_{CO}}{dt}$

The Attempt at a Solution

Hello, I want to give a little context for this problem so that you can understand what is going on. $F_{CO,S}$ is the flow rate of carbon monoxide entering the basin with a concentration of carbon monoxide $C_{CO,s}$ from the Santa Ana winds. $F_{CO,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\frac{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 \ge 72$, and I'm not sure how to remedy this in the ode45 built in function. It is that $Q$ is a piecewise function of $t$, as well as $C_{CO} = f(t)$. Anyways, here is what I wrote so far. I know the function handle is supposed to have @(t,x) for a derivative $\frac {dx}{dt}$, so I used in my function handle $@(t,C_{CO})$ for $\frac {dC_{CO}}{dt}$. I keep getting the error Undefined function or variable 't'.

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');

 

[DrClaude]

What you should do is define your function using "function", such that you can incorporate the code defining Q

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

 

[gfd43tg]

I'm confused how to make it work how you are describing

function dC_COdt = fHan (t, C_CO)

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

   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

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

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

Will not work, I get this error

fHan
Error: File: fHan.m Line: 21 Column: 1
This statement is not inside any function.
 (It follows the END that terminates the definition of the
 function "fHan".)

 

[DrClaude]

The file fHan.m should contain only the function, starting from the "function" statement and terminating with the "end".

This function will then be called by your main program. You still need to initialize the time and set the initial conditions. Note also that the function is called "fHan", not "dC_COdt". You simply need to remove some lines from your original posting, leaving:

tSpan = [0 100]; % hr

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

[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');

 

[gfd43tg]

Something isn't right still, I just copy paste your code in #2

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 make a separate script, and I put in arguments for fHan in my ode45 call,

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 the 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, ...

the ode45 function doesn't seem to like the way I input the first argument

 

[DrClaude]

Something isn't right still, I just copy paste your code in #2
Then make a separate script, and I put in arguments for fHan in my ode45 call,

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');

Do not put arguments, it's ode45 which supplies the correct arguments.

What was missing in my post is to indicate that fHan is a function:

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

 

[gfd43tg]

Thanks, ode45 arguments behave somewhat mysteriously to me. It just magically knows what to put it

 

[DrClaude]

Thanks, ode45 arguments behave somewhat mysteriously to me. It just magically knows what to put it

Of course, there is nothing magical about it

It is basically a stepper method: starting from some initial conditions $y_0$ at $x_0$, it uses a discrete form of the differential equation to find the solution $y_1$ at $x_1$. For that, it needs to know the differential equations $dy/dx$, which is what you supply. In the simplest case, you would have
$$
y_1 = y_0 + (x_1 - x_0) \left. \frac{dy}{dx} \right|_{x_0}
$$
which the algorithm can calculate by feeding $x_0$ and $y_0$ into your function.

In reality, the algorithm used by ode45, the Runge-Kutta-Fehlberg 4-5 method, is a more sophisticated version of the equation above with a smaller error. In addition, it can make an estimate of the error on $y_1$, and it reduces it by subdividing the interval between $x_0$ and $x_1$ into smaller intervals, until a desired accuracy is obtained. For that, it will need to call your function many times to get $dy/dx$ at different values of $x$. In the end, ode45 returns the values of $x$ for all the intermediate steps it took, along with the corresponding value of $y(x)$.

A good introduction to the subject can be found in the book Numerical Recipes in C, see http://apps.nrbook.com/c/index.html, starting at page 707.