Trouble with ODE45 for an array of values

Hi,

I have an equation udC/dx=R, where C and R are arrays. I've written a code but its giving a linear plot. Need help

% Solve ODE for Multi-Species Isothermally
% Constants
u = 1.5;
A  = 2.2e10;
Ea = 130;
Ru = 8.3145;
T  = 750;
L  = 0.12;
N  = 39;

% Initial Conditions

Cco(1)   = 0.52; 
Co2(1)  = 0.99; 
Ch2o(1)  = 0; 
delx  = L/(N-1);
x     = 0:delx:0.12;
C0    = [Cco(1) Co2(1) Ch2o(1)];
C     = zeros(3,1);

Rco   = -A * (exp(-Ea / (Ru*T))) * Cco * sqrt(Co2); 
Ro2   = 0.5 * Rco;  
Rh2o  = -Rco;       
dC    = [Rco Ro2 Rh2o] / u;
R     = [Rco Ro2 Rh2o];
R     = zeros(3,1);

% ODE45 solver
[x,C] = ode45(@(x,C) R/u, x, C0);

plot(x,C,'-+')
title('Molar Concentration vs. Distance');
legend('CO','O2','CO2');
xlabel('Axial(x) Direction [m]');
ylabel('Concentrations [mol/m3]');

Thank you.

Answer 1

[DIP]

% Solution 4: Chemical Species Concentration as a Function of Axial Distance

clc
clear

% Constants
L  = 0.12;
N  = 39;
T(1:N) = 750;
delx     = L/(N-1);
xspan    = 0:delx:0.12;
C0       = [0.52 0.99 0.0];

% ODE 45 Code
[x,C]=ode45('hw2_4',xspan,C0);

% Plots
subplot(2,1,1);
plot(x,C(:,1),'b--o',x,C(:,2),'g--+',x,C(:,3),'r--s')
legend('C_{CO}','C_{O2}','C_{CO2}')
xlabel('Axial (x) Direction [m]')
ylabel('Concentrations [mol/m^3]')
title('Molar Concentration vs Distance')

Function File:

function dcdx=hw2_4(x,C)

% Constants
u = 1.5;
A  = 2.2e10;
Ea = 130000;
Ru = 8.3145;
T  = 750;

% Rate Constants

R(1)     = -A * (exp(-Ea / (Ru*T))) * C(1) * sqrt(C(2)); 
R(2)     = 0.5 * R(1);  
R(3)     = -R(1);
dcdx     = [R(1);R(2);R(3)]/u;

[DIP]

Thank You Torsten, Walter.

Answer 2

[Walter Roberson]

You have

[x,C] = ode45(@(x,C) R/u, x, C0);

and R and u are constants, 3 x 1 and 1 x 1. The integral of a constant is the straight line a*x + b so all of your lines are linear, not curved.

Note that you have

R     = [Rco Ro2 Rh2o];
R     = zeros(3,1);

which constructs R and then overwrites it with 0. You might have wanted

R     = [Rco Ro2 Rh2o] .';

which would give you lines with different slopes, not curves.

You cannot get a curve unless your function depends upon the time or the input derivatives.

[DIP]

How can I correct my code ?

i get a plot like this which is still wrong.

[Walter Roberson]

Your line

[x,C] = ode45(@(x,C) R/u, x, C0);

makes the derivative values independent of x or C. Your function needs to depend upon at least one of the two of those in order to get a curve.

When I glance over your code, my eye is drawn to

Rco   = -A * (exp(-Ea / (Ru*T))) * Cco * sqrt(Co2);

which gives me the impression that Rco should depend upon time, since T often means time. But possibly in your case T means Temperature; if so then I would wonder whether maybe temperature should depend upon time?

It would help if you posted the equations so we could look to see whether you correctly went from equation to code.

[DIP]

Reposting the question,

[DIP]

even I have the feeling that I have made an error in defining the arrays for R.

[Walter Roberson]

Notice in the line after the one that defines u, the numeric values for the C_CO and so on are given with "i=1". Those are initial conditions.

Notice in the definition of R_CO that the C_CO and C_O2 do not have the "i=1". They are not constants: they are functions that could also be written as C_CO(x) and C_O2(x) and then those relate back to the initial u*dC/dx = R where C is a vector with [C_CO, C_O2, C_CO2] and C so is really the derivative of [C_CO, C_O2, C_CO2] with respect to x.

If you go back and write in the (x) everywhere appropriate then it becomes much more clear where the differential equations are.

[DIP]

Walter, do I need to use a loop for that ?

[Torsten]

[x,C] = ode45(@(x,C) [-A * (exp(-Ea / (Ru*T))) * C(1) * sqrt(C(2)) ; 0.5*(-A * (exp(-Ea / (Ru*T))) * C(1) * sqrt(C(2))) ; -(-A * (exp(-Ea / (Ru*T))) * C(1) * sqrt(C(2)))]/u, x, C0);

Best wishes

Torsten.

[Walter Roberson]

Thanks, Torsten.

[DIP]

Hi Walter, i get an error saying that x is not an appropriate index, when i index the R calculations. When I modify the code as per Torsten's suggestions, it solves forever.

[DIP]

are you able to run the code and get any answers with the modifications ?

[Torsten]

Use a stiff solver, e.g. ODE15S, instead of ODE45.

Best wishes

Torsten.