How to solve a ODE that is a function of a cubic polynomial f(x) = 0?

1 visualizzazione (ultimi 30 giorni)
Many greetings to anyone who stumbles upon this thread. I'm currently trying to solve this differential equation for a chemical reaction rate regression analysis.
Anotación 2019-12-04 132452.png
As you can see, I have 15 parameters that I've already fitted and real & known values. The differential equation has DCF as y and O3 and FeOOH are also known constant values. My only problem is writing a code that continously finds the roots of [.OH] to input inside the differential equation. My time span is 20 minutes and my initial value of DCF is 32.439. I've tried with the following code:
function dy = ozecat(t,y,x,g,ox)
OH = [x(4) (x(5)*y+x(6)*g(1)+x(7)*g(1)/y-x(8)*ox(1)) (x(9)*ox(1)*y+x(10)*(ox(1)^2)-x(11)*ox(1)*g(1)+x(12)*(ox(1)*g(1)/y)) (-x(13)*ox(1)^2*y-x(14)*ox(1)^2*g(1)+x(15)*(ox(1)^2*g(1)/y))];
dy = -x(1)*ox(1)*y-x(2)*y*roots(OH)-x(3)*g(1)*ox(1);
end
where g stands for FeOOH data, ox for O3 data and x are my array of 16 parameters (from a to l and the k's). I've also tried to select one single root of OH by adding this two lines of code to the above:
rroot_OH = roots(OH);
rroot_OH = rroot_OH(imag(rroot_OH) == 0;
But to no avail. I get the following errors. This first for the original code:
Error using horzcat
Dimensions of matrices being concatenated are not consistent.
Error in ozecat (line 2)
OH = [x(4) (x(5)*y+x(6)*g(1)+x(7)*g(1)/y-x(8)*ox(1)) (x(9)*ox(1)*y+x(10)*(ox(1)^2)-x(11)*ox(1)*g(1)+x(12)*(ox(1)*g(1)/y))
(-x(13)*ox(1)^2*y-x(14)*ox(1)^2*g(1)+x(15)*(ox(1)^2*g(1)/y))];
Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.
and this one for the modified with the two lines added:
Error using odearguments (line 95)
OZECAT returns a vector of length 3, but the length of initial conditions vector is 1. The vector returned by OZECAT and the initial
conditions vector must have the same number of elements.
Error in ode45 (line 115)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
Error in parametrows (line 34)
[t_calc y] = ode45(@ozecat, tspan, y0,[],x,g,ox);
So how exactly is this done? Or is it not possible by ode45?
  3 Commenti
Miguel Figueroa
Miguel Figueroa il 4 Dic 2019
Yes, I'm aware. Thanks for the input. However, if I try to give each root a place inside the differential equation by means of making the initial condition [y0 y*0 y**0] it still gives me a horzcat error. ¿what can I do in that case?

Accedi per commentare.

Risposte (1)

Steven Lord
Steven Lord il 4 Dic 2019
This doesn't look like a system of ODEs. The presence of the variable d in the first set of curly braces threw me for a second, but I think this is a system of one differential equation and one algebraic equation. That makes it a system of differential algebraic equations (DAEs.) You may want to use the technique shown for the Robertson problem on that page as a model for your DAE function.
If you do this, you may want to choose names for your variables inside your DAE function carefully so they look closer to your mathematical equations. This may help you if you need to debug or present your work to someone more familiar with chemical reaction rate regression analysis than I am.
function dDCF = ozecat(t,DCF,params,FeOOH,O3)
a = params(1); % params used to be called x
b = params(2);
% etc
% Not sure how OH is calculated
OH = ...;
term(1) = a*OH.^3;
term(2) = (OH.^2).* ...;
% ...
dDCF = ...;
end
  1 Commento
Miguel Figueroa
Miguel Figueroa il 6 Dic 2019
Hello, sorry for the late reply. I've been trying to follow your advice and use ode15s solving it as a DAE system. However, I get an error of 'Index exceeds matrix dimensions' when I restructure my function to look just like the Robertson problem you mentioned. I guess it might have something to do with the polynomial giving 3 roots to the Differential equation below?
Here's how it looks now:
y(1) is the x of a cubic polynomial f(x)
y(2) is the solution to my ODE dy/dt. y(1) is inside the ODE.
function dy = ozecat(t,y,k_o,k_oh,k_9,ka,kb,kc,kd,ke,kf,kg,kh,ki,kj,kk,kl,feooh,gotehita_1_set)
dy = [ka*y(1)^3+(y(1)^2)*(kb*y(2)+kc*gotehita_1_set+kd*gotehita_1_set/y(2)-ke*feooh)...
+y(1)*(kf*feooh*y(2)+kg*(feooh^2)-kh*feooh*gotehita_1_set+ki*(feooh*gotehita_1_set/y(2)))...
+(-kj*feooh^2*y(2)-kk*feooh^2*gotehita_1_set+kl*(feooh^2*gotehita_1_set/y(2)));
-k_o*feooh*y(2)-k_oh*y(2)*y(1)-k_9*gotehita_1_set*feooh];
end

Accedi per commentare.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by