Trouble with ODE45 for an array of values

9 views (last 30 days)
DIP on 27 Feb 2017
Edited: DIP on 8 Apr 2017
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.

DIP on 6 Mar 2017
% 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 on 6 Mar 2017
Thank You Torsten, Walter.

Walter Roberson on 28 Feb 2017
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.
Torsten on 1 Mar 2017
Use a stiff solver, e.g. ODE15S, instead of ODE45.
Best wishes
Torsten.