How to solve an ODE with an array using ODE45?

I have an ODE, udY/dx=S , where Y and S are arrays , how do I solve it using ODE45 ???

[Jan]

What is "u" in "udY"?

[DIP]

velocity. its constant

[Jan]

@DIP: Please provide more details. "velocity" does not help, because Matlab does not care about the physical meaning of variables. Post the function you want to integrate, preferrably in valid Matlab code. If it is written correctly, all you have to do is calling ode45.

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 Jan, Star Strider,

Answer 2

[Jan]

ode45 works with vectors, therefor if Y is a vector, the provided examples in the docs should help: https://www.mathworks.com/help/matlab/ref/ode45.html#examples . If not, please post more details by editing the question.

[EDITED, Considering your comments] Perhaps you mean this:

function main
% Initial Conditions
Cco(1)  = 0.52; %Y1
Co2(1)  = 0.99; %Y2
Ch2o(1) = 0;    %Y3

L     = 0.12;
N     = 39;
delx  = L/(N-1);
x     = 0:delx:0.12;
C0    = [Cco; Co2; Ch2o];
[x,C] = ode45(@DiffEQ, x, C0);

plot (x,C)
end

function dC = DiffEQ(x, C)
u = 1.5;

A  = 2.2e10;
Ea = 130;
Ru = 8.3145;
T  = 750;

Cco  = C(1);
Co2  = C(2);
Ch2o = C(3);

Rco  = -A * (exp(-Ea / (Ru*T))) * Cco * sqrt(Co2); %S1
Ro2  = 0.5 * Rco;  %S2
Rh2o = -Rco;       %S3
dC   = [Rco; Ro2; Rh2o] / u;
end

This needs a lot of processing time. Using ode15s was successful in less than a second (are you sure the system is not stiff). The diagram looks nicer, if you define:

x = [0, 0.12];

[DIP]

Hi Jan, I'll repost the whole question in detail.

I have an ODE, u*dY/dx=S

Y is an array [Y1 Y2 Y3] of concentrations of different species in a chemical reaction. S is an array for the homogeneous kinetic expressions. [S1 S2 S3]. S1, S2, S3 themselves are exponential equations whose values are known. I have to solve this using ODE45. u is a constant. Is it a straightforward solution with ode45?

[Jan]

If the system is non-stiff, this should be a straightforward solution with ODE45. Try it and post what you have done so far, if you get any troubles.

[DIP]

Hi Jan, the system is non stiff. The values [Y1 Y2 Y3] are different and the values [S1 S2 S3] are different too. How do i formulate the ODE 45 syntax for this ?

% Constants

u=1.5;
L=0.12;
N=39;

A=2.2*10e10;
Ea=130;
Ru=8.3145;
T=750;

% Initial Conditions
Cco(1)=0.52; %Y1
Co2(1)=0.99; %Y2
Ch2o(1)=0;   %Y3

Rco=-A*(exp((-Ea)/(Ru*T)))*Cco*sqrt(Co2); %S1
Ro2=0.5*Rco; %S2
Rh2o=-Rco; %S3

x(1)=0; %Initial Condition

delx=L/(N-1);
x = 0:delx:0.12;
C=[Cco Co2 Ch2o];
R=[Rco Ro2 Rh2o];
C0=[Cco(1) Co2(1) Ch2o(1)];
[x,C] = ode45(@(x,C) R/u, x(1), C0);

plot (x,C)

[DIP]

Hi Jan, Im supposed to use the Ode45 solver only.

[DIP]

The ode 15s does not give the accurate plot.

[DIP]

Hi Jan , do you have any documentation for the last line of your main function ?? What does it mean ?

[DIP]

Hi Jan, I am extremely thankful for all the help.

I modified your code, im able to speed up the calculation using ode45, but I cant get it to plot like a curve, it gives me a straight line.

% 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 Co2 Ch2o];
Cco  = C(1);
Co2  = C(2);
Ch2o = C(3);
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,'-o')
title('Molar Concentration vs. Distance');
legend('CO','CO2','H2O');
xlabel('Axial(x) Direction [m]');
ylabel('Concentrations [mol/m3]');

[DIP]

If i get rid of the indexing for the constants and initial conditions, I get the same plot. Why am I not able to get the correct plot mentioned a few comments above ?

% 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   = 0.52; 
Co2   = 0.99; 
Ch2o  = 0; 
delx  = L/(N-1);
C0    = [Cco Co2 Ch2o];
Rco   = -A * (exp(-Ea / (Ru*T))) * Cco * sqrt(Co2); 
Ro2   = 0.5 * Rco;  
Rh2o  = -Rco;  

x     = [0:delx:0.12];
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,'-o')
title('Molar Concentration vs. Distance');
legend('CO','O2','CO2');
xlabel('Axial(x) Direction [m]');
ylabel('Concentrations [mol/m3]');

[DIP]

I think its got to do with defining the arrays of C and R correctly.

[Jan]

@DIP: @(x,C) R/u is constant. Changes in the coordinates are not considered, because all calls of "R/u" reply the same value. Better use a normal function instead of a anonymous function. I've showed you already, how to do this.

Concerning:

 documentation for the last line of your main function

do you mean "plot(x,C)"?

Answer 3

[Star Strider]

See if the approach in Monod kinetics and curve fitting will do what you want.

You will have to adapt it to your data and differential equations. The essential design should work.