ODE45 problem with concentration of CO

1 次查看（过去 30 天）
Rick
Rick 2014-9-10
评论： Star Strider 2014-9-11
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?

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

采纳的回答

Star Strider
Star Strider 2014-9-10
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.
  3 个评论
显示 1更早的评论
Rick
Rick 2014-9-11
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
Star Strider 2014-9-11
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);

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Lengths and Angles 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by