How to solve an ODE with an array using ODE45?

8 次查看（过去 30 天）

显示更早的评论

DIP 2017-2-20

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/325911-how-to-solve-an-ode-with-an-array-using-ode45

评论： DIP 2017-3-6

采纳的回答： DIP

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

3 个评论
显示 1更早的评论隐藏 1更早的评论

DIP 2017-2-20

编辑：DIP 2017-2-20

velocity. its constant

Jan 2017-2-20

@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.

请先登录，再进行评论。

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

采纳的回答

DIP 2017-3-6

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/325911-how-to-solve-an-ode-with-an-array-using-ode45#answer_257476

在 MATLAB Online 中打开

% 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;

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

DIP 2017-3-6

Thank you Jan, Star Strider,

请先登录，再进行评论。

更多回答（2 个）

Jan 2017-2-20

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/325911-how-to-solve-an-ode-with-an-array-using-ode45#answer_255511

编辑：Jan 2017-2-21

在 MATLAB Online 中打开

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];

10 个评论
显示 8更早的评论隐藏 8更早的评论

DIP 2017-2-21

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

Jan 2017-2-21

在 MATLAB Online 中打开

@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)"?

请先登录，再进行评论。