Differential model won't calculate with powers <1
조회 수: 1 (최근 30일)
이전 댓글 표시
Hi everyone,
I am working on a model which describes multiple reactions over time, using simple power law models (R = k*[C]^n).
ode1 = diff(Ca) == -k1*Ca;
cond1 = Ca(0) == 5;
ode2 = diff(Cb) == k1*Ca - k2*Cb;
cond2 = Cb(0) == 0;
ode3 = diff(Cc) == k2*Cb - k3*Cc;
cond3 = Cc(0) == 0;
ode4 = diff(Cd) == k3*Cc - k4*Cd;
cond4 = Cd(0) == 0;
ode5 = diff(Ce) == k4*Cd;
cond5 = Ce(0) == 0;
t works just like I want to with n = 1, however, our data suggests that n < 1. I tried adding powers to my concentrations, but then, Matlab has a hard time calculating it, and it never finishes.
I want to calculate the concentrations of all components over time. All constants (k1, k2, k3, k4) and reaction orders will be added later. However, the model should be able to calculate the concentrations when the order is lower than 1.
Do you have any suggestions on how to tackle this problem? Do I have to use other functions?
댓글 수: 2
채택된 답변
Bjorn Gustavsson
2021년 3월 8일
It is not given that a general nonlinear system of ODEs have analytical solutions. If you cannot get the symbolic toolbox to find one for you, the best (possibly?) advice is to turn to numerical solutions. For that have a look at ode45 and its siblings, and the large number of ODE-examples. This should be very straightforward to implement for numerical integration. The ODE-function might look something like this:
function dCdt = your_power_reactions_ODE(t,C,k,g)
Ca = C(1);
Cb = C(2);
Cc = C(3);
Cd = C(4);
Ce = C(5);
dCadt = -k(1)*Ca^g(1);
dCbdt = k(1)*Ca^g(1) - k(2)*Cb^g(2);
dCcdt = k(2)*Cb^g(2) - k(3)*Cc^g(3);
dCddt = k(3)*Cc^g(3) - k(4)*Cd^g(4);
dCedt = k(4)*Cd^g(4);
% etc...
end
That should be pretty standard to integrate. It might be a numerical stiff set of ODEs - I've not spent any time to judge that. You then simply set the initial conditions and the time-period and then run:
C0 = [5 0 0 0 0];
t_span = [0 17];
k_all = [1 log(2),2^pi,12];
g_all = 1./[1:4];
[t_out,C_out] = ode45(@(t,C) your_power_reactions_ODE(t,C,k_all,g_all),...
t_span,C0);
If it is too stiff, try ode15s, ode113, or ode23s...
HTH
댓글 수: 3
Bjorn Gustavsson
2021년 3월 8일
But to my inderstanding these are all "decay-type reactions" - like for the concentrations of elements in a radioactive decay-chain?
추가 답변 (0개)
참고 항목
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!